【趣学】贝叶斯统计数据如何说服我去健身房- 专知

会员服务 ·

【趣学】贝叶斯统计数据如何说服我去健身房

专知出品

【导读】你是否真的了解贝叶斯方法？本文通过简单有趣的案例为大家详解如何利用贝叶斯方法学习线性回归理论。代码附在文章最后，试试看get一份关于你自己的身体素质报告，让数据告诉你，少年你是否该去健身房了。

作者 | Tuan Doan Nguyen

编译 | 专知

整理 | Xiaowen

How Bayesian statistics convinced me to hit the gym？

这是一段有趣的旅程，我们用贝叶斯的方法学习线性回归理论。

正如这篇文章的标题诉述的，我将对我的身体进行科学的调查。

“你是认真的吗……”

是的。所以我有一点背景资料需要介绍：我最初来自越南，在新加坡上高中，目前在美国上大学。我经常听到人们取笑我看起来有多弱小，我体重有多轻，我应该如何锻炼，去健身房增加体重，以便拥有“更好的体格”……我经常听到这些评论。对身高1m69(5‘6)，体重58公斤(127磅)的人来说，我几乎有接近完的BMI指数(20.3)。

考虑到这一点，人们可能有一个观点：考虑到我比一个越南男人高一点，体重和一般的越南男人完全一样(维基百科的简介是1m62和58公斤)，我可能会“看起来”有点瘦。“看”这是关键词。逻辑很简单，不是吗？如果你伸展一下身体，有同样的体重，你应该看上去更瘦。我认为这是一个很严肃的科学问题，需要进一步研究：

“对于一个身高1m69的越南男性来说，我有多小(体重)？”

我们需要一种方法来研究这个问题。一个很好的方法是尽可能多地找到越南男子身高和体重的数据，看看我在他们中的位置。

越南人口概况

在浏览网络后，我发现了一项调查研究数据，其中包含了1万多名越南人的人口统计资料。我把抽样范围限制在18-29岁年龄组的男性，这样我得到了383名年龄在18-29岁左右的越南男性的样本，这应该足以进行分析。

首先，让我们想象一下人口的重量，看看我在年轻的越南男人中的表现。

红线显示样本的中数，橙色线显示平均值。

这个图表明我略低于那383名越南年轻人的平均体重和中位体重。令人鼓舞的消息？不过，问题不在于我的体重与样本本身的体重比较如何，但考虑到1m68的高度，假设越南男性人口是由这383人所代表的，我们可以推断出我的体重与整个越南人口相比，是怎样的呢？要做到这一点，我们需要进行回归分析。

首先，让我们绘制一个高度和权重的二维散点图。

嗯，我只是看平均。事实上，如果我们只看168厘米高的人的数据(想象一下168厘米的垂直线，穿过红点)，我的体重就比这些人少了一点。

另一项重要的观察是，图中越南男性的身高与体重呈正线性关系，我们将做定量分析，以找出这一关系的根源。

首先，让我们快速添加一个“标准最小二乘”线。稍后我会说更多关于这行的内容，现在就先展示出来。

我们可以把我们的最小二乘线y=−86.32+0.889x解释为：在我这个年龄，越南男性的一般轮廓是，身高增加1厘米，体重就会增加0.88公斤。

然而，这并没有回答我们的问题，那就是，考虑到我的身高(1m68)，58公斤的体重被认为过高、太低或几乎平均吗？用更定量的方式重新表述这个问题，如果我们有1m68身高的人的分布，那么我体重下降到25%，50%或75%的几率有多大呢？为了做到这一点，让我们来挖掘一下。深入了解回归背后的理论。

线性回归理论

在线性回归模型中，Y变量的期望值(在我们的例子中是人的体重)是x(高度)的线性函数，让我们调用这个线性关系β0和β1的截距和斜率，即我们假定E(Y|X=x)=β0+β1×x，我们不知道β0和β1，所以我们把它们当作未知的参数。

在最标准的线性回归模型中，我们进一步假定Y的条件分布X=x是正态分布。

这意味着简单的线性回归模型：

可以重写为：

（1）

值得注意的是，在许多模型中，我们可以用精确参数τ代替方差参数σ，其中τ=1/σ。

总结如下：因变量Y服从正态分布，其参数为均值μi，精度参数τ。μi为X参数的线性函数，由β0和β1进行参数化。

最后，我们还假设未知方差不依赖于x；这个假设称为同方差性（homoscedasticity）。

在一个实际的数据分析问题中，我们得到的只是黑点-数据。利用这些数据，我们的目标是对我们不知道的事物进行推断，包括β0、β1(确定阴影、蓝色虚线)和σ(确定绘制y数据值的红色正态密度的宽度)。每个点周围的分布看起来完全相同。这是同方差性（homoscedasticity）的本质。

参数估计

现在，有几种方法可以用来估计β0和β1。如果你使用最小二乘法估计这种模型，就不必担心概率公式，因为你正在通过最小化拟合值的平方误差来寻找β0和β1的最优值。另一方面，你可以使用极大似然估计来估计这类模型，其中你将通过最大化似然函数来寻找参数的最优值。

