Quantify Monte Carlo uncertainty
Estimate a circle area with seeded random points, use Wilson intervals and inspect repeated-experiment coverage.
Goal and how to run
This lesson assumes basic Python knowledge. Read each method and caption before changing parameters. AI exercises are optional; no AI account is needed for the core workflow. Data lives in notebook metadata, so there is no companion data download. Keep the notebook saved in the active working folder.
Your goal is to understand the research question, execute transparent calculations, check the results, and state the limits of the conclusion.
Setup and parameters
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
Load the embedded snapshot
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']
Read the method functions
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
Data and model boundary
Analytic functions and Monte Carlo points are mathematical teaching models, not experimental measurements. When Iris observations are used, row exchangeability and the sampling design remain assumptions; convenience records alone do not establish population representativeness. Numerical discretization error is distinct from sampling uncertainty.
Steps
1. Sample a known geometric probability
Uniform points in the unit square fall inside the quarter-unit circle with probability pi/4. Multiplying the fraction by four estimates pi. The points are explicitly simulated; seed controls reproducibility, not scientific representativeness.
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. Follow a single experiment
At checkpoints, calculate the binomial Wilson interval before scaling its endpoints by four. Point estimates can fluctuate as sample size grows; convergence does not require every successive error to decrease.
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. Inspect repeated-experiment coverage
Repeat 40 independent simulated runs at n=2000. The displayed intervals illustrate finite-simulation coverage: around 95% is a long-run target, not a guarantee that exactly 38 of 40 will contain pi.
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. Compare with an independent implementation and export
Use SciPy’s binomial test object to obtain a Wilson proportion interval and compare before scaling. Distinguish probability, point estimate, interval confidence level and observed simulation coverage.
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
Exercises and reference explanation
- Why multiply the proportion interval by four?
- Does a 95% confidence level guarantee coverage in a particular run?
Reference explanation. The probability target is pi/4. Long-run confidence does not guarantee that one interval or exactly 95% of a small simulation batch covers the parameter.
Keep a copy before changing parameters. Compare old and new figures and record which inputs, calculations and conclusions changed.
Optional AI coding exercises
Prompt 1
Using NumPy/SciPy/Matplotlib, compare n=500, 2000 and 8000 across independent seeded simulations.
Prompt 2
With the same packages, compare Wilson with the normal Wald interval, including success counts 0 and n.
Human acceptance. Check probability endpoints before scaling, reproduce the SciPy Wilson reference and report all simulation seeds.
Review the edited source, then manually run affected cells and subsequent checks. Ask the assistant to explain the purpose, packages and data provenance of its changes, and reconcile that explanation with actual output.
Optional online extension
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
Sources, snapshots and next steps
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.
Next steps. Transfer these methods to your own public or authorized data while retaining input checks, recorded parameters and result validation. Document the new experimental design and method limits; this example does not validate a new experiment. Full data licenses and generation instructions accompany the cookbook.