企业官网建设流程全解析

连续粉体进料与混合建模教程（10）：RTD 模型①罐串联与分数 n-CSTR——从整数 N 到分数 n：把混合质量变成一个数

版本声明块

• 工具/软件：Python 3.9 + numpy（<2）+ scipy（<1.14，用scipy.special.gammaln）+ matplotlib；E_tis、fit_n_from_variance为本系列二次开发示例函数，非官方 API。
• 语言/环境：Python 3.9；时间一律以秒为基准；θ = t/τ 为无量纲时间，τ 为平均停留时间（s）。
• 本文目标：把"这台连续混合器混得好不好"从一句形容词，变成一个可标定、可比较、可反推操作条件的实数 n。

一句话结论：N 个等体积全混流（continuous stirred tank reactor，CSTR）串联的无量纲 RTD 是 E(θ)=N^{N/(N−1)!·θ}(N−1)·exp(−Nθ)（θ=t/τ），其方差 σ²_θ=1/N 直接给出 N=1/σ²_θ、偏度 γ₁=2/N，把 N 放宽为正实数 n 后，n=1 为理想全混、n>1 为过混合、n<1 为欠混合、n=0.5 出现短路（short-circuiting），某连续混合技术（continuous mixing technology，CMT）器件实测 n 介于 0.5 与 1.1 之间、示例拟合 n=0.95，而数值实现必须用scipy.special.gammaln计算 ln Γ(n) 以防 n 小时的阶乘溢出。

〇、本篇要解决的认知问题

1. N 个等体积 CSTR 串联的 RTD 公式 E(θ)=N^{N/(N−1)!·θ}(N−1)·exp(−Nθ) 是怎么来的？为什么它是伽马分布？
2. 为什么"实验测到的方差"就能反推出罐数 N？这个换算的适用边界在哪？
3. 教科书说"N 必须是整数"，可工业上为什么这个默认值不够用？分数 n 到底在描述什么？
4. n < 1 时 RTD 早期出现的那个尖峰代表什么？为什么叫短路（short-circuiting）？
5. 用scipy.special.gammaln算 Γ(n) 和直接算阶乘，差别真的有那么大吗？n 取 0.5 时会炸吗？

一、机制解析

1.1 从"一个釜"到"N 个釜"：公式的来源

先从一个理想全混流（CSTR）出发。底账登记它的 RTD 是E_tank(t) = exp(−t/τ)/τ——指数衰减，均值 τ，方差 τ²。这个"指数尾巴"意味着：总有一部分物料在釜里待得特别久，这正是返混（back-mixing）的特征。

现在把 N 个等体积釜串起来。因为等体积，单釜停留时间为 τ₁ = τ/N（τ 是整串的平均停留时间）。按第 9 篇的卷积规则，整串 RTD 是 N 个相同指数分布的卷积。数学上，指数分布是伽马分布的形状参数为 1 的特例，而独立同分布伽马变量的和仍是伽马分布、形状参数相加——所以 N 个卷积的结果是形状参数 N、尺度参数 τ₁ 的伽马分布。写成无量纲形式（θ = t/τ），就得到底账登记的标准式：

E(θ) = N^N / (N − 1)! · θ^(N−1) · exp(−N θ) θ = t/τ

底账同时点明：这就是形状参数 N、尺度参数 1/N 的伽马分布。

N = 1（单釜）： N = 4： N → 大： E|\\ E| E| ┌─┐ | \\___ | /\ | / \ | \______ | / \___ | / \___ | \____ |/ \_____ | / \____ +------------------→ θ +----------------→ θ +----------------→ θ 单调下降、长尾（强返混） 单峰、右偏 窄峰（近平推流）

1.2 累积分布与矩：三个能直接用的数

底账给出的累积分布是

F(θ) = 1 − exp(−N θ) · Σ_{k=0}^{N−1} (N θ)^k / k!

括号里的求和正是泊松分布的累积项——这也解释了为什么 TIS 模型与"等待 N 次事件"的过程同构。F(θ) 是第 18 篇做物料追溯与批定义时直接要用的函数：F(θ)=0.95 的那个 θ 就是"95% 的物料已出料"的时间点。

真正让 TIS 变成工程工具的是它的三个矩（底账登记）：

N = 1/σ²_θ 是本章最有价值的一行。它意味着你不需要做复杂的曲线拟合，只要测出 E(t) 的一阶与二阶矩，就立刻得到一个"等效罐数"。这在现场调试时极其高效：一次脉冲试验、一个方差、一个除法。

