# 教材位置：第2章，第2.9节；仅终端输出。
# 单位、变量与形状：points 为 (100000,2)，inside 为 (100000,)；坐标、圆周率估计和标准误
# 均无量纲。
# 学习目标：通过几何概率估计 π，并理解蒙特卡洛估计的标准误和近似置信区间。
# 核心关系：在单位正方形均匀抽点，落入第一象限单位圆内的概率 q₀=π/4；
# 若 I 是是否落入圆内的指示变量，则 π估计=4×样本均值，理论标准误为 4sqrt(q₀(1-q₀)/n)。
# 数据与输出：抽取 100000 个二维点，输出 π 的估计、估计标准误、近似 95% 区间和实际绝对误
# 差。
# 运行方式：安装 NumPy 后执行 python e05_monte_carlo.py；不保存外部数据或图片。

import numpy as np

# 固定种子用于复现实验；n 是独立随机点的个数，而不是每个坐标轴的网格点数。
rng = np.random.default_rng(3)
n = 100000
# random((n, 2)) 生成形状为 (100000, 2) 的数组，每行是一个 (x, y)。
# 两个坐标各自均匀分布于 [0,1)，因此这些点均匀覆盖单位正方形。
points = rng.random((n, 2))
# 沿 axis=1 对每行的 x²、y² 求和，得到到原点距离的平方，无需求平方根。
# <=1 判断点是否在四分之一单位圆内，inside 是形状 (n,) 的布尔数组。
# 连续均匀模型下边界概率为 0，是否包含圆弧不改变理论面积比。
inside = np.sum(points ** 2, axis=1) <= 1
# mean 将 True 当 1、False 当 0，得到落入圆内的样本比例；q 是对 q₀ 的估计。
q = inside.mean()
# 四分之一圆的面积乘以 4 得到整圆面积，因此将样本比例乘以 4 来估计 π。
pi_hat = 4 * q
# 由独立伯努利样本的方差得到比例标准误，再利用常数倍的标准差缩放规则乘以 4。
# 用 q 代替未知 q₀ 得到代入估计；它描述重复抽样时估计量的波动，不等于本次实际误差。
se_hat = 4 * np.sqrt(q * (1 - q) / n)
print("圆周率估计值：", pi_hat)
print("估计标准误：", se_hat)
# 大样本下利用正态近似构造 Wald 区间，1.96 是标准正态约 97.5% 分位数。
# [-1,1] 生成左右两个端点；约 95% 指重复实验中这种区间的覆盖率，而非必然包含真值。
# 当命中比例非常接近 0 或 1、或样本过少时，这种近似可能不好；本例命中比例适中且 n 大。
print("近似 95% 置信区间：", pi_hat + np.array([-1, 1]) * 1.96 * se_hat)
# π 在本例中已知，所以能直接计算实现后的绝对误差；一般蒙特卡洛问题未必知道真值。
print("绝对误差：", abs(pi_hat - np.pi))
# 确保每个二维点对应一个判断结果，防止求和方向错误。
assert inside.shape == (n,)
# 将 q₀=π/4 代入理论标准误并化简，检查两个等价公式；这里未要求模拟估计恰好等于 π。
assert np.isclose(4 * np.sqrt((np.pi / 4) * (1 - np.pi / 4) / n),
                  np.sqrt(np.pi * (4 - np.pi) / n))