Bayesian Statistics

HJ. S.

2023-12-11

贝叶斯统计（Bayesian Statistics）定义

贝叶斯统计是一种统计学方法，它基于贝叶斯定理，用于更新随着证据累积而变化的概率估计。与频率主义统计不同，贝叶斯统计将概率视为不确定性的度量，并允许在分析中使用先验知识。

核心：贝叶斯定理，它描述了先验概率（先验信念）如何通过新的证据（数据）更新为后验概率。
公式：
- \(P(\theta | \text{data}) = \frac{P(\text{data} | \theta) \times P(\theta)}{P(\text{data})}\)
- 其中 \(P(\theta | \text{data})\) 是后验概率，\(P(\text{data} | \theta)\) 是似然，\(P(\theta)\) 是先验概率。

故事：霍格沃茨的魔法石猜想

在霍格沃茨，有一个关于魔法石的古老传说。据说，魔法石被藏在学校的某个秘密房间里。教授们对这个传说的真实性持怀疑态度，但一些学生坚信它的存在。

教授的先验信念

教授们基于历史记载和逻辑推理，认为找到魔法石的可能性非常低（低先验概率）。

新证据出现

然而，在一次考古发掘中，学生们发现了一张古老的地图，上面标记着一个神秘的房间。

更新信念：贝叶斯更新

教授们用贝叶斯统计来更新他们关于魔法石存在的信念。他们考虑了新的证据（地图）对原有信念的影响。

似然：地图的出现显著增加了魔法石存在的可能性。
后验概率：结合先验信念和新的证据，教授们更新了他们对魔法石存在可能性的估计，后验概率变得更高。

结论

最终，虽然没有找到确凿的证据证明魔法石的存在，但通过贝叶斯方法，教授们对魔法石可能真实存在的信念有所增强。

实践中的应用

在现实世界的应用中，贝叶斯统计可以用于处理各种不确定性问题，包括疾病诊断、风险评估和机器学习。它允许我们将先验知识与新数据结合起来，形成对问题的全面理解。

贝叶斯统计的一个关键优势是其灵活性和包容性，能够在分析中考虑到不确定性和主观判断。然而，选择合适的先验和正确解释后验概率需要专业知识和实践经验。

贝叶斯统计进阶

继续深入理解贝叶斯统计之后，你可以探索更多与贝叶斯方法相关的高级主题和应用领域。以下是一些建议的后续学习方向：

高级贝叶斯方法

贝叶斯网络（Bayesian Networks）:
- 学习如何使用贝叶斯网络建模复杂的概率关系，尤其在存在不确定性和条件依赖时。
MCMC（Markov Chain Monte Carlo）方法:
- 理解并应用MCMC方法进行复杂后验分布的抽样，这在计算后验概率分布时非常有用。
贝叶斯层次模型（Hierarchical Bayesian Models）:
- 探索如何在数据具有自然层次结构时应用贝叶斯方法。
贝叶斯模型比较:
- 学习如何比较不同贝叶斯模型的优劣，包括使用贝叶斯因子（Bayes Factor）。

应用领域

机器学习和数据科学:
- 贝叶斯方法在机器学习领域有广泛应用，包括分类、回归和聚类等任务。
生物统计学和流行病学:
- 在医学研究中，贝叶斯统计用于药物试验、疾病风险评估和基因数据分析。
金融模型:
- 在金融市场分析和风险管理中应用贝叶斯方法进行预测和决策。

工具和技术

R和Python中的贝叶斯分析包:
- 学习如何使用如Stan、JAGS、PyMC3等工具进行贝叶斯分析。
实际案例研究:
- 通过分析真实数据集来应用和巩固贝叶斯统计的知识。

贝叶斯网络（Bayesian Networks）定义

贝叶斯网络，也称为信念网络或因果网络，是一种图形模型，用于表示变量间的概率关系。它由节点（代表随机变量）和有向边（代表这些变量之间的依赖关系）组成。贝叶斯网络使我们能够进行复杂的概率推理和决策。

贝叶斯网络的关键特点

因果关系:
- 节点间的边代表因果关系或影响。
条件概率:
- 每个节点都有一个条件概率表，它定义了在给定其父节点状态下该节点状态的概率。
推理:
- 贝叶斯网络支持基于新证据更新信念的推理过程。

贝叶斯与霍格沃茨密室

背景

霍格沃茨学校陷入了密室事件的恐慌中。几位学生和图书馆管理员的猫神秘消失了。校方希望使用贝叶斯网络来揭示背后的真相。

构建贝叶斯网络

教授们构建了一个贝叶斯网络来表示不同因素之间的概率关系。网络中的节点包括：

学生活动（如上课、在图书馆等）
目击证言（如听到奇怪的声音、看到可疑人物）
密室相关证据（如密室的开启、怪异的蛇形符号）
消失事件（学生和猫的失踪）

收集证据

图书馆管理员报告说在失踪前夜听到了低语。
部分学生称在图书馆附近看到了蛇形符号。