但要立刻附上边界条件：这条换算假设了 RTD 确实是等体积 TIS 形态。若实测曲线存在明显的双峰或早期尖峰，用 σ²_θ 折出来的 N 会是一个"等效值"，它可能与任何整数都差很远——这恰恰引出 1.4 节。

1.3 工程判读：N 到底多大才算"接近平推流"

底账给出的判读口径非常实用：

注意"N > 10 足够接近平推流"这句话里的限定词：对多数工程目的。这不是数学上的极限，而是工程精度上的取舍。如果你要做的是毫秒级的时间反算（第 18 篇），N = 10 与 N = ∞ 的差别仍然足以影响分流判定——判读要跟着目的走，不能拿一句经验法则通吃。

还要分清两个"方差"：

• 绝对方差σ² 随串联段数增加而增加（第 9 篇已证）；
• 无量纲方差σ²_θ = σ²/τ² = 1/N 随 N 增加而减小。

"更接近平推流"说的是后者。把这两个混为一谈，是最常见的一处概念事故。

1.4 为什么整数 N 不够用：从"罐数"到"混合质量"

标准 TIS 只允许整数N。这个限制来自模型的物理叙事：你有 1 个釜、2 个釜、3 个釜……没有"2.7 个釜"这种设备。

但连续混合器不是釜列。它是一根带螺旋转子的筒，物料在里面的流动既不是一次次"倒进倒出"，也不是一根干净的活塞流——它介于两者之间，而且随操作条件连续变化（底账 §6.1 明确记录：即使平均停留时间 MRT 被出口阀主动控制，RTD 的形状仍随叶轮转速、持料量而变化）。设备形态没变、罐数没变，但 RTD 变了。用一个不变的整数 N 去描述它，等于把连续变量硬塞进整数格子。

于是把 N 放宽到任意正实数 n：

n = 1 → 理想 CSTR，指数型 RTD（完全返混） n > 1 → 过混合（over-mixing），比同等 MRT 的 TIS 更"细"的粒度 n < 1 → 欠混合（under-mixing） n = 0.5 → 混合过程存在短路（short-circuiting），RTD 早期出现尖峰

n = 0.5 时 RTD 为什么会出现早期尖峰？代入公式看形状：θ^(n−1) = θ^(−0.5)，在 θ → 0 时趋于无穷；乘以 exp(−0.5θ) 后，曲线在零点附近被抬成一个陡峭的峰。物理上这代表相当一部分物料走了"快通道"，几乎不停留就出去了——这正是短路的定义。第 13 篇讲偏析、第 16 篇讲扰动传播时，"短路"意味着扰动还没来得及被混合器平滑，就直接穿透到出口，这是最危险的一类流动形态。

而 n > 1 为什么叫"过混合"？因为同样的 MRT 下，n 越大 RTD 越窄，物料停留时间的离散度越小——相当于在同一台设备上做出了更多级串联的效果。注意"过"不等于"好"：过混合往往意味着更高的剪切与能耗，对易碎颗粒（底账记录振动式 LIW 进料器专用于易碎物料，就是为避开机械剪切）反而是伤害。

1.5 一手锚点：CMT 器件的实测 n

底账登记了一组一手实测锚点（★★★ AIChE Annual Meeting 2018, paper 557b）：

请特别注意这组数的含义：一台真实连续混合器的 n，落在 1 附近甚至低于 1。这意味着它远没有达到"多级串联"的混合能力——n = 0.95 甚至略低于理想全混釜。这与很多人的直觉（“工业混合器肯定比一个理想搅拌釜强”）正好相反，也是"MRT 恒定 ≠ RTD 恒定"这条铁律的实证。

更有价值的是它的反推用法：因为不同操作点对应不同 RTD，就可以把 n 反向关联到操作条件（转速、持料量、吞吐），进而用 RTD 建模来划定工艺操作空间，保证混合组成达标。这条线在第 19 篇会展开成完整的设计空间扫描。

操作条件 (CPP) RTD 形状 混合质量 (CQA) ┌────────────────────┐ ┌──────────────────┐ ┌────────────────────┐ │ 螺杆/叶轮转速 ω │ │ │ │ │ │ 持料量 hold-up │ ─标定─▶│ n（或完整 E(t)） │─判据─▶│ 出口 RSD / Lacey │ │ 质量吞吐 ṁ │ │ │ │ 组成是否达标 │ └────────────────────┘ └──────────────────┘ └────────────────────┘ ▲ │ └─────── 反推操作空间 ◀───────┘

