贝叶斯估计求取分布参数Python

贝叶斯估计是一种求取分布参数的方法，通过先验分布和观测数据的似然度来计算后验分布，并以参数在后验分布下的期望来估算参数值。在Python中，可以使用scipy库中的stats模块来进行贝叶斯估计。具体步骤如下：
1. 导入所需的库：```
import numpy as np
from scipy import stats

2. 定义先验分布：根据你的先验知识和问题需求，选择合适的先验分布，并设置其参数。
3. 生成观测数据：根据实际问题，生成观测数据集。
4. 计算后验分布：根据观测数据和先验分布，使用贝叶斯公式计算后验分布。这里可以使用stats模块中的相应函数来计算。
5. 估算参数值：根据后验分布，计算参数在后验分布下的期望值作为参数的估计值。

相关问题

贝叶斯推理基于gamma分布的先验

### 贝叶斯推理中使用 Gamma 分布作为先验的方法

在贝叶斯统计框架内，参数被视作随机变量，并通过先验分布表达对其初始信念。当考虑Gamma分布作为先验时，这通常适用于那些取正值的未知参数估计场景。

#### 为什么选择Gamma分布？

Gamma分布在建模非负连续型数据方面非常有用，尤其适合用于描述事件发生率或等待时间等问题中的不确定性。其灵活性来源于两个形状参数\(k\)（有时也称为\(\alpha\)）和尺度参数\(\theta\)（或逆尺度参数\(\beta=1/\theta\)），允许调整分布形态以匹配具体应用场景的需求[^3]。

对于某些特定类型的似然函数来说，比如泊松过程中的强度参数λ,选用gamma(α,β)形式的共轭先验能极大简化后续计算流程；因为此时后验也会保持同样的分布族特性，从而使得解析解成为可能。

#### 实际操作案例

假设有这样一个简单例子：已知某网站平均每小时访问次数服从Poisson分布，则该平均速率μ可由Gamma分布来充当合适的先验模型：

\[ \mu|\alpha,\beta \sim \text{Gamma} (\alpha ,\beta ) \]

其中，
- \(E[\mu]=\frac{\alpha}{\beta}\),
- Var(\(\mu)=\frac{\alpha}{\beta ^2}\).

现在观察到了一组样本观测值y_1,...,y_n (每个小时内的实际点击量)，基于最大似然原则更新后的后验分布同样遵循伽玛分布规律:

\[ p(\mu|y_{1:n},\alpha',\beta')=\prod _i^{n}[p(y_i|\mu)p(\mu)]\propto [\mu^{\sum y_i+\alpha -1}]e^{-\mu(n+\beta)} \]

因此新的超参数变为：
- 更新后的shape:\(\alpha'=\alpha +\sum y_i\)
- 更新后的rate:\(\beta'= n+\beta\)

这意味着随着更多信息的到来，原始设定下的不确定度逐渐减小，最终趋向于更加精确的结果。

import numpy as np
from scipy.stats import gamma

# 假设初始先验为 shape=2, rate=0.5
prior_shape = 2  
prior_rate = 0.5 

# 模拟一些观测数据
observed_data = [np.random.poisson(lam=4., size=1)[0] for _ in range(10)]

# 计算后验分布的新参数
posterior_shape = prior_shape + sum(observed_data)
posterior_rate = len(observed_data) + prior_rate

print(f"Posterior Shape: {posterior_shape}")
print(f"Posterior Rate : {posterior_rate}")

# 绘制前后对比图
import matplotlib.pyplot as plt

x = np.linspace(gamma.ppf(0.01, a=prior_shape, scale=1/prior_rate), 
                gamma.ppf(0.99, a=posterior_shape, scale=1/posterior_rate), 100)

plt.plot(x, gamma.pdf(x, a=prior_shape, scale=1/prior_rate),'r-', lw=5, alpha=0.6, label='Prior')
plt.plot(x, gamma.pdf(x, a=posterior_shape, scale=1/posterior_rate),'b-', lw=5, alpha=0.6, label='Posterior')

plt.legend()
plt.show()