审计模拟吸光度校准

生成明确声明的 Beer–Lambert 教学观测，拟合校准、检查残差，并向未知量传播不确定性。

化学 · 第 06 课🔋 离线核心分析⏱ 60 分钟教学模拟数据📦 numpy · pandas · 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
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='ch06-calibration'
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 linear_calibration(concentration,absorbance):
    """OLS with an intercept; concentrations must state their units in the calling cell."""
    x=np.asarray(concentration,dtype=float);y=np.asarray(absorbance,dtype=float)
    if x.ndim!=1 or x.shape!=y.shape or len(x)<3 or not np.isfinite(x+y).all() or np.ptp(x)==0:
        raise ValueError('Need finite, varying calibration concentrations and >2 observations')
    design=np.column_stack([np.ones(len(x)),x]);coefficient=np.linalg.lstsq(design,y,rcond=None)[0]
    residual=y-design@coefficient;variance=float(residual@residual/(len(x)-2))
    covariance=variance*np.linalg.inv(design.T@design)
    return dict(coefficients=coefficient,residual=residual,variance=variance,covariance=covariance,fitted=design@coefficient)

print("Method ready: linear_calibration")

Method ready: linear_calibration

数据与模型边界

PubChem 性质与构象是数据库计算字段，不是新实验测量。本课 V2000 解析范围明确限定为原子、来源键型与形式电荷。滴定及校准使用已声明的理想/模拟教学假设；水蒸气压使用署名的 IAPWS 参考关联式。描述符、距离和视觉邻近不能证明生物活性或结合强度。

分步分析

1. 生成明确标记的模拟校准观测

模型为 A=blank+epsilon·l·c，再加声明的高斯噪声。九个浓度水平各有三次模拟重复。这是方法基准，不是实验室数据。浓度显示为 μmol/L，吸光度无量纲。

In [4] · Python

models=snapshot("models");p=models["calibration"];rng=np.random.default_rng(models["seed"])
concentration=np.repeat(np.linspace(0,80,9),3)
truth=p["blank_absorbance"]+p["absorptivity_L_per_mol_cm"]*p["path_length_cm"]*concentration*1e-6
absorbance=truth+rng.normal(0,p["noise_sd_absorbance"],len(truth))
fitted=linear_calibration(concentration,absorbance);intercept,slope=fitted["coefficients"]
fig,ax=plt.subplots(figsize=(8,4));ax.scatter(concentration,absorbance,label="Simulated observations")
grid=np.linspace(0,80,101);ax.plot(grid,intercept+slope*grid,label="OLS calibration")
ax.set(xlabel="Declared concentration (µmol/L)",ylabel="Simulated absorbance",title="Beer–Lambert teaching calibration");ax.legend();plt.tight_layout();plt.show()

2. 反算前检查残差

残差为观测减去拟合吸光度。声明的噪声为恒方差高斯噪声；修改模拟后，可用图检查偏离。残差标准差不自动等于可追溯仪器不确定性。

In [5] · Python

fig,ax=plt.subplots(figsize=(8,4));ax.scatter(fitted["fitted"],fitted["residual"])
ax.axhline(0,color="black",linestyle="--");ax.set(xlabel="Fitted absorbance",ylabel="Observed − fitted absorbance",title="Simulated calibration residual audit");plt.tight_layout();plt.show()
print("Fitted absorptivity (L/mol/cm):",slope*1e6/p["path_length_cm"],"; fitted blank:",intercept)

Fitted absorptivity (L/mol/cm): 1526.046117109752 ; fitted blank: 0.020375277723509527

3. 传播校准及未知测量的变异

模拟三次未知样本吸光度。按估计协方差抽取校准系数，并按条件高斯模型抽取未知均值。所得区间依赖校准及噪声假设，不构成普遍实验保证。

In [6] · Python