1.6 分数 n 的数值陷阱：阶乘必须走 gammaln

公式里有一个(N − 1)!。整数时好办，n = 0.5 时怎么办？答案是阶乘的连续推广就是伽马函数：(n − 1)! = Γ(n)。于是：

E(θ) = n^n / Γ(n) · θ^(n−1) · exp(−n θ)

但直接算n**n / gamma(n)有两个坑：

• 溢出：n 较大时 n^n 与 Γ(n) 都是天文数字，浮点直接 inf/inf = nan；
• n 很小时 Γ(n) → ∞：Γ(0.5) ≈ 1.7725 还好，但 n = 0.01 时 Γ(n) ≈ 99.4，n = 0.001 时约 999，继续减小会溢出；而 n^n → 1、θ^(n−1) 的负幂次又会爆炸。

标准解法是全程在对数域运算：用scipy.special.gammaln(n)得到 ln Γ(n)，先算ln E = n·ln n − lnΓ(n) + (n−1)·ln θ − n·θ，最后再exp回来。这样既避免了中间量溢出，也把"负幂次 × 大分母"的相消合并到一次指数运算里。

顺带一个极易被忽略的除法陷阱：无量纲 E(θ) 的积分是 1（对 θ 积分），但有量纲的 E(t) 还要除以 τ：

E(t) = E(θ) / τ

漏掉这个 1/τ，整条 RTD 的面积会变成 τ 而不是 1，第 9 篇的归一化断言就会报警。

二、完整代码与逐行剖析

第一段代码实现E_tis(t, tau, n)、F_tis与fit_n_from_variance，并对多个 n 做面积归一化断言。

# -*- coding: utf-8 -*-""" 第 10 篇（代码一）：罐串联（TIS）与广义分数 n-CSTR 的 RTD ======================================================= 本脚本属于"连续粉体进料与混合建模教程"系列的二次开发示例，非任何官方 API。 时间基准：秒（s）。theta = t / tau 为无量纲时间。 """importnumpyasnpfromscipy.specialimportgammalndefE_tis(t,tau,n):"""等体积罐串联 RTD（支持整数与分数 n）。 E(theta) = n^n / Gamma(n) * theta^(n-1) * exp(-n*theta)，theta = t / tau 有量纲形式：E(t) = E(theta) / tau 参数 ---- t : 时间（s），标量或数组 tau : 平均停留时间（s），> 0 n : 等效罐数，正实数（n=1 理想 CSTR；n>1 过混合；n<1 欠混合） 说明 ---- 全程在对数域运算，用 gammaln 求 ln Gamma(n)，避免 n**n 与 Gamma(n) 各自溢出。 """ifn<=0:raiseValueError("n 必须为正实数")iftau<=0:raiseValueError("tau 必须为正")theta=np.asarray(t,float)/tau ln_pre=n*np.log(n)-gammaln(n)# ln[ n^n / Gamma(n) ]E_theta=np.zeros_like(theta,dtype=float)ok=theta>0# theta=0 时下面会出现 -inf/nanwithnp.errstate(divide="ignore",invalid="ignore",over="ignore"):ln_E=ln_pre+(n-1.0)*np.log(theta[ok])-n*theta[ok]E_theta[ok]=np.exp(ln_E)returnE_theta/tau# 有量纲：除以 taudefF_tis(t,tau,n):"""等体积罐串联的累积分布 F(theta) = 1 - exp(-n*theta) * sum_{k=0}^{n-1}(n*theta)^k/k!。 注意：k 的上限对分数 n 没有定义，故分数 n 时改用伽马函数的下不完全形式； 此处实现对整数 n 用泊松求和，对分数 n 用数值积分 E_tis（稳健且不引入新公式）。 """t=np.asarray(t,float)iffloat(n).is_integer()andn>=1:N=int(n)theta=t/tau s=np.zeros_like(theta,dtype=float)term=np.ones_like(theta,dtype=float)# k = 0 项forkinrange(N):ifk>0:term=term*(N*theta)/k s=s+termreturn1.0-np.exp(-N*theta)*s# 分数 n：对 E_tis 做累计数值积分（梯形，要求 t 为等步长网格）dt=t[1]-t[0]E=E_tis(t,tau,n)returnnp.cumsum(0.5*(np.concatenate([[0.0],E[:-1]])+E)*dt)deffit_n_from_variance(sigma2_theta):"""由无量纲方差反推等效罐数：N = 1 / sigma2_theta。"""ifsigma2_theta<=0:raiseValueError("sigma2_theta 必须为正")return1.0/sigma2_thetadef_area(t,E):"""梯形积分求面积，用于归一化断言。"""returnfloat(np.sum(0.5*(E[1:]+E[:-1])*np.diff(t)))defmain():tau=30.0# 平均停留时间 30 s# n < 1 时 theta→0 处有陡峰，用小步长 + 大范围的非均匀网格兼顾尖峰与长尾t=np.unique(np.concatenate([np.arange(0.0,3.0,0.005),np.arange(3.0,400.0,0.1)]))print("=== 不同 n 的 TIS 形态与矩自检（tau = %.0f s）==="%tau)fornin(0.5,0.8,1.0,2.0,5.0,10.0):E=E_tis(t,tau,n)a=_area(t,E)MRT=float(np.sum(0.5*(t[1:]*E[1:]+t[:-1]*E[:-1])*np.diff(t)))/a var=float(np.sum(0.5*((t[1:]-MRT)**2*E[1:]+(t[:-1]-MRT)**2*E[:-1])*np.diff(t)))/aprint(f"n={n:4.1f}面积={a:.6f}MRT={MRT:7.3f}s "f"sigma2_theta={var/tau**2:.4f}反推 n={fit_n_from_variance(var/tau**2):6.3f}")assertabs(a-1.0)<1.5e-2,f"n={n}面积{a}归一化失败"assertabs(MRT-tau)<0.2*tau,f"n={n}MRT 偏离 tau 过多"print("\n[自检] 所有 n 的面积归一化与均值 theta=1 均通过")if__name__=="__main__":main()