使用贝叶斯网络进行推理

教授们输入了上述证据到网络中。
网络分析显示，在某些特定条件下，密室的存在与消失事件之间存在较高的概率联系。

R语言实现

这里我们用R语言简化地模拟构建这样一个贝叶斯网络。

# 安装和加载bnlearn包
if (!require("bnlearn")) install.packages("bnlearn")

## Loading required package: bnlearn

library(bnlearn)

我们首先创建了一个包含所有需要的节点的空网络。随后，我们使用set.arc函数添加了从各个因素（学生活动、目击证言和密室相关证据）到消失事件的边。

# 创建一个空的网络
net <- empty.graph(nodes = c("StudentActivity", "Disappearance", "WitnessTestimony", "ChamberEvidence"))

# 逐步添加边来定义网络结构
net <- set.arc(net, from = "StudentActivity", to = "Disappearance")
net <- set.arc(net, from = "WitnessTestimony", to = "Disappearance")
net <- set.arc(net, from = "ChamberEvidence", to = "Disappearance")

# 查看网络结构
print(net)

## 
##   Random/Generated Bayesian network
## 
##   model:
##    [StudentActivity][WitnessTestimony][ChamberEvidence]
##    [Disappearance|StudentActivity:WitnessTestimony:ChamberEvidence]
##   nodes:                                 4 
##   arcs:                                  3 
##     undirected arcs:                     0 
##     directed arcs:                       3 
##   average markov blanket size:           3.00 
##   average neighbourhood size:            1.50 
##   average branching factor:              0.75 
## 
##   generation algorithm:                  Empty

bnlearn包不支持字符型（character）变量。在bnlearn中，节点变量通常需要是因子（factor）类型。在这段代码中，我们使用as.factor函数将所有字符型列转换为因子。这样，bn.fit函数就可以正确处理数据，并根据提供的网络结构和数据来学习网络参数。

贝叶斯网络结果分析

# 模拟一些数据
set.seed(123)  # 为了结果的可重复性
data <- data.frame(
  StudentActivity = sample(c("Class", "Library"), 100, replace = TRUE),
  WitnessTestimony = sample(c("HeardWhispers", "HeardNothing"), 100, replace = TRUE),
  ChamberEvidence = sample(c("SawSymbols", "SawNothing"), 100, replace = TRUE),
  Disappearance = sample(c("Occurred", "NotOccurred"), 100, replace = TRUE)
)

# 学习网络的参数
# 确保所有列都是因子类型
data$StudentActivity <- as.factor(data$StudentActivity)
data$WitnessTestimony <- as.factor(data$WitnessTestimony)
data$ChamberEvidence <- as.factor(data$ChamberEvidence)
data$Disappearance <- as.factor(data$Disappearance)

# 使用转换后的数据拟合贝叶斯网络
fitted_net <- bn.fit(net, data)

print(fitted_net)

## 
##   Bayesian network parameters
## 
##   Parameters of node StudentActivity (multinomial distribution)
## 
## Conditional probability table:
##    Class Library 
##    0.57    0.43 
## 
##   Parameters of node Disappearance (multinomial distribution)
## 
## Conditional probability table:
##  
## , , WitnessTestimony = HeardNothing, ChamberEvidence = SawNothing
## 
##              StudentActivity
## Disappearance     Class   Library
##   NotOccurred 0.3529412 0.4000000
##   Occurred    0.6470588 0.6000000
## 
## , , WitnessTestimony = HeardWhispers, ChamberEvidence = SawNothing
## 
##              StudentActivity
## Disappearance     Class   Library
##   NotOccurred 0.4666667 0.5714286
##   Occurred    0.5333333 0.4285714
## 
## , , WitnessTestimony = HeardNothing, ChamberEvidence = SawSymbols
## 
##              StudentActivity
## Disappearance     Class   Library
##   NotOccurred 0.4285714 0.3846154
##   Occurred    0.5714286 0.6153846
## 
## , , WitnessTestimony = HeardWhispers, ChamberEvidence = SawSymbols
## 
##              StudentActivity
## Disappearance     Class   Library
##   NotOccurred 0.2727273 0.6153846
##   Occurred    0.7272727 0.3846154
## 
## 
##   Parameters of node WitnessTestimony (multinomial distribution)
## 
## Conditional probability table:
##   HeardNothing HeardWhispers 
##          0.54          0.46 
## 
##   Parameters of node ChamberEvidence (multinomial distribution)
## 
## Conditional probability table:
##  SawNothing SawSymbols 
##       0.49       0.51

这段结果是关于一个贝叶斯网络（Bayesian network）的参数描述。贝叶斯网络是一种用于建模随机变量之间的依赖关系的图结构，其中节点表示随机变量，边表示它们之间的依赖关系。每个节点的参数描述了给定其父节点的条件下，该节点的概率分布。

以下是对这段结果的分析：

