2017-06-04 66 views
0

我在网上发现这个问题。有人可以详细解释一下,为什么使用OLS更好?仅仅是因为样本数量不够?另外,为什么不使用全部1000个样本来估计先前的分布?贝叶斯vs OLS

我们有1000个随机抽样的数据点。目标是尝试构建一个回归模型,其中包含来自k个回归因子 变量的一个响应变量。哪个更好? 1.(贝叶斯回归)使用第一个 500个样本来估计假设的先验分布的参数,然后使用最后的500个样本更新 之前的后验分布,后验分布用于 最终回归模型。 2.(OLS回归)使用一个简单的普通 最小二乘回归模型与所有1000个回归变量

回答

3

“更好”永远是见仁见智,这在很大程度上取决于上下文。

对于频率主义OLS方法的优点:更简单,更快,更容易被更广泛的受众使用(因此解释更少)。我的一位明智的教授曾经说过:“当苍蝇拍将会做到这一点时,你不需要制造原子粉碎机。”

优势为等效的贝叶斯方法:更多的灵活性,以进一步模型开发,可直接模型导出/计算量的后验(有更多,但这些一直是我的动机去贝叶斯给定分析)。请注意单词“等效” - 您可以在贝叶斯框架中执行某些您无法在频率方法中执行的操作。

嘿,这是R的一次探索,首先模拟数据,然后使用典型的OLS方法。

N <- 1000 
x <- 1:N 
epsilon <- rnorm(N, 0, 1) 
y <- x + epsilon 

summary(lm(y ~ x)) 
## 
## Call: 
## lm(formula = y ~ x) 
## 
## Residuals: 
##  Min  1Q Median  3Q  Max 
## -2.9053 -0.6723 0.0116 0.6937 3.7880 
## 
## Coefficients: 
##    Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 0.0573955 0.0641910 0.894 0.371  
## x   0.9999997 0.0001111 9000.996 <2e-16 *** 
## --- 
## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
## 
## Residual standard error: 1.014 on 998 degrees of freedom 
## Multiple R-squared:  1, Adjusted R-squared:  1 
## F-statistic: 8.102e+07 on 1 and 998 DF, p-value: < 2.2e-16 

...这里是一个等效的贝叶斯回归,在回归参数和所有1000个数据点上使用非信息性先验。

library(R2jags) 
cat('model { 
    for (i in 1:N){ 
    y[i] ~ dnorm(y.hat[i], tau) 
    y.hat[i] <- a + b * x[i] 
    } 
    a ~ dnorm(0, .0001) 
    b ~ dnorm(0, .0001) 
    tau <- pow(sigma, -2) 
    sigma ~ dunif(0, 100) 
}', file="test.jags") 

test.data <- list(x=x,y=y,N=1000) 
test.jags.out <- jags(model.file="test.jags", data=test.data, 
        parameters.to.save=c("a","b","tau","sigma"), n.chains=3, n.iter=10000) 
test.jags.out$BUGSoutput$mean$a 
## [1] 0.05842661 
test.jags.out$BUGSoutput$sd$a 
## [1] 0.06606705 
test.jags.out$BUGSoutput$mean$b 
## [1] 0.9999976 
test.jags.out$BUGSoutput$sd$b 
## [1] 0.0001122533 

请注意,参数估计值和标准误差/标准偏差基本相等!

现在这里是另一个贝叶斯回归,它使用前500个数据点来估计先验,然后是最后500个估计后验。

test.data <- list(x=x[1:500],y=y[1:500],N=500) 
test.jags.out <- jags(model.file="test.jags", data=test.data, 
        parameters.to.save=c("a","b","tau","sigma"), n.chains=3, n.iter=10000) 

cat('model { 
    for (i in 1:N){ 
    y[i] ~ dnorm(y.hat[i], tau) 
    y.hat[i] <- a + b * x[i] 
    } 
    a ~ dnorm(a_mn, a_prec) 
    b ~ dnorm(b_mn, b_prec) 
    a_prec <- pow(a_sd, -2) 
    b_prec <- pow(b_sd, -2) 
    tau <- pow(sigma, -2) 
    sigma ~ dunif(0, 100) 
}', file="test.jags1") 

test.data1 <- list(x=x[501:1000],y=y[501:1000],N=500, 
        a_mn=test.jags.out$BUGSoutput$mean$a,a_sd=test.jags.out$BUGSoutput$sd$a, 
        b_mn=test.jags.out$BUGSoutput$mean$b,b_sd=test.jags.out$BUGSoutput$sd$b) 
test.jags.out1 <- jags(model.file="test.jags1", data=test.data1, 
        parameters.to.save=c("a","b","tau","sigma"), n.chains=3, n.iter=10000) 

test.jags.out1$BUGSoutput$mean$a 
## [1] 0.01491162 
test.jags.out1$BUGSoutput$sd$a 
## [1] 0.08513474 
test.jags.out1$BUGSoutput$mean$b 
## [1] 1.000054 
test.jags.out1$BUGSoutput$sd$b 
## [1] 0.0001201778 

有趣的是,推论与OLS结果相似,但差不多如此。这导致我怀疑用于训练先前的500个数据点在分析中没有像过去的500个那样重要,并且先前实际上已经被淘汰,尽管我不确定这一点。无论如何,我想不出所有1000个数据点(以及无法提供信息的先验)的原因,特别是因为我怀疑500 + 500会使用前500和后500的不同方法。

或许,所有这一切的答案是:我相信OLS和1000点贝叶斯结果超过500 + 500,而OLS更简单。

+0

非常感谢马特! – Erin