逐段剖析：

• ln_pre = n * np.log(n) - gammaln(n)是整段代码的心脏。当 n = 100 时，n^n ≈ 1e200、Γ(100) ≈ 9.3e155，两者都在浮点边缘；取对数后分别是 460 与 359，安全地相减得到 101。这就是"对数域运算"的价值：把溢出风险消灭在减法里。
• theta[ok]的掩码必不可少。当 n < 1 时(n−1)为负，θ = 0 处log(0) = −inf，−inf × 负数 = +inf，再exp(+inf) = inf。用掩码把 θ = 0 排除、令 E(0) = 0，是对"点值无定义"的诚实处理；若你的应用需要 E(0) 的极限值，应显式单独计算而不是让它污染数组。
• return E_theta / tau是把无量纲分布折算成有量纲 RTD 的唯一一步，也是归一化断言能否通过的关键。
• F_tis对分数 n 走数值积分而非硬套泊松求和。原因很直接：Σ_{k=0}^{N−1}的上限对 n = 0.5 无意义。底账只登记了整数形式的 F(θ)，因此分数情形的闭式以下不完全伽马函数的原文口径为准；本篇用数值积分规避，不编造新公式。
• fit_n_from_variance只有一行除法，但它把"实验测方差 → 等效罐数"这条工程捷径固化成了可复用接口，第 17 篇标定时会反复调用。
• 断言用 1.5e-2 的面积阈值：因为 n = 0.5 的曲线在 θ = 0 附近有陡峭尖峰，即便前 3 s 用了 5 ms 的细步长，面积仍会偏低约 0.8%。这不是 bug，是离散化的必然代价——想更准就继续加密网格。实测参考值：n=0.5 面积约 0.992、n=0.8 约 0.9995、n≥1 优于 0.9999。

第二段代码画 n 扫描族图，并把 n = 0.5 的短路尖峰单独标出来。