节点 “StudentActivity” 的参数（多项分布）：
- 条件概率表显示了两个可能的状态 “Class” 和 “Library”，以及它们的概率分布。
- 对于 “Class”，概率为 0.57，对于 “Library”，概率为 0.43。
- 这表示在给定其他节点的情况下，“StudentActivity” 节点有两种可能的状态，“Class” 和 “Library”，并且这些状态的概率分布是给定的。
节点 “Disappearance” 的参数（多项分布）：
- 这个节点的条件概率表根据两个父节点的不同状态进行了分组，父节点是 “WitnessTestimony” 和 “ChamberEvidence”。
- 具体地，有四种可能的组合：HeardNothing & SawNothing，HeardWhispers & SawNothing，HeardNothing & SawSymbols，HeardWhispers & SawSymbols。
- 每种组合下，都有 “Disappearance” 节点的 “NotOccurred” 和 “Occurred” 两种状态，以及它们的概率分布。
节点 “WitnessTestimony” 的参数（多项分布）：
- 这个节点表示 “HeardNothing” 和 “HeardWhispers” 两种可能的状态，以及它们的概率分布。
- “HeardNothing” 的概率为 0.54，“HeardWhispers” 的概率为 0.46。
- 这表示 “WitnessTestimony” 节点的状态取决于这两种可能的状态，并且它们的概率分布是给定的。
节点 “ChamberEvidence” 的参数（多项分布）：
- 这个节点表示 “SawNothing” 和 “SawSymbols” 两种可能的状态，以及它们的概率分布。
- “SawNothing” 的概率为 0.49，“SawSymbols” 的概率为 0.51。
- 这表示 “ChamberEvidence” 节点的状态取决于这两种可能的状态，并且它们的概率分布是给定的。

这些参数描述了一个复杂的贝叶斯网络，可以用于推断在不同条件下各个节点的状态的概率。这个网络的结构和参数可以用于解决特定的概率推断问题，例如根据观察到的证据来估计学生是否消失的概率。

构建逻辑回归模型

我们首先创建和拟合了贝叶斯网络。然后，我们使用逻辑回归模型来预测Disappearance节点的状态。最后，我们定义了一组新的观测数据，并使用predict函数来预测Disappearance发生的概率。

使用predict函数对您的贝叶斯网络进行预测，我们需要先将贝叶斯网络模型转换为一种可以应用预测模型的格式。由于bnlearn包并不直接支持使用predict函数进行预测，我们需要采取替代方法，比如使用逻辑回归。

我们可以基于贝叶斯网络结构来构建一个逻辑回归模型，以预测Disappearance节点的状态。请注意，这种方法实际上并不涉及使用贝叶斯网络进行预测，而是使用网络结构来指导传统的统计模型的建立。

# 构建逻辑回归模型
model <- glm(Disappearance ~ StudentActivity + WitnessTestimony + ChamberEvidence, 
             data = data, family = binomial())

进行预测

# 定义新的观测数据用于预测
new_data <- data.frame(
  StudentActivity = factor("Library", levels = c("Class", "Library")),
  WitnessTestimony = factor("HeardWhispers", levels = c("HeardWhispers", "HeardNothing")),
  ChamberEvidence = factor("SawSymbols", levels = c("SawSymbols", "SawNothing"))
)

# 使用模型进行预测
predicted_probs <- predict(model, newdata = new_data, type = "response")

# 输出预测的概率
print(predicted_probs)

##         1 
## 0.4693946

在这个例子中，首先使用glm函数构建了一个逻辑回归模型，其中Disappearance是因变量，而StudentActivity、WitnessTestimony和ChamberEvidence是自变量。然后，我们定义了一组新的观测数据，并使用predict函数来预测Disappearance发生的概率。经过贝叶斯网络模型对新的观测数据进行了预测，预测的概率为 0.4693946。这个概率表示，在给定新的观测数据的情况下，学生发生消失的概率。

请注意，这种方法并不是直接使用贝叶斯网络进行预测，而是利用网络结构来指导传统统计模型的构建。

MCMC（Markov Chain Monte Carlo）马尔可夫链

MCMC（Markov Chain Monte Carlo）是一种在统计学中使用的算法，用于从复杂的概率分布中抽取样本。在贝叶斯统计中，MCMC特别重要，因为它可以用来近似复杂后验概率分布，特别是当这些分布难以解析计算时。它特别适用于那些难以用解析方法直接计算的复杂概率分布。MCMC通过构建马尔可夫链来生成样本，这些样本的分布逐渐逼近目标概率分布，从而允许我们进行有效的近似推断。

MCMC通过构建一个马尔可夫链来工作。这个链在可能解的空间中随机“游走”，其特点是每一步的位置只依赖于前一步的位置。随着时间的推移，这个链的分布趋于稳定，从而近似目标分布（例如后验分布）。

MCMC的关键概念

马尔可夫链:
- 一种随机过程，当前状态只依赖于前一个状态，不依赖于更早的历史。
蒙特卡罗方法:
- 一种利用随机抽样来近似复杂计算的方法。

