量化 Monte Carlo 不确定性

用固定随机点估计圆面积，计算 Wilson 区间，并检查重复实验的覆盖率。

数学与统计 · 第 04 课🔋 离线核心分析⏱ 50 分钟教学模拟数据📦 numpy · pandas · scipy · matplotlib · IPython

目标与运行方法

本课适合具备 Python 基础的本科生。先阅读每一步的方法与图注，再修改参数。AI 练习为可选部分，不需要 AI 账号即可完成核心课程。数据写在 notebook metadata 中，不需要另下载数据文件。请保持 notebook 已保存到当前工作文件夹。

学习任务是理解本课的研究问题、执行透明的计算、校验结果，并说明结论的边界。

环境与参数

In [1] · Python

from pathlib import Path
import base64, gzip, hashlib, json, platform
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from IPython.display import display, HTML
import math
from scipy.stats import binomtest
from IPython.display import display as _display
def display(value):
    if isinstance(value,pd.DataFrame):
        table=value.to_html(max_rows=12,escape=True)
        style="<style>.vibeit-readable-table{max-width:100%;overflow-x:auto}.vibeit-readable-table table{width:max-content;max-width:none;border-collapse:collapse}.vibeit-readable-table td,.vibeit-readable-table th{white-space:nowrap!important;word-break:normal!important;overflow-wrap:normal!important}</style>"
        return _display(HTML(style+'<div class="vibeit-readable-table">'+table+'</div>'))
    return _display(value)
RECIPE_ID='ma04-monte-carlo'
OUTPUT_DIR=Path.cwd()/"results"/RECIPE_ID
plt.rcParams.update({"figure.dpi":120,"font.size":10,"axes.spines.top":False,"axes.spines.right":False,
                      "axes.prop_cycle":plt.cycler(color=["#18785f","#416aa6","#bc6541","#9974a5"])})
print("Python:",platform.python_version(),"; NumPy:",np.__version__,"; Pandas:",pd.__version__)
print("Working folder:",Path.cwd())
print("Core data mode: embedded snapshot; no network requests")

Python: 3.13.13 ; NumPy: 2.5.3 ; Pandas: 3.0.6
Working folder: /Users/zhiluo/Developer/vibeit-cookbook/cookbook/downloads
Core data mode: embedded snapshot; no network requests

读取内嵌快照

In [2] · Python

def load_snapshot(recipe_id, expected_sha256, snapshot_key):
    """Read embedded data from this notebook, including a renamed copy."""
    for path in sorted(Path.cwd().glob("*.ipynb")):
        try:
            notebook = json.loads(path.read_text(encoding="utf-8"))
            metadata = notebook.get("metadata", {}).get("vibeit_cookbook", {})
            blob = metadata.get("snapshots", {}).get(snapshot_key)
            if metadata.get("recipe_id") != recipe_id or not blob:
                continue
            if blob.get("sha256") != expected_sha256:
                continue
            raw = gzip.decompress(base64.b64decode(blob["gzip_base64"], validate=True))
            if hashlib.sha256(raw).hexdigest() != expected_sha256:
                raise ValueError("Embedded snapshot checksum mismatch: " + snapshot_key)
            return json.loads(raw)
        except (OSError, json.JSONDecodeError):
            continue
    raise FileNotFoundError(
        "Save/open the complete .ipynb in the active working folder. "
        "The embedded snapshot could not be located for " + recipe_id)

SNAPSHOT_HASHES={'models': '48fd2e6e706650fe2ee0c0ad80b6d794a546a1dec78dd9d370ac09db3cc3adbf'}
def snapshot(key):
    return load_snapshot(RECIPE_ID,SNAPSHOT_HASHES[key],key)
print("Embedded snapshots:",list(SNAPSHOT_HASHES))

Embedded snapshots: ['models']

方法函数

In [3] · Python

