数学公式渲染测试

数学公式显示测试

如果炸弹每秒钟爆炸概率提高一点,从数学期望上来看最有可能在哪一秒爆炸?

结论是：在第 10 秒爆炸概率最大,爆炸时间的期望是 12.5331 秒。
其实可以给出更为普适的答案：假设有很多个互不影响的炸弹,每个未炸炸弹在下一秒被引爆的概率期望是

πm2⋅erf⁡(m2)−mem/2 \sqrt{\frac{\pi m}{2}} \cdot \operatorname{erf}\left( \sqrt{\frac{m}{2}} \right) - \frac{m}{e^{m/2}}

这题在很多物理化学问题中都有实际应用：假设有某个化学反应,其反应速率正比于反应物剩余浓度C(相当于炸弹未爆炸的概率)和温度T(实际上随T指数增长,不过原理差不多)。那么在升温时,随时间的增长,一方面反应速率会先因温度的升高而升高；但另一方面,随着反应的进行,反应物浓度越来越低,逐渐限制反应进行,因此反应速率在某个温度达到峰值后,会开始随时间降低。典型的热脱附谱 (Thermal desorption spectroscopy) 峰值的温度,来推算出反应的能量学信息。
典型的热脱附结果如下：

先来看离散情形吧，也就是把时间拆成一秒一秒来看。这种情形最简单，高中数学就能解决。
根据题设可知,某个炸弹在1~n-1秒都没爆炸的概率为

∏i=1n−1(1−im) \prod_{i=1}^{n-1} \left(1 - \frac{i}{m}\right)

令炸弹在第 ( n ) 秒爆炸的概率为 ( P(n) ) ,有：

P(n)=nm∏i=1n−1(1−im) P(n) = \frac{n}{m} \prod_{i=1}^{n-1} \left(1 - \frac{i}{m}\right)

显然

P(n+1)P(n)=n+1m⋅m−nm⋅∏i=1n−1(1−im)nm⋅∏i=1n−1(1−im)=(n+1)(m−n)mn \frac{P(n+1)}{P(n)} = \frac{\frac{n+1}{m} \cdot \frac{m-n}{m} \cdot \prod_{i=1}^{n-1} \left(1 - \frac{i}{m}\right)}{\frac{n}{m} \cdot \prod_{i=1}^{n-1} \left(1 - \frac{i}{m}\right)} = \frac{(n+1)(m-n)}{mn}

该比值函数在研究范围内是个单调递减函数,所以当 P(n+1)P(n)=1 \frac{P(n+1)}{P(n)} = 1 时,P(n)取得最大值。
对应的解为 n=m+14−12 n = \sqrt{m + \frac{1}{4}} - \frac{1}{2} 。当m比较大时,可以近似认为 n=m n = \sqrt{m} 。
简单写了个程序(代码见答案底部),当 ( m = 100 )时,数值解出来的爆炸期望为 12.21,概率随时间变化的函数长这样子,是不是和上面的热脱附谱长得很像？
再来看连续情形。假设炸弹在 0→t 秒内不爆炸的概率为 C(t),那么将炸弹在 t→t+dt时间内爆炸的概率，记为 p(t)，它等于p(t)=C(t)−C(t+dt)=−dC(t) p(t) = C(t) - C(t+dt) = -dC(t)
其实题中问的是单位时间内的爆炸概率,即概率密度,记为 ( P(t) )：

P(t)=p(t)dt=−dC(t)dt P(t) = \frac{p(t)}{dt} = -\frac{dC(t)}{dt}

按照题中定义,它正比于时间 ( t ) ,以及未爆炸概率 ( C(t) ),

P(t)=tmC(t) P(t) = \frac{t}{m} C(t)

将这两个式子结合起来并积分,得到:
如果炸弹每秒钟爆炸概率提高一点,从数学期望上来看最有可能在..
可以看出 ( C(t) ) 是个正态分布, ( P(t) ) 是个瑞利分布 ( ^{+} ) (Rayleigh distribution)
对 ( P(t) ) 求导,使其等于零,得爆炸概率密度的最大值

dP(t)dt=C(t)m(1−t2m)=0 \frac{dP(t)}{dt} = \frac{C(t)}{m} \left(1 - \frac{t^{2}}{m}\right) = 0

解得 t=m t = \sqrt{m} ,与m较大时的离散解是一致的。
积分可得爆炸的时间期望:

E=∫0mtP(t)dt=πm2⋅erf⁡(m2)−mem/2 E = \int_{0}^{m} t P(t) dt = \sqrt{\frac{\pi m}{2}} \cdot \operatorname{erf}\left( \sqrt{\frac{m}{2}} \right) - \frac{m}{e^{m/2}}

当 ( m = 100 ) 时,时间期望 ( E = 12.5331 ) ,与离散情形稍有出入,应该是离散/连续模型中时间分割细致度不一样导致的。
相关收藏：认真写科普

上面数值计算的 MATLAB 代码为:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
clear;clc;close all;

days=100;
P_explode(1:days)=0;
P_explode(1)=1/days;

P_unexplode(1:days)=0;
P_unexplode(1)=1-1/days;

E=P_explode(1);

for n=2:1:days;
    P_unexplode(n)=P_unexplode(n-1)*(1-n/days);
    %第n秒不爆炸的概率
    P_explode(n)=P_unexplode(n-1)*(n/days);
    %第n秒爆炸的概率
    E=E+P_explode(n)*n;
end

E=E

set(0,'defaultfigurecolor','w');
figure
set(gcf,'Position',[300 300 800 600]);
set(gca,'Position',[.10 .10 .85 .85]);
plot(1:days,P_explode,'b-','linewidth',1.5);
lx=xlabel('Time (s)','FontSize',16);
ly=ylabel('Explosion rate','FontSize',16);