MCMC在贝叶斯统计中的应用

后验分布的抽样:
- 在贝叶斯分析中，MCMC用于从后验分布中抽取样本，特别是当这个分布难以直接计算时。
参数估计:
- MCMC方法使我们能够估计模型参数的后验分布，从而了解参数的不确定性。

故事：霍格沃茨火焰杯与马尔可夫链

一年一度三强争霸赛比赛时，哈利波特需要找到特殊的药水，使他能在水下呼吸。这种药水的配方被隐藏起来，只有通过解决一系列复杂的谜题才能找到

谜题的贝叶斯解读

问题: 药水配方由多个成分组成，每个成分的正确比例是未知的。
目标: 使用MCMC来估计这些成分的最佳比例。

MCMC在寻找答案中的角色

构建马尔可夫链: 三人组通过实验和收集线索，逐步构建了一条马尔可夫链，每一步都基于前一步的信息进行调整，逐渐逼近药水的正确配方。
探索可能的解决方案: 通过MCMC，他们能够探索一系列可能的配方，每个配方都有一定的概率是正确的
找到最可能的配方: 随着时间的推移，他们发现某些配方比其他配方更可能是正确的。MCMC帮助他们确定哪些成分比例最有可能产生有效的药水。

药水的成功制作

最终，通过MCMC的帮助，赫敏成功地锁定了最有可能的药水配方，并让哈利波特在水下任务中成功救出芙蓉的妹妹和罗恩。

MCMC的实际应用

在现实世界中，MCMC被广泛用于估计复杂模型的参数，特别是在贝叶斯统计学中。它允许研究人员从概率模型的后验分布中抽取样本，从而进行参数估计和预测。虽然MCMC是一种强大的工具，但它也需要仔细的设计和验证，以确保得到的样本是代表性的并且收敛于真实的后验分布。

在实际应用中，MCMC被用于各种统计模型和问题中，例如在经济学、生物统计学和机器学习等领域。它特别适合于那些后验分布难以直接计算，或者模型过于复杂无法解析求解的情况。

MCMC是一种强大但计算密集的方法，它要求仔细地选择和调整算法参数，比如步长、迭代次数等，以确保有效和准确的抽样。

R 语言中的MCMC实现

在R中，有多个包可用于实现MCMC，如rjags、Stan和BUGS等。这些包提供了工具和函数来定义模型、执行MCMC抽样并分析结果。

我使用R语言实现MCMC的完整示例。在这个例子中，我们将使用rjags包，它是一个接口到JAGS（Just Another Gibbs Sampler）的库，JAGS是一个用于贝叶斯模型的MCMC抽样的程序。

注意！！MAC上无法install的小伙伴看过来：

Easiest way to install native arm64 jags: brew install pkg-config jags /Applications/RStudio.app/Contents/MacOS/RStudio install.packages(“rjags”, type=“source”)

JAGS简介

JAGS是一个独立的程序，用于执行MCMC模拟。它接受贝叶斯模型描述（先验加似然）和数据作为输入，返回后验分布的MCMC样本。JAGS使用Metropolis采样、Gibbs采样和其他MCMC算法的组合。

加载数据
- 数据可以从文件加载，或通过足够的摘要统计信息指定。例如，可以指定样本大小和成功次数。
指定模型：似然和先验
- JAGS模型规范以model开始。模型提供了似然和先验的文本描述。这个文本字符串随后会传递给JAGS进行翻译。
- 以Beta-Binomial模型为例，先验分布是Beta分布，似然函数是二项分布。
在JAGS中编译模型
- 将模型（即文本字符串）和数据传递给JAGS进行编译。可以通过textConnection函数定义文本字符串。模型也可以保存在单独的文件中，文件名被传递给JAGS。数据以列表形式传递给JAGS。
从后验分布中模拟值
- 在JAGS中模拟值分为两个步骤。首先，使用update命令运行模拟一个“烧入”期。然后，使用coda.samples从后验分布中模拟我们实际保留的值。使用coda.samples可以方便地使用coda包来总结和诊断MCMC模拟。
总结模拟值和诊断检查
- 使用R的标准函数（如summary和plot）来总结coda.samples的结果。可以对theta的模拟值进行总结，以近似后验分布。

if (!require("rjags")) {
    install.packages("rjags", type="source")
}

## Loading required package: rjags

## Loading required package: coda

## Linked to JAGS 4.3.2

## Loaded modules: basemod,bugs

library(rjags)

n = 35 # sample size
y = 31 # number of successes

# 编写模型
model_string <- "model{

  # Likelihood
  y ~ dbinom(theta, n)

  # Prior
  theta ~ dbeta(alpha, beta)
  alpha <- 3 # prior successes
  beta <- 1 # prior failures

}"

dataList = list(y = y, n = n)

model <- jags.model(file = textConnection(model_string), 
                    data = dataList)

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 1
##    Unobserved stochastic nodes: 1
##    Total graph size: 5
## 
## Initializing model