注意：一个有趣的结果(这里没有证明)是，如果我们进一步假设误差也属于正态分布，最小二乘估计也是最大似然估计。

利用贝叶斯方法进行线性回归

贝叶斯方法不只是最大化似然函数，而是假定参数的先验分布，并使用Bayes定理：

似然函数与上面的相同，不同的是假设参数β0，β1，τ的先验分布，并将它们包含在方程中：

（2）

“等等，这是什么？为什么它会使我们的公式看起来更复杂？”

相信我，先验概率感觉有点奇怪，但它很直观。事实是，对于为什么我们可以用一些看似任意的分布来分配一个未知参数(在我们的例子中是β0、β1、τ)，有一个非常强烈的哲学推理。这些先验分布是为了在看到数据之前捕捉我们对情况的信念（captureour beliefs about the situation before seeing the data）。在观察了一些数据之后，我们应用Bayes准则得到了这些未知数的后验分布，它既考虑了先验分布，又考虑了数据，由此可以计算出未来观测的预测分布。

最终的估计将取决于1)你的数据和2)你的先验信息，但是你的数据中包含的信息越多，先验的影响力就越小。

“这样我就可以选择任何先验分布？”

这是一个很好的问题，因为选择的数量是无限的。(理论上)只有一种正确的先验，即捕捉你(主观)先验信念的先验。然而，在实践中，先验分布的选择可能相当主观，有时是任意的。我们只需选择一个具有较大标准差(小精度)的正态分布先验；例如，我们可以假设β0和β1来自正态分布，平均值为0和SD 10，000。这被称为无信息的先验（uninformative prior），因为这个分布本质上是相当平坦的(也就是说，它将几乎相等的分布分配到特定范围内的任何值)。

接下来，如果这类分布是我们的选择，我们不需要担心哪种分布可能更好，因为它们在可能性不可忽略的区域几乎都是平坦的，而后验分布并不关心在可能性可以忽略的区域中先验分布是什么样子。

同样，对于精度τ，我们知道它们必须是非负的，所以选择一个限制于非负值的分布是有意义的，例如，我们可以使用一个形状和尺度参数较低的伽马分布（Gamma distribution）。

另一个有用的信息不足的选择是均匀分布。如果你为σ或τ选择均匀分布，你可能会得到一个由JohnK.Kruschke演示的模型：

用R和JAGS进行仿真

到目前为止对理论有好处。问题是，求解方程(2)在数学上是有挑战性的。在绝大多数情况下，后验分布是不可直接获得的(看看正态分布和Gamma分布有多复杂，你必须把这些分布相乘起来)。马尔可夫链蒙特卡罗方法（Markov Chain Monte Carlomethods）常被用来估计模型的参数。JAGS帮助我们做到了这一点。

JGS模型是一个基于马尔可夫链蒙特卡罗(MCMC)的模拟过程，它将产生参数空间θ=(β0；β1；τ)的多次迭代，该参数空间中每个参数生成的样本分布近似于参数的总体分布。

为什么是这样呢？这比这篇文章所能提供的解释要复杂得多。一个线性解释是：MCMC通过构造一个具有目标后验分布的马尔可夫链(？)来生成样本。

我告诉过你这并不有趣。快速的结论是：我们可以用数学证明我们的样本的分布是β0，β1，τ的真实分布，而不是解方程(2)，这通常在分析上不可能的，我们可以做一些聪明的抽样。

我们在R中运行JAGS，如下所示：

n <- nrow(data) #383 data points

mymodel <- "
model{
for(i in 1:n){
y[i] ~ dnorm(a + b*x[i], tau)
}
a ~ dnorm(0, 1e-6)
b ~ dnorm(0, 1e-6)
tau ~ dgamma(.01,.01)
sig <- 1/sqrt(tau)
}
"

然后，我们命令JAGS进行模拟。这里我要求JAGS模拟参数空间θ值10000次。

jm <- jags.model(file = textConnection(mymodel), 
data=list(n=n, x=data$height, y=data$weight))
cs <- coda.samples(jm, c("a","b","sig"), 11000)
sample_data <- as.data.frame(cs[[1]][-(1:1000),])

经过这样的采样，我们得到了θ=(β0；β1；τ)的采样数据：

现在我们有参数空间θ的10，000次迭代，请记住，通过方程(1)：

这意味着，如果我们用x=168 cm来代替每一次迭代，我们将得到10，000个权重值，因此，在高度为168cm的情况下，权重的分布：

cmean <- sample_data$a + sample_data$b*m_height 
# "conditional mean"

我们对体重的百分比感兴趣，这取决于我的身高。

m_perc <- pnorm(q = m_weight, mean = cmean, 
sd = sample_data$sig) 
truehist(m_perc, 
main = "Posterior distribution for 
my weight percentile", xlab = "percentile", 
ylab = "Frequency")