def wilson_interval(successes,trials,z=1.959963984540054):
    """A binomial Wilson score interval; the returned endpoints are probabilities."""
    if not isinstance(trials,(int,np.integer)) or trials<=0 or not 0<=successes<=trials or z<=0:
        raise ValueError('Invalid success count, trial count or z')
    proportion=successes/trials;denominator=1+z*z/trials
    center=(proportion+z*z/(2*trials))/denominator
    half=z*np.sqrt(proportion*(1-proportion)/trials+z*z/(4*trials*trials))/denominator
    return float(max(0,center-half)),float(min(1,center+half))

print("Method ready: wilson_interval")

Method ready: wilson_interval

数据与模型边界

解析函数与 Monte Carlo 点属于数学教学模型，不是实验测量。使用 Iris 观测时，可交换性和采样设计仍是统计假设；方便样本不自动代表总体。数值离散误差与采样不确定性应分别解释。

分步分析

1. 采样已知几何概率

单位正方形的均匀点落在四分之一单位圆内的概率为 pi/4，比例乘四可估计 pi。这些点明确属于模拟，随机种子控制复现，不代表科学采样代表性。

In [4] · Python

seed=int(snapshot("models")["seed"]);rng=np.random.default_rng(seed)
TRIALS=20000;points=rng.uniform(0,1,size=(TRIALS,2));inside=np.sum(points**2,axis=1)<=1
estimate=4*inside.mean();interval=4*np.asarray(wilson_interval(int(inside.sum()),TRIALS))
print("Estimated pi:",estimate,"; Wilson-derived 95% interval:",interval)
part=points[:1500];mask=inside[:1500]
fig,ax=plt.subplots(figsize=(5,5));ax.scatter(*part[mask].T,s=4,label="Inside");ax.scatter(*part[~mask].T,s=4,label="Outside")
theta=np.linspace(0,np.pi/2,101);ax.plot(np.cos(theta),np.sin(theta),color="black")
ax.set(xlabel="Simulated x",ylabel="Simulated y",xlim=(0,1),ylim=(0,1),aspect="equal",title="Simulated quarter-circle experiment");ax.legend();plt.tight_layout();plt.show()

Estimated pi: 3.146 ; Wilson-derived 95% interval: [3.12306454 3.16849531]

2. 追踪一次实验

在检查点先计算二项 Wilson 区间，再将端点乘四。估计随样本量增加仍会波动；收敛不要求每一步误差都下降。

In [5] · Python

sizes=np.array([100,200,500,1000,2000,5000,10000,20000]);cumulative=np.cumsum(inside)
estimates=4*cumulative[sizes-1]/sizes
intervals=np.array([4*np.asarray(wilson_interval(int(cumulative[n-1]),int(n))) for n in sizes])
fig,ax=plt.subplots(figsize=(8,4));ax.plot(sizes,estimates,"o-");ax.fill_between(sizes,*intervals.T,alpha=.2)
ax.axhline(math.pi,color="black",linestyle="--",label="Known pi")
ax.set(xscale="log",xlabel="Simulated trial count",ylabel="Pi estimate",title="Single-run estimates and 95% Wilson intervals");ax.legend();plt.tight_layout();plt.show()

3. 检查重复实验的覆盖率

在 n=2000 下重复 40 次独立模拟。图中的区间展示有限模拟覆盖率：约 95% 是长期目标，不保证本次恰有 38 个区间包含 pi。

In [6] · Python

repeated=[]
for run in range(40):
 sample=rng.uniform(0,1,size=(2000,2));success=int(np.sum(np.sum(sample**2,axis=1)<=1))
 lo,hi=4*np.asarray(wilson_interval(success,2000));repeated.append([run,4*success/2000,lo,hi,lo<=math.pi<=hi])
repeat=pd.DataFrame(repeated,columns=["run","estimate","lower","upper","contains_pi"])
fig,ax=plt.subplots(figsize=(8,5))
for row in repeat.itertuples():ax.plot([row.lower,row.upper],[row.run,row.run],color="green" if row.contains_pi else "orange",lw=1)
ax.axvline(math.pi,color="black",linestyle="--");ax.set(xlabel="Pi interval",ylabel="Simulated run",title="Repeated Wilson intervals");plt.tight_layout();plt.show()
print("Observed repeated-run coverage:",float(repeat.contains_pi.mean()))