# 从后验分布中模拟值
update(model, n.iter = 1000) # 模拟的值的数量

Nrep = 10000 # number of values to simulate

posterior_sample <- coda.samples(model,
                       variable.names = c("theta"),
                       n.iter = Nrep)

summary(posterior_sample)

## 
## Iterations = 2001:12000
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##           Mean             SD       Naive SE Time-series SE 
##      0.8717344      0.0522225      0.0005222      0.0007031 
## 
## 2. Quantiles for each variable:
## 
##   2.5%    25%    50%    75%  97.5% 
## 0.7532 0.8387 0.8779 0.9104 0.9544

plot(posterior_sample)

thetas = as.matrix(posterior_sample)
head(thetas)

##          theta
## [1,] 0.8904315
## [2,] 0.9017134
## [3,] 0.8613904
## [4,] 0.8596023
## [5,] 0.9038987
## [6,] 0.8702215

hist(thetas)

ynew = rbinom(Nrep, n, thetas)

plot(table(ynew),
     main = "Posterior Predictive Distribution for samples of size 35",
     xlab = "y")

The JAGS (Just Another Gibbs Sampler) is a program used for performing Markov Chain Monte Carlo (MCMC) simulations in Bayesian statistics. JAGS is particularly useful for Bayesian model estimation, where it takes a model description (prior and likelihood) and data to return an MCMC sample from the posterior distribution. The process involves various MCMC algorithms, including Metropolis and Gibbs sampling.

Key steps in using JAGS are:

Load the Data: Data can be loaded from a file or specified via summary statistics. An example used is a “data is singular” context with a specified sample size and number of successes.
Specify the Model: This involves writing a model statement that describes the likelihood and prior distributions. For example, in a Beta-Binomial model, you specify the Beta distribution as the prior and Binomial as the likelihood function.
Compile the Model in JAGS: The model, written as a text string, is passed to JAGS along with the data for compilation.
Simulate Values from the Posterior Distribution: After a ‘burn-in’ period using the update function, values from the posterior distribution are simulated using coda.samples. This step is crucial for generating the posterior samples needed for inference.
Summarize Simulated Values and Check Diagnostics: Standard functions from R (like summary and plot) are used to summarize and visually inspect the results from the MCMC simulation. This step is vital for understanding the posterior distribution and checking the convergence and other diagnostic measures of the MCMC simulation.

The section emphasizes the ease of integrating JAGS with R using the rjags package, allowing for a seamless workflow within the R environment. This integration is particularly beneficial for those who are already familiar with R and wish to conduct Bayesian data analysis without switching between different software environments.

贝叶斯层次模型（Bayesian Hierarchical Models）简介

贝叶斯层次模型，也称为多级模型，是一种统计模型，它允许数据在不同层次上进行建模，从而更好地处理复杂的数据结构。在这种模型中，参数本身也可以是随机变量，其分布可以依赖于更高层次的参数。

关键特点

多层次结构:
- 数据可以在不同层次上进行建模，例如个体层次和群体层次。
灵活性:
- 能够适应复杂的数据结构，如嵌套数据或分组数据。
共享信息:
- 不同组或个体之间的信息可以共享，提高估计的准确性。

故事：霍格沃茨魔法考试成绩分析

在霍格沃茨，每年都会举行魔法考试。考试成绩受多种因素影响，比如学生的勤奋程度、所属学院以及教授的教学方式。

场景设置

数据层次:
- 学生层次：每个学生的勤奋程度和考试成绩。
- 学院层次：每个学院的教学资源和声誉。
- 教授层次：教授的教学风格和经验。

贝叶斯层次模型的应用

分析目标:
- 分析学生的魔法考试成绩，了解不同因素如何影响成绩。
模型构建:
- 在学生层次，模型考虑个人勤奋程度对成绩的影响。
- 在学院层次，模型考虑学院资源对学生成绩的影响。
- 在教授层次，模型考虑教授教学风格的影响。
共享信息:
- 学生的成绩信息可以帮助了解学院和教授层次的效果。
结果应用:
- 分析结果可以用来改善教学方法，提高学院的教学资源配置。

当然，我可以提供一个使用R语言来实现贝叶斯层次模型的简单示例。在这个例子中，我们将使用rstan包，这是一个R接口到Stan，后者是一个用于贝叶斯统计推断的C++库。

R 语言示例场景

假设我们有一组学生的考试成绩数据，这些学生来自不同的班级。我们的目标是估计每个班级的平均成绩，同时考虑到班级间的差异。

安装和加载必要的包

if (!require("rstan")) {
    install.packages("rstan")
}

## Loading required package: rstan

## Loading required package: StanHeaders

## 
## rstan version 2.32.3 (Stan version 2.26.1)

## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
## For within-chain threading using `reduce_sum()` or `map_rect()` Stan functions,
## change `threads_per_chain` option:
## rstan_options(threads_per_chain = 1)

## 
## Attaching package: 'rstan'