unknown_truth=p["unknown_concentration_mol_per_L"]*1e6
unknown_replicates=p["blank_absorbance"]+p["absorptivity_L_per_mol_cm"]*p["path_length_cm"]*unknown_truth*1e-6+rng.normal(0,p["noise_sd_absorbance"],3)
unknown_mean=float(unknown_replicates.mean());point=(unknown_mean-intercept)/slope
coefficient_draws=rng.multivariate_normal(fitted["coefficients"],fitted["covariance"],size=5000)
measurement_draws=rng.normal(unknown_mean,np.sqrt(fitted["variance"]/3),5000)
valid=coefficient_draws[:,1]>0
concentration_draws=(measurement_draws[valid]-coefficient_draws[valid,0])/coefficient_draws[valid,1]
interval=np.quantile(concentration_draws,[.025,.975])
fig,ax=plt.subplots(figsize=(8,4));ax.hist(concentration_draws,bins=40,edgecolor="white")
ax.axvspan(*interval,color="orange",alpha=.2,label="Conditional 95% range");ax.axvline(unknown_truth,color="black",linestyle="--",label="Known simulation truth")
ax.set(xlabel="Inverted concentration (µmol/L)",ylabel="Propagation draws",title="Conditional calibration uncertainty");ax.legend();plt.tight_layout();plt.show()
print("Unknown estimate (µmol/L):",point,"; conditional interval:",interval,"; rejected nonpositive slopes:",int((~valid).sum()))

Unknown estimate (µmol/L): 29.025732555623645 ; conditional interval: [24.44653939 33.6977996 ] ; rejected nonpositive slopes: 0

4. 核对拟合并导出全部假设

检查残差正交性，并用直接中心化公式核对斜率。不静默裁剪负浓度抽样，也不强制区间包含已知模拟真值。保存种子、噪声、单位及完整系数协方差。

In [7] · Python

direct_slope=np.sum((concentration-concentration.mean())*(absorbance-absorbance.mean()))/np.sum((concentration-concentration.mean())**2)
assert np.isclose(slope,direct_slope)
assert np.linalg.norm(np.column_stack([np.ones(len(concentration)),concentration]).T@fitted["residual"])<1e-9
assert np.isfinite(concentration_draws).all() and interval[0]<=interval[1]
OUTPUT_DIR.mkdir(parents=True,exist_ok=True)
pd.DataFrame({"concentration_micromol_per_L":concentration,"simulated_absorbance":absorbance,"fitted_absorbance":fitted["fitted"],"residual":fitted["residual"]}).to_csv(OUTPUT_DIR/"simulated-calibration.csv",index=False)
pd.DataFrame({"concentration_draw_micromol_per_L":concentration_draws}).to_csv(OUTPUT_DIR/"conditional-propagation.csv",index=False)
(OUTPUT_DIR/"calibration-model.json").write_text(json.dumps({"seed":models["seed"],"declared_simulation":p,"coefficient_units":["absorbance","absorbance per micromol/L"],"coefficients":fitted["coefficients"].tolist(),"covariance":fitted["covariance"].tolist(),"unknown_point_micromol_per_L":point,"conditional_95_micromol_per_L":interval.tolist()},indent=2))
print("Calibration checks passed; exports:",OUTPUT_DIR)

Calibration checks passed; exports: /Users/zhiluo/Developer/vibeit-cookbook/cookbook/downloads/results/ch06-calibration

自主练习与参考答案

为何校准需要保留截距？
为何传播区间具有条件性？

参考解释。 空白/背景贡献截距。区间依赖线性/噪声模型及估计系数协方差，不保证每次真实实验。

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

可选的 AI 编程练习

提示词 1

仅用 NumPy、Pandas、Matplotlib，在相同浓度设计下比较两个声明噪声水平，记录独立种子。

提示词 2

使用同样的包导出未知重复，分别做仅测量与校准加测量的传播。

人工验收。 明确模拟数据、浓度/吸光度单位、种子及完整协方差，不强制覆盖真值。

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

可选联网扩展

In [8] · Python

RUN_ONLINE_UPDATE=False
if RUN_ONLINE_UPDATE:
    import requests
    try:
        response=requests.get('https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.brentq.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 网页。