Observed repeated-run coverage: 0.975

4. 对照独立实现并导出

使用 SciPy 二项检验对象取得 Wilson 比例区间，先核对再缩放。区分概率、点估计、置信水平和本次模拟覆盖率。

In [7] · Python

reference=binomtest(int(inside.sum()),TRIALS).proportion_ci(confidence_level=.95,method="wilson")
assert np.allclose(interval,4*np.array([reference.low,reference.high]))
assert np.all((intervals>=0)&(intervals<=4))
OUTPUT_DIR.mkdir(parents=True,exist_ok=True)
pd.DataFrame({"n":sizes,"estimate":estimates,"lower":intervals[:,0],"upper":intervals[:,1]}).to_csv(OUTPUT_DIR/"monte-carlo-checkpoints.csv",index=False)
repeat.to_csv(OUTPUT_DIR/"repeated-intervals.csv",index=False)
(OUTPUT_DIR/"simulation.json").write_text(json.dumps({"seed":seed,"trials":TRIALS,"known_probability":math.pi/4,"known_pi":math.pi,"interval":"binomial Wilson, scaled by four"},indent=2))
print("Independent interval check passed; exports:",OUTPUT_DIR)

Independent interval check passed; exports: /Users/zhiluo/Developer/vibeit-cookbook/cookbook/downloads/results/ma04-monte-carlo

自主练习与参考答案

为何把比例区间乘四？
95% 置信水平是否保证本次覆盖？

参考解释。 概率目标为 pi/4。长期置信水平不能保证某一区间或一小批模拟恰有 95% 覆盖参数。

先保留当前 notebook 副本，再修改参数。把新图与原图一起比较，并记录改变了哪些输入、哪些计算和哪些结论。

可选的 AI 编程练习

提示词 1

仅用 NumPy、SciPy、Matplotlib，在独立固定种子模拟中比较 n=500、2000、8000。

提示词 2

使用同样的包比较 Wilson 与正态 Wald 区间，包含成功次数为零或 n 的情况。

人工验收。 缩放前检查概率端点，复现 SciPy Wilson 参考，并记录所有种子。

检查修改后的源代码，再手动运行受影响单元及后续校验。要求助手解释修改的目的、使用的包和数据来源，并用实际输出核对解释。

可选联网扩展

In [8] · Python

RUN_ONLINE_UPDATE=False
if RUN_ONLINE_UPDATE:
    import requests
    try:
        response=requests.get('https://docs.scipy.org/doc/scipy/tutorial/integrate.html',timeout=20)
        response.raise_for_status()
        print("Source:",response.url)
        print(response.text[:700])
    except requests.RequestException as error:
        print("Online extension unavailable:",type(error).__name__,str(error)[:160])
else:
    print("Online extension skipped; embedded core data unchanged")

Online extension skipped; embedded core data unchanged

来源、快照与下一步

IAPWS saturation reference release — IAPWS SR1-86(1992): unrestricted publication allowed in all countries (cover page); retrieved 2026-10-09. SHA-256 a10d38d73339b0713c86cfac88f9adcba137202167ccdfa1cd58238c438361d7. Reference equations and numerical coefficients transcribed from release pages 2–3; pedagogical model values are not new experiments.
Authored numerical teaching model parameters — Original teaching parameters/code MIT; reference equation coefficients attributed to IAPWS SR1-86(1992); retrieved 2026-10-09. SHA-256 4350c4305a83c772fa295a8be0ef86a247f2bfaae930a081720d3cd0c753abe0. Declared parameters and seed; generated trajectories and calibration observations must be labelled simulated.

下一步。 在保留本课的数据检查、参数记录和结果校验之后，将相同方法用于自己的公开或已获授权数据。记录研究设计与方法限制，不把示例结果当成新实验的验证。完整数据许可与生成说明见 cookbook 网页。