## The following object is masked from 'package:coda':
## 
##     traceplot

library(rstan)

准备数据

# 模拟数据
set.seed(123)
n_classes <- 4
n_students <- 100

# 随机分配学生到班级
class_id <- sample(1:n_classes, n_students, replace = TRUE)

# 生成班级的真实平均成绩
class_means <- rnorm(n_classes, mean = 70, sd = 5)

# 生成学生的成绩
grades <- class_means[class_id] + rnorm(n_students, mean = 0, sd = 10)

# 数据列表
data_list <- list(J = n_classes,
                  N = n_students,
                  class_id = class_id,
                  y = grades)

定义贝叶斯层次模型

model_string <- "
data {
  int<lower=0> J;                 // 班级数量
  int<lower=0> N;                 // 学生数量
  int<lower=1, upper=J> class_id[N]; // 班级索引
  vector[N] y;                    // 成绩
}

parameters {
  vector[J] mu;                   // 班级平均成绩
  real<lower=0> sigma;            // 班级间的标准差
  real<lower=0> sigma_y;          // 学生成绩的标准差
}

model {
  sigma ~ uniform(0, 30);
  sigma_y ~ uniform(0, 30);
  mu ~ normal(70, sigma);         // 班级平均成绩的先验分布
  for (n in 1:N)
    y[n] ~ normal(mu[class_id[n]], sigma_y); // 学生成绩的分布
}
"

运行模型

# 编译Stan模型
stan_model <- stan_model(model_code = model_string)

## Trying to compile a simple C file

## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 15.0.0 (clang-1500.0.40.1)’
## using SDK: ‘’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
##                ^
##                ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
##          ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1

# 运行MCMC
fit <- sampling(stan_model, data = data_list, iter = 2000, chains = 4)

## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 2.4e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.24 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.056 seconds (Warm-up)
## Chain 1:                0.042 seconds (Sampling)
## Chain 1:                0.098 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 6e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.067 seconds (Warm-up)
## Chain 2:                0.039 seconds (Sampling)
## Chain 2:                0.106 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 7e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 0.064 seconds (Warm-up)
## Chain 3:                0.027 seconds (Sampling)
## Chain 3:                0.091 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 6e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 0.058 seconds (Warm-up)
## Chain 4:                0.078 seconds (Sampling)
## Chain 4:                0.136 seconds (Total)
## Chain 4:

## Warning: There were 3996 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.

## Warning: Examine the pairs() plot to diagnose sampling problems

## Warning: The largest R-hat is 3.96, indicating chains have not mixed.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#r-hat

## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess

## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-ess

分析结果

# 摘要统计
print(fit)

## Inference for Stan model: anon_model.
## 4 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##            mean se_mean   sd    2.5%     25%     50%     75%   97.5% n_eff
## mu[1]     -0.71    0.70 1.00   -1.71   -1.27   -1.06   -0.52    1.10     2
## mu[2]      0.24    0.76 1.08   -1.66   -0.20    0.58    1.02    1.40     2
## mu[3]     -0.66    0.53 0.76   -1.60   -1.33   -0.82    0.06    0.54     2
## mu[4]     -0.60    0.67 0.95   -1.72   -1.17   -0.86   -0.20    0.95     2
## sigma     28.88    0.07 0.93   26.47   28.39   29.11   29.63   29.96   188
## sigma_y   29.93    0.00 0.07   29.76   29.91   29.95   29.98   30.00   206
## lp__    -647.38    3.15 4.62 -654.08 -651.15 -649.27 -644.18 -639.72     2
##          Rhat
## mu[1]   32.12
## mu[2]   19.39
## mu[3]   13.41
## mu[4]   24.21
## sigma    1.03
## sigma_y  1.03
## lp__     3.82
## 
## Samples were drawn using NUTS(diag_e) at Mon Dec 11 11:28:34 2023.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

# 提取和查看参数
posterior_samples <- extract(fit)
print(posterior_samples$mu[3])

## [1] -0.9727095

使用诊断工具（如traceplot、pairs等）来诊断问题并调试模型

traceplot(fit, pars = c("mu", "sigma", "sigma_y"))

在traceplot中，你希望看到每个参数的稳定和相对平坦的线条，这表明了良好的混合和收敛。如果traceplot显示出非随机的行为（如趋势或循环），这可能是一个问题的迹象。

这张图是MCMC的追踪图（traceplot），展示了Stan模型运行后每个参数的迭代值。每行表示一个参数，每个参数对应四条线，每条线代表一个MCMC链。理想情况下，这些线应当比较平坦且在一个恒定的水平范围内波动，表明链已经达到稳态并且在后验分布中进行有效的随机游走。现在让我们逐一解释这些追踪图：

