这是一个很好的问题，似乎经常出现：链接 1，链接 2。Cam.Davidson.Pilon指出的 Bayesian Estimation Superseeds the T-Test论文是该主题的绝佳资源。它也是最近的，发表于 2012 年，我认为部分原因是当前对该领域的兴趣。

我将尝试总结贝叶斯替代两个样本 t 检验的数学解释。这个总结类似于 BEST 论文，它通过比较两个样本的后验分布的差异来评估两个样本的差异（在下面的 R 中解释）。

set.seed(7)

#create samples
sample.1 <- rnorm(8, 100, 3)
sample.2 <- rnorm(10, 103, 7)

#we need a pooled data set for estimating parameters in the prior.
pooled <- c(sample.1, sample.2)
par(mfrow=c(1, 2))

hist(sample.1)
hist(sample.2)

为了比较样本意味着我们需要估计它们是什么。这样做的贝叶斯方法使用贝叶斯定理：P(A|B) = P(B|A) * P(A)/P(B)（P(A|B) 的语法被解读为A 给定 B)

感谢现代数值方法，我们可以忽略 B 的概率，P(B)，并使用比例语句：P(A|B)$\propto$P(B|A)*P(A) 在贝叶斯白话 中，后验与似然乘以先验成正比

将贝叶斯理论应用于我们的问题，在给定一些数据的情况下，我们想知道样本的均值$P(mean.1 | sample.1)$ $\propto$ $P(sample.1 | mean.1) * P(mean.1)$. 右边的第一项是可能性，$P(sample.1 | mean.1)$，即在给定均值的情况下观察样本数据的概率。 1．第二项是前项，$P(mean.1)$, 这只是 mean.1 的概率。找出合适的先验仍然是一门艺术，也是对贝叶斯方法的最大批评之一。

让我们把它放在代码中。代码让一切变得更好。

likelihood <- function(parameters){
  mu1=parameters[1]; sig1=parameters[2]; mu2=parameters[3]; sig2=parameters[4]
  prod(dnorm(sample.1, mu1, sig1)) * prod(dnorm(sample.2, mu2, sig2))
}

prior <- function(parameters){
  mu1=parameters[1]; sig1=parameters[2]; mu2=parameters[3]; sig2=parameters[4]
  dnorm(mu1, mean(pooled), 1000*sd(pooled)) * dnorm(mu2, mean(pooled), 1000*sd(pooled)) * dexp(sig1, rate=0.1) * dexp(sig2, 0.1)
}

我在前面做了一些需要证明的假设。为了防止先验影响估计的平均值，我想让它们在合理的值上变得广泛和统一，目的是让数据产生后验的特征。我使用了 BEST 推荐的设置，并以均值 = 均值（合并）和广泛的标准差 = 1000*sd（合并）来正常分配 mu。我将标准差设置为广泛的指数分布，因为我想要一个广泛的无界分布。

现在我们可以使后

posterior <- function(parameters) {likelihood(parameters) * prior(parameters)}

我们将使用带有 Metropolis Hastings 修正的马尔可夫链蒙特卡罗(MCMC) 对后验分布进行采样。用代码最容易理解。

#starting values
mu1 = 100; sig1 = 10; mu2 = 100; sig2 = 10
parameters <- c(mu1, sig1, mu2, sig2)

#this is the MCMC /w Metropolis method
n.iter <- 10000
results <- matrix(0, nrow=n.iter, ncol=4)
results[1, ] <- parameters
for (iteration in 2:n.iter){
  candidate <- parameters + rnorm(4, sd=0.5)
  ratio <- posterior(candidate)/posterior(parameters)
  if (runif(1) < ratio) parameters <- candidate #Metropolis modification
  results[iteration, ] <- parameters
}

结果矩阵是来自每个参数的后验分布的样本列表，我们可以用它来回答我们最初的问题：sample.1 与 sample.2 不同吗？但首先为了避免初始值的影响，我们将“烧入”链的前 500 个值。

#burn-in
results <- results[500:n.iter,]

现在，sample.1 与 sample.2 不同吗？

mu1 <- results[,1]
mu2 <- results[,3]

hist(mu1 - mu2)

mean(mu1 - mu2 < 0)
[1] 0.9953689

从这个分析中，我得出结论，sample.1 的平均值小于 sample.2 的平均值的可能性为 99.5%。

正如 BEST 论文中所指出的，贝叶斯方法的一个优点是它可以提出强有力的理论。EG sample.2 比 sample.1 大 5 个单位的概率是多少。

mean(mu2 - mu1 > 5)
[1] 0.9321124

我们可以得出结论，sample.2 的平均值比 sample.1 大 5 个单位的可能性为 93%。细心的读者会发现这很有趣，因为我们知道真实总体的平均值分别为 100 和 103。这很可能是由于样本量小，以及选择使用正态分布的可能性。

