问题提出
假设一次实验结果只有两种，还没有去做实验收集数据，这时有理由相信每个事件发生的概率是一个均匀分布
（这个问题中，事件发生的概率本身是一个有待估计的随机变量）
先验分布是这样的：
\(\left\{ \begin{array}{l}
P(x=0)=\theta\\
P(x=1)=1-\theta\\
\theta\sim U(0,1)
\end{array}\right.\)
(有理由相信，未得到实验数据之前，先验的参数服从均匀分布)
问题解决
根据贝叶斯公式
$P(\theta\mid X)$
$=\dfrac{P(X\mid\theta)P(\theta)}{\int_{-\infty}^{+\infty} P(X\mid\theta )P(\theta)d\theta}$
假设经过两次实验，结果是(1,1)，那么
$P(X\mid\theta)=(1-\theta)^2$
$P(\theta)=1,\forall \theta\in(0,1)$
结果是$P(\theta\mid X)=3(1-\theta)^2$
数值计算
def func_prob_theta(theta):
'''
p(\theta) 先验分布
'''
if 0 < theta < 1:
return 1
else:
return 0
def func_prob_X_given_theta(X, theta):
'''
$P(X\mid\theta)$
:param X: list，每次实验的结果，例如[1,0,0,1,1]
:param theta:
:return:
例如，X=[1,0,0,1,0,0]
func_prob_X_given_theta(X,0.5)，代表均匀0-1分布下，X这个实验结果出现的概率
'''
prob_X_given_theta = func_prob_theta(theta)
for x in X:
prob_X_given_theta *= theta ** (1 - x) * (1 - theta) ** x
return prob_X_given_theta
# 假如某次实验结果是这样的：
X = [1, 1, 0, 0, 1, 1]
# 先计算分母，这是一个积分：
from scipy import integrate
func = lambda theta: func_prob_X_given_theta(X, theta)
S, _ = integrate.quad(func, 0, 1)
# 计算并画出Theta的概率密度分布图
import numpy as np
list_p_theta = []
Theta = np.arange(0, 1, 0.01)
for theta in Theta:
list_p_theta.append(func(theta) / S)
import matplotlib.pyplot as plt
plt.plot(Theta, list_p_theta)
plt.title(X)
plt.show()
上结果：(标题是此次实验结果)
这是20个1和10个0的实验结果
这是200个1和100个0的实验结果
火车头问题
选自弗雷德里克·莫斯泰勒《五十个概率难题的解法》
铁路上火车头以1到N命名，你看到3个火车头，标号分别为30,60,90. 请估计铁路上有多少火车头
1. 假设先验为均匀分布
$\theta \sim U(90,N)$
X:看到标号 30,60,90
所以，\(P(X\mid \theta)=\left\{
\begin{array}{ll}
\dfrac{1}{\theta^3}&\theta\geq 90\\
0&o/w
\end{array}\right.\)
$P(\theta)=\dfrac{1}{N-90}$
计算得到 $P(\theta \mid X)=\dfrac{P(X\mid \theta)P(\theta)}{\int_{-\infty}^{+\infty} P(X\mid \theta)P(\theta)}=\dfrac{1/\theta^3}{1/(2\cdot 90^2)-1/(2N^2)}$
我们想知道平均可能有多少个火车头，也就是 $E(\theta \mid X)$
用期望的计算公式$E(\theta\mid X)=\int \theta P(\theta \mid X)d\theta=\dfrac{1/N-1/90}{1/(2\cdot 90^2)-1/(2N^2)}$
2. 假设为指数分布
我们相信，企业的规模服从幂律分布，也就是说，概率密度函数是$f(\theta)=1/x^\alpha$
进一步，我们近似假设\(f(\theta)=\left\{\begin{array}{l}
\lambda \cdot 1/\theta & 90\leq \theta \leq N\\
0
\end{array}\right.\)，其中 $\lambda$是使得全事件概率为1的因子。这个近似假设在实际中是合理的，因为总有一个规模的上限。
用代码模拟计算$E(\theta\mid X)$，其中贝叶斯公式$P(\theta \mid X)=\dfrac{P(X\mid \theta)P(\theta)}{\int_{-\infty}^{+\infty} P(X\mid \theta)P(\theta)}$
N=500 # 试试增加N，结果收敛
def func_prior_prob(theta):
return 1/theta if 90<=theta<=N else 0
def func_X_given_theta(theta):
return 1/theta**3
prob_X=sum([func_X_given_theta(i)*func_prior_prob(i) for i in range(90,N+1)])
sum([func_X_given_theta(i)*func_prior_prob(i)*i/prob_X for i in range(90,N+1)])
(二战期间，英美统计了德国坦克编码规则和缴获坦克的编码，使用类似方法，计算得到了有巨大价值的情报)