# -*- coding: utf-8 -*-""" 第 10 篇（代码二）：n 扫描 RTD 族与"短路尖峰"可视化 =============================================== """importnumpyasnpfromscipy.specialimportgammalnimportmatplotlib matplotlib.use("Agg")importmatplotlib.pyplotaspltdefE_tis(t,tau,n):"""同代码一：分数 n 罐串联 RTD（对数域运算防溢出）。"""ifn<=0ortau<=0:raiseValueError("n 与 tau 必须为正")theta=np.asarray(t,float)/tau ln_pre=n*np.log(n)-gammaln(n)out=np.zeros_like(theta,dtype=float)ok=theta>0withnp.errstate(divide="ignore",invalid="ignore",over="ignore"):out[ok]=np.exp(ln_pre+(n-1.0)*np.log(theta[ok])-n*theta[ok])returnout/taudefmain():tau=30.0# n = 0.5 在 theta 很小时有陡峰，故用非线性加密网格：前 5 s 步长 0.002 st=np.concatenate([np.arange(0.0,5.0,0.002),np.arange(5.0,400.0,0.2)])t=np.unique(t)# 去重并升序n_list=[0.5,0.8,1.0,2.0,5.0,10.0]fig,ax=plt.subplots(figsize=(8.4,4.4))forninn_list:E=E_tis(t,tau,n)# 面积自检：n=0.5 的尖峰需要足够密的网格a=float(np.sum(0.5*(E[1:]+E[:-1])*np.diff(t)))assertabs(a-1.0)<0.02,f"n={n}面积{a}偏离过大"ax.plot(t,E*tau,label=f"n ={n}")# 画无量纲 E(theta)print(f"n={n:4.1f}面积={a:.5f}E(theta=0.05)={np.interp(0.05,t/tau,E*tau):.4f}")ax.set_xlim(0,3.0)# 只看 theta ∈ [0, 3]ax.set_xlabel("无量纲时间 theta = t / tau")ax.set_ylabel("E(theta)")ax.set_title("分数 n-CSTR 族：n<1 欠混合(短路)、n=1 全混、n>1 过混合")ax.legend(fontsize=9)ax.grid(alpha=0.3)fig.tight_layout()fig.savefig("fig10_n_scan.png",dpi=130)print("已输出 fig10_n_scan.png")print("[提示] n=0.5 在 theta→0 处抬升，即短路(short-circuiting)的早期尖峰")if__name__=="__main__":main()

这段代码的两个要点：

1. 网格加密是 n = 0.5 的必要条件。θ^(−0.5) 在 θ → 0 处发散，若用均匀网格 dt = 0.2 s，前几个点会严重漏采，面积甚至可能偏离 3% 以上。代码用np.arange(0.0, 5.0, 0.002)在前 5 s 加密，把尖峰采出来——这是"分数 n 的数值陷阱"在实际计算中的具体形态。
2. 画的是 E(θ) 而不是 E(t)。n 扫描的意义在于比较形状，而无量纲化才能把不同 τ 的设备放到同一张图上比较。若画E(t)，曲线的横向位置会随 τ 漂移，掩盖形状差异。

三、常见报错与排查

现象 1：E_tis在 n = 0.5 时返回大量nan或inf。
根因：θ = 0 处的log(0) × (n−1)产生inf，再exp溢出。
解法：用theta > 0掩码把零点排除；配合np.errstate(..., over="ignore")抑制告警；若确需 E(0) 的极限，单独解析计算。

现象 2：面积断言失败，n 越大偏得越多。
根因：漏掉E_theta / tau这一步（面积变成 τ），或时间上限 t_max 太小导致右尾截断。
解法：确认返回除 τ；把 t 上限扩到 ≥ 8τ 复测面积是否收敛到 1。

现象 3：RuntimeWarning: overflow encountered in power（用n**n / gamma(n)直接算时）。
根因：中间量溢出，n 稍大就 inf/inf。
解法：改用gammaln在对数域运算；这是本篇必须掌握的写法，scipy.special.gammaln是标准工具。

现象 4：用实验方差反推的 n 是个奇怪的数（比如 0.37 或 3.8）。
根因：N = 1/σ²_θ假设 RTD 是等体积 TIS 形态；若实测曲线有双峰、早期尖峰或长拖尾，这个换算只是"等效值"。
解法：先画 E(t) 目视检查形状；若有明显短路尖峰，说明应走分数 n 而非整数 N；再用第 17 篇的curve_fit做完整标定并报告 1σ 误差。

现象 5：把不同文献的 n 值直接横比。
根因：n 与 τ（或 MRT）是耦合的，且不同文献的归一化与截断口径不同。
解法：比较前先统一：时间单位（s）、τ 的确定方式（持料量/吞吐 还是拟合值）、尾部截断位置。只比"是否显著小于 1（短路倾向）"这类方向性结论，不比数值。

四、动手练习

练习 1（基础）：面积与均值自检。
任务：对 n = 0.5、1、2、5、10 各算一遍E_tis(t, tau=30, n)，输出面积与 MRT。
判定标准：所有 n 的面积相对误差 < 0.5%（n = 0.5 放宽到 2%），MRT 与 τ = 30 s 的相对误差 < 2%。

