## Inference for Stan model: 6dcfbccbf2f063595ccc9b93f383e221. ## 4 chains, each with iter=5000; warmup=2500; thin=1; ## post-warmup draws per chain=2500, total post-warmup draws=10000. ## ## mean se_mean sd 10% 90% n_eff Rhat ## theta 0.27 0.00 0.09 0.16 0.39 3821 1 ## lp__ -13.40 0.01 0.73 -14.25 -12.90 3998 1 ##
# 提取后验抽样 # 计算后均值（估计） mean(theta_draws)
## [1] 0.2715866
# 计算后验区间 ## 10% 90%## 0.1569165 0.3934832 ggplot(theta_draws_df, aes(x=theta)) +geom_histogram(bins=20, color="gray")
RStan：MAP，MLE
Stan的估算优化；两种观点：
- 最大后验估计(MAP)。
- 最大似然估计（MLE）。
optimizing(model, data=c("N", "y")) ## $par ## theta ## 0.4 ## ## $value ## [1] -3.4 ## ## $return_code ## [1] 0
- 种群竞争模型 ---Lotka-Volterra模型
- 洛特卡（Lotka，1925）和沃尔泰拉（Volterra，1926）制定了参数化微分方程，描述了食肉动物和猎物的竞争种群。
- 完整的贝叶斯推断可用于估计未来（或过去）的种群数量。
- Stan用于对统计模型进行编码并执行完整的贝叶斯推理，以解决从噪声数据中推断参数的逆问题。
在此示例中，我们希望根据公司每年收集的毛皮数量，将模型拟合到1900年至1920年之间各自种群的加拿大猫科食肉动物和野兔猎物。
数学模型
我们表示U（t）和V（t）作为猎物和捕食者种群数量 分别。与它们相关的微分方程为：
这里：
- α：猎物增长速度。
- β：捕食引起的猎物减少速度。
- γ：自然的捕食者减少速度。
- δ：捕食者从捕食中增长速度。
stan中的Lotka-Volterra
real[] dz_dt(data real t, // 时间 real[] z, // 系统状态 real[] theta, // 参数 data real[] x_r, // 数值数据 data int[] x_i) // 整数数据 { real u = z[1]; // 提取状态 real v = z[2];
观察到已知变量：
-
：表示在时间 
的 
物种数量
必须推断未知变量）：
- 初始状态：
：k的初始物种数量。 - 后续状态
：在时间t的物种数量k。 - 参量
。
假设误差是成比例的（而不是相加的）：
等效：
与
建立模型
已知常数和观测数据的变量。
data { int<lower = 0> N; // 数量测量 real ts[N]; // 测量次数>0 real y0[2]; // 初始数量 real<lower=0> y[N,2]; // 后续数量 }
未知参数的变量。
parameters { real<lower=0> theta[4]; // alpha, beta, gamma, delta real<lower=0> z0[2]; // 原始种群 real<lower=0> sigma[2]; // 预测误差 }
先验分布和概率。
model { // 先验 sigma ~ lognormal(0, 0.5); theta[{1, 3}] ~ normal(1, 0.5); // 似然（对数正态） for (k in 1:2) { y0[k] ~ lognormal(log(z0[k]), sigma[k]);
我们必须为预测的总体定义变量 ：
- 初始种群（
z0）。 - 初始时间（
0.0），时间（ts）。 - 参数（
theta）。 - 最大迭代次数（
1e3）。
Lotka-Volterra参数估计
print(fit, c("theta", "sigma"), probs=c(0.1, 0.5, 0.9))
获得结果：
mean se_mean sd 10% 50% 90% n_eff Rhat ## theta[1] 0.55 0 0.07 0.46 0.54 0.64 1168 1 ## theta[2] 0.03 0 0.00 0.02 0.03 0.03 1305 1 ## theta[3] 0.80 0 0.10 0.68 0.80 0.94 1117 1 ## theta[4] 0.02 0 0.00 0.02 0.02 0.03 1230 1 ## sigma[1] 0.29 0 0.05 0.23 0.28 0.36 2673 1 ## sigma[2] 0.29 0 0.06 0.23 0.29 0.37 2821 1
分析所得结果：
- Rhat接近1表示收敛；n_eff是有效样本大小。
- 10％，后验分位数；例如
。 - 后验均值是贝叶斯点估计：α=0.55。
- 后验平均估计的标准误为0。
- α的后验标准偏差为0.07。