mu[1]到mu[4]:
- 这些是班级平均成绩的追踪图。理想情况下，每个班级的mu值在每条链中都应稳定在一个特定水平，且四条链之间应该是相互重叠的，表明链间一致性。
- 在您的图中，mu[1]和mu[2]看起来有稳定的追踪，但是mu[3]和mu[4]在某些链中有上升或下降的趋势，这可能表示链还没有完全达到稳态。
sigma:
- 这个参数表示班级间平均成绩的标准差。我们希望看到的是所有链在某个水平附近波动，没有明显的趋势。
- 在您的图中，sigma的四条链看起来非常混乱且有很大的波动，这表明参数可能没有很好的收敛。
sigma_y:
- 这个参数代表学生成绩的标准差。理想的追踪图会显示出稳定的波动，没有趋势。
- 在您的图中，sigma_y的追踪图显示出所有链都在一个非常狭窄的范围内波动，这看起来似乎收敛良好。

诊断

如果追踪图中的线条呈现出明显的非随机趋势或没有稳定在一个常数水平，这通常是缺乏收敛的迹象。
如果不同的链在不同的水平上波动，这可能表明链之间存在不一致，也是缺乏收敛的标志。
如果所有链看起来很像，并且在一定范围内随机波动，则可能表明参数有好的收敛性。

调整模型

根据追踪图的显示结果，您可能需要更多的迭代次数，尤其是在预热（warmup）阶段，以便链达到稳定状态。
您可能还需要调整采样参数，如步长（step size）或调整模型本身以改善收敛性。
重新审视先验分布的选择也可能有助于改善模型的行为。

追踪图是了解MCMC采样过程的非常重要的工具，因为它们直观地显示了参数值随迭代过程的变化情况。在您的模型中，尤其是mu[3]、mu[4]和sigma似乎有一些潜在的收敛问题。这可能需要进一步的调查和潜在的模型调整。

# 使用pairs函数
pairs(fit, pars = c("mu", "sigma", "sigma_y"))

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

## Warning in par(usr): argument 1 does not name a graphical parameter

在pairs图中，你希望看到参数之间没有明显的奇怪模式或强相关性。如果pairs图显示出一些参数之间有强烈的相关性，可能需要重新考虑模型的参数化。

这张图是Stan模型结果的散点图矩阵，也称为pairs图。它展示了模型参数的双变量关系和单变量边缘分布。对于每个参数对，散点图描绘了一个参数相对于另一个的值，而直方图则显示了单个参数的边缘分布。这些图有助于识别参数间的潜在相关性和模型的行为。让我们逐一解释这些图：

对角线上的直方图:
- 对角线上的图显示的是各个参数的边缘后验分布。这有助于我们了解每个参数的不确定性和可能的值范围。理想情况下，我们希望看到的是类似钟形的正态分布，表示参数估计的稳定性。
散点图:
- 非对角线上的图显示的是两个参数间的散点图。蓝色点表示散点的密度，红色点表示实际样本。这些图可以揭示参数之间的潜在相关性。没有明显模式的随机分布表明参数间相互独立，而清晰的模式或趋势则可能表示参数间存在相关性。
参数之间的关系:
- 对于mu参数，我们希望看到参数之间没有明显的相关性，这意味着班级平均成绩的估计是独立的。
- sigma和sigma_y之间的图显示了它们的联合分布。由于sigma代表班级间的差异，而sigma_y代表个体差异，我们希望看到这两个参数之间没有强烈的相关性。

图中的问题

mu[3]和mu[4]的边缘分布看起来偏离了正态分布，表明可能存在模型设定的问题。
sigma的边缘分布非常奇异，表明了潜在的模型问题，可能是由于先验分布的选择不当，或者模型没有很好地捕捉数据的结构。
散点图中出现的非随机模式可能表明参数之间存在意料之外的相关性，这需要进一步调查。

调整模型

需要重新审视模型的结构和先验设置，特别是针对sigma的先验。
检查数据是否有异常值或数据处理的问题，这可能会影响模型参数的估计。
如果模型仍然表现不佳，可能需要考虑使用不同的模型或统计方法。

在使用这些诊断工具时，重要的是要结合参数的实际意义和模型的上下文来解释结果，并做出相应的调整。在贝叶斯建模中，可能需要多次迭代来找到合适的模型和参数设置。

模型组成部分解释

这个例子中的模型考虑了班级间的平均成绩差异（mu），以及学生成绩的个体差异（sigma_y）。通过贝叶斯层次模型，我们能够同时估计班级层面和个体层面的参数

这个模型是一个贝叶斯层次线性模型，用于分析学生的成绩数据，这些数据被分组到不同的班级中。模型的目标是估计每个班级的平均成绩，同时考虑到班级间的差异以及学生间的个体差异。