我将以警告结束这个答案：此代码用于教学目的。对于实际分析，请使用 RJAGS 并根据您的样本量拟合可能性的 t 分布。如果有兴趣，我将使用 RJAGS 发布 t 检验。

编辑：这里要求的是一个 JAGS 模型。

model.str <- 'model {
    for (i in 1:Ntotal) {
        y[i] ~ dt(mu[x[i]], tau[x[i]], nu)
    }
    for (j in 1:2) {
        mu[j] ~ dnorm(mu_pooled, tau_pooled)
        tau[j] <- 1 / pow(sigma[j], 2)
        sigma[j] ~ dunif(sigma_low, sigma_high)
    }
    nu <- nu_minus_one + 1
    nu_minus_one ~ dexp(1 / 29)
}'

# Indicator variable
x <- c(rep(1, length(sample.1)), rep(2, length(sample.2)))

cpd.model <- jags.model(textConnection(model.str),
                        data=list(y=pooled,
                                  x=x,
                                  mu_pooled=mean(pooled),
                                  tau_pooled=1/(1000 * sd(pooled))^2,
                                  sigma_low=sd(pooled) / 1000,
                                  sigma_high=sd(pooled) * 1000,
                                  Ntotal=length(pooled)))
update(cpd.model, 1000)
chain <- coda.samples(model = cpd.model, n.iter = 100000,
                      variable.names = c('mu', 'sigma'))
rchain <- as.matrix(chain)
hist(rchain[, 'mu[1]'] - rchain[, 'mu[2]'])
mean(rchain[, 'mu[1]'] - rchain[, 'mu[2]'] < 0)
mean(rchain[, 'mu[2]'] - rchain[, 'mu[1]'] > 5)

用 Python 实现的 user1068430 的出色回答

import numpy as np
from pylab import plt

def dnorm(x, mu, sig):
    return 1/(sig * np.sqrt(2 * np.pi)) * np.exp(-(x - mu)**2 / (2 * sig**2))

def dexp(x, l):
    return l * np.exp(- l*x)

def like(parameters):
    [mu1, sig1, mu2, sig2] = parameters
    return dnorm(sample1, mu1, sig1).prod()*dnorm(sample2, mu2, sig2).prod()

def prior(parameters):
    [mu1, sig1, mu2, sig2] = parameters
    return dnorm(mu1, pooled.mean(), 1000*pooled.std()) * dnorm(mu2, pooled.mean(), 1000*pooled.std()) * dexp(sig1, 0.1) * dexp(sig2, 0.1)

def posterior(parameters):
    [mu1, sig1, mu2, sig2] = parameters
    return like([mu1, sig1, mu2, sig2])*prior([mu1, sig1, mu2, sig2])


#create samples
sample1 = np.random.normal(100, 3, 8)
sample2 = np.random.normal(100, 7, 10)

pooled= np.append(sample1, sample2)

plt.figure(0)
plt.hist(sample1)
plt.hold(True)
plt.hist(sample2)
plt.show(block=False)

mu1 = 100 
sig1 = 10
mu2 = 100
sig2 = 10
parameters = np.array([mu1, sig1, mu2, sig2])

niter = 10000

results = np.zeros([niter, 4])
results[1,:] = parameters

for iteration in np.arange(2,niter):
    candidate = parameters + np.random.normal(0,0.5,4)
    ratio = posterior(candidate)/posterior(parameters)
    if np.random.uniform() < ratio:
        parameters = candidate
    results[iteration,:] = parameters

#burn-in
results = results[499:niter-1,:]

mu1 = results[:,1]
mu2 = results[:,3]

d = (mu1 - mu2)
p_value = np.mean(d > 0)

plt.figure(1)
plt.hist(d,normed = 1)
plt.show()

通过贝叶斯分析，您可以指定更多的东西（这实际上是一件好事，因为它提供了更多的灵活性和能力来模拟您认为的真相）。您是否假设可能性的正常值？两组的方差是否相同？

一种直接的方法是对 2 个均值（和 1 个或 2 个方差/离差）进行建模，然后查看 2 个均值差异的后验和/或 2 个均值差异的可信区间。

我可以使用哪些贝叶斯方法来测试两个样本的平均值之间的差异的数学解释。

有几种方法可以“测试”这一点。我会提到一对：

• 如果你想要一个明确的决定，你可以看看决策理论。

• 有时要做的一件非常简单的事情是找到均值差异的区间，并考虑它是否包含 0。这将涉及从观察模型开始，先验参数和计算以数据为条件的均值差异的后验分布。

  您需要说出您的模型是什么（例如，正常的、恒定的方差），然后（至少）一些关于均值差异的先验和方差的先验。您可能依次对这些先验的参数有先验。或者你可能不会假设恒定的方差。或者你可能会假设正常以外的其他东西。