练习 2（进阶）：由方差反推 n，闭环验证。
任务：对每个 n 用数值算出 σ²_θ，再喂给fit_n_from_variance反推 n̂。
判定标准：n̂ 与真值 n 的相对误差 < 5%；n ≥ 1 时应优于 2%。把 n 与 n̂ 的对照表打印出来，并指出误差随 n 增大而上升的原因（网格分辨率与尾部截断）。

练习 3（综合）：短路尖峰的网格敏感性。
任务：对 n = 0.5，分别用均匀网格 dt = 0.5、0.2、0.05、0.01 s 计算面积。
判定标准：列出"步长 → 面积"表格，明确指出要让面积误差 < 1%，步长需要小于多少；并解释为什么 n = 0.5 比 n = 10 对步长敏感得多（θ^(n−1) 的幂指数）。

五、小结与下一篇预告

罐串联模型把"混合"这件事压缩成两个数：τ 刻画"平均待多久"，n 刻画"待多久的离散程度"。整数 N 来自"釜列"的物理叙事，而工业连续混合器需要的是一根连续变化的旋钮，于是有了分数 n：n = 1 理想全混、n > 1 过混合、n < 1 欠混合、n = 0.5 出现短路尖峰。一手锚点提醒我们别高估设备——某 CMT 器件的实测 n 只有 0.5 到 1.1（示例 0.95）。数值上必须用scipy.special.gammaln在对数域算n^n/Γ(n)，并在 n < 1 时对 θ 的小值区加密网格。

TIS 是一个"离散化"的世界观：把设备切成一个个理想釜。下一篇换一个完全不同的视角——第 11 篇的轴向扩散（axial dispersion）模型：不切釜，而是把非理想性归给一个类 Fick 扩散项，用 Peclet 数 Pe = uL/D_a 衡量对流与扩散的比值，讨论开-开边界（Taylor 1953）与闭-闭边界（Danckwerts 条件）怎么选，并给出N ≈ Pe/2 + 1这条只在 N > 5 时成立的换算。TIS 的 n 与轴向扩散的 Pe 会在这条换算上相遇——但请记住，N 小时两者对不上，这正是第 12 篇要用 GEA CDC-50 的组合模型去跨越的鸿沟。

本篇认知问题回显（FAQ）

Q1：N 个等体积 CSTR 串联的 RTD 公式是怎么来的，为什么是伽马分布？
A：它等于 N 个相同指数分布（单釜 E(t)=exp(−t/τ₁)/τ₁，τ₁=τ/N）的卷积；独立伽马变量相加时形状参数相加，故结果为形状参数 N、尺度参数 1/N 的伽马分布，无量纲形式即 E(θ)=N^{N/(N−1)!·θ}(N−1)·exp(−Nθ)，θ=t/τ。

Q2：怎么由实验方差反推罐数 N，边界在哪？
A：等体积 TIS 的无量纲方差 σ²_θ=1/N，故 N=1/σ²_θ，测出一二阶矩即可估计；该换算只在实测 RTD 接近等体积罐串联形态时成立，若有双峰、早期尖峰或长拖尾，得到的只是"等效罐数"而非真实罐数。

Q3：教科书要求 N 是整数，为什么工业上要用分数 n？
A：连续混合器的 RTD 形状随转速、持料量等操作条件连续变化，用固定整数 N 无法刻画；放宽到正实数 n 后，n 本身就成为描述混合质量的可标定参数，n 可反推操作条件并用于划定工艺操作空间。

Q4：n-CSTR 模型里 n 小于 1 时 RTD 为什么出现早期尖峰？
A：n<1 时 E(θ) 含 θ^(n−1) 的负幂次，在 θ→0 处发散并被 exp(−nθ) 压成陡峰，代表一部分物料几乎不停留就出料，即短路（short-circuiting）；此时 n=0.5 是典型值，扰动能直接穿透混合器。

Q5：用 scipy.special.gammaln 算 Γ(n) 和直接算阶乘差多少，n=0.5 会溢出吗？
A：差别在于中间量：直接算 n^n/Γ(n) 在 n 稍大时 n^n 与 Γ(n) 都已接近浮点上限而溢出成 nan，而先取对数相减再 exp 全程安全；n=0.5 时 Γ(0.5)≈1.7725 不会溢出，但 θ=0 处 log(0) 会产生 inf，必须用掩码单独处理。