数据部分:
- int<lower=0> J;：定义了一个整数变量J，表示班级的总数，其值必须非负。
- int<lower=0> N;：定义了一个整数变量N，表示学生的总数，其值也必须非负。
- int<lower=1, upper=J> class_id[N];：这是一个整数数组，每个元素对应一个学生，值介于1和J之间，表示学生所在的班级。
- vector[N] y;：一个实数向量，包含所有学生的成绩。
参数部分:
- vector[J] mu;：一个实数向量，长度为J，表示每个班级的平均成绩。
- real<lower=0> sigma;：一个非负实数，表示班级平均成绩之间的标准差，用于描述班级间成绩的变异性。
- real<lower=0> sigma_y;：一个非负实数，表示每个学生成绩的标准差，用于描述个体学生成绩的变异性。
模型部分:
- sigma ~ uniform(0, 30);：为sigma参数指定了一个0到30的均匀先验分布。
- sigma_y ~ uniform(0, 30);：为sigma_y参数指定了一个0到30的均匀先验分布。
- mu ~ normal(70, sigma);：为每个班级的平均成绩mu指定了一个以70为均值，sigma为标准差的正态先验分布。
- 对于每个学生n，y[n] ~ normal(mu[class_id[n]], sigma_y);：这表示每个学生的成绩y[n]服从以其班级平均成绩mu[class_id[n]]为均值，以sigma_y为标准差的正态分布。

模型的目的和应用

这个模型的主要目的是估计不同班级的平均成绩，并了解班级间成绩的差异。通过对mu的估计，可以了解哪些班级的表现更好或更差；而通过sigma和sigma_y的估计，可以了解班级间成绩差异的大小以及个体成绩的波动范围。这种模型在教育研究中非常有用，尤其是在评估和比较不同班级或教学方法的效果时。

模型结果

是Stan模型的结果摘要包含了模型参数的后验分布统计摘要，包括每个参数的均值（mean）、标准误（se_mean）、标准差（sd）、以及不同百分位数的估计值（如2.5%、25%、50%、75%和97.5%）。此外，还包括了每个参数的有效样本大小（n_eff）和收敛诊断统计量Rhat。

分析模型结果的关键点

参数估计:
- 查看mu[1]到mu[4]的均值，这些是不同班级平均成绩的估计。
- sigma的均值表示班级间平均成绩的变异性。
- sigma_y的均值表示个体学生成绩的标准差。
不确定性估计:
- 查看每个参数的置信区间（2.5%至97.5%），了解参数估计的不确定性。
- 参数的标准差（sd）也提供了关于估计不确定性的信息。
模型诊断:
- n_eff是有效样本大小，值越大表示估计越稳定。
- Rhat是收敛诊断指标，理想值为1。值远大于1可能表明模型未收敛。

注意事项

有效样本大小 (n_eff):
- 您的结果中mu的有效样本大小非常低（例如，mu[1]的n_eff为2），这表明抽样效率可能很低。
收敛问题 (Rhat):
- mu的Rhat值远大于1（例如，mu[1]的Rhat为31.55），这是模型未收敛的强烈迹象。通常，Rhat值应接近1。
可能的问题:
- 低n_eff和高Rhat值表明模型可能有问题。可能是由于先验设置不当、模型过于复杂或数据不足等原因。

建议

增加迭代次数:
- 增加MCMC的迭代次数可能有助于改善有效样本大小和收敛。
调整模型或先验:
- 重新审视模型结构和先验分布，确保它们适合您的数据。
诊断和调试:
- 使用诊断工具（如traceplot、pairs等）来诊断问题并调试模型。

贝叶斯层次模型的实际应用

在实际应用中，贝叶斯层次模型广泛应用于教育、生态学、医学等领域，尤其是在数据具有自然分层结构时。它们允许研究人员对复杂的数据进行深入分析，同时考虑数据中的不确定性和变异性。

Bayesian Statistics

贝叶斯统计（Bayesian Statistics）定义

故事：霍格沃茨的魔法石猜想

教授的先验信念

新证据出现

更新信念：贝叶斯更新

结论

实践中的应用

贝叶斯统计进阶

高级贝叶斯方法

应用领域

工具和技术

贝叶斯网络（Bayesian Networks）定义

贝叶斯网络的关键特点

贝叶斯与霍格沃茨密室

背景

构建贝叶斯网络

收集证据

使用贝叶斯网络进行推理

R语言实现

贝叶斯网络结果分析

构建逻辑回归模型

进行预测

MCMC（Markov Chain Monte Carlo）马尔可夫链

MCMC的关键概念

MCMC在贝叶斯统计中的应用

故事： 霍格沃茨火焰杯与马尔可夫链

MCMC的实际应用

R 语言中的MCMC实现

JAGS简介

贝叶斯层次模型（Bayesian Hierarchical Models）简介

关键特点

故事：霍格沃茨魔法考试成绩分析

场景设置

贝叶斯层次模型的应用

R 语言示例场景

安装和加载必要的包

准备数据

定义贝叶斯层次模型

运行模型

分析结果

使用诊断工具（如traceplot、pairs等）来诊断问题并调试模型

诊断

调整模型

图中的问题

调整模型

模型组成部分解释

模型的目的和应用

模型结果

分析模型结果的关键点

注意事项

建议

贝叶斯层次模型的实际应用

故事：霍格沃茨火焰杯与马尔可夫链