2012-05-14 75 views
3

我看着an excellent post on Bayesian Linear Regression (MHadaptive)置信區間

給人一種輸出 後置信區間

BCI(mcmc_r) 
#            0.025    0.975 
# slope      -5.3345970 6.841016 
# intercept    0.4216079 1.690075 
# epsilon   3.8863393 6.660037 

什麼功能,我現在用它來構建與置信區間的模型從這些參數?

+4

我不understand-後置信區間* *是置信區間的貝葉斯等同。爲什麼不是你想要的? –

+0

您是否在尋求截距和斜率a和b的點估計?如果是這樣,沒有自己運行代碼,我懷疑你只是使用mean(a)和mean(b)。 –

+0

您是否下載並查看發佈在本頁底部的文件?我懷疑在那裏提供的R代碼會返回點估計值和可信區間。 https://sites.google.com/site/mcgillbgsa/workshops/bayesian –

回答

8

爲什麼不使用從MCMC獲得的分佈來預測y從任何點x的分佈?在您使用的例子,這裏有有關章節,其中eggmass = y和長度= X

##@ 3.1 @## 

## Function to compute predictions from the posterior 
## distribution of the salmon regression model 
predict_eggmass<-function(pars,length) 
{ 
    a <- pars[, 1]  #intercept 
    b <- pars[, 2]  #slope 
    sigma <- pars[, 3] #error  
    pred_mass <- a + b * length 
    pred_mass <- rnorm(length(a), pred_mass, sigma) 
    return(pred_mass) 
} 

### -- ### 
##@ 3.2 @## 

## generate prediction 
pred_length <- 80  # predict for an 80cm individual 
pred <- predict_eggmass(mcmc_salmon$trace, length=pred_length) 
## Plot prediction distribution 
hist(pred, breaks=30, main='', probability=TRUE) 

## What is the 95% BCI of the prediction? 
pred_BCI <- quantile(pred, p=c(0.025, 0.975)) 
    2.5% 97.5% 
33.61029 43.16795 

我想你指的是在您的評論的分佈可這裏pred和置信區間爲pred_BCI

+1

我認爲OP可能需要置信區間(即在'predict_eggmass'函數中跳過'rnorm'步驟),而不是預測區間(這裏就是這個)。 –

1

如果您想要查看每個參數的後驗邊際密度,可以使用density()以及存儲在mcmc_r對象的trace組件中的樣本。

library(MHadaptive) 
data(mcmc_r) 

BCI(mcmc_r) 
#    0.025 0.975 post_mean 
# a  -6.6113522 7.038858 0.001852978 
# b  0.2217377 1.543519 0.902057671 
# epsilon 3.8094802 6.550360 4.957292114 

head(mcmc_r$trace) 
#   [,1]  [,2]  [,3] 
# [1,] 3.1448136 0.7211228 5.449728 
# [2,] 2.2287645 0.7155189 4.602004 
# [3,] 2.0812509 0.8035820 4.224071 
# [4,] 1.2444855 0.8448825 4.737466 
# [5,] 3.2765630 0.5947548 4.740052 
# [6,] 0.4271876 0.9014841 5.333821 

plot(density(mcmc_r$trace[,3]), main=mcmc_r$par_names[3]) 

enter image description here