使用R2WinBUGS拟合后怎么绘制后验概率密度

相似文件 换一批

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

使用R2WinBUGS进行线性拟合后，
如何绘制验概率密度（Posterior probability density）
请问大神用哪个语句？
非常感谢

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

关键词：R2WinBUGS winbugs WINBUG 后验概率 概率密度 如何

没有在R中调用过，直接使用的WinBUGS，设定监测点，给出后验概率参数。

用coda包。

soccy 发表于 2014-12-28 09:41
用coda包。

谢谢，能给一个示例吗？看看怎么做的

http://www.jstatsoft.org/v12/i03/
http://people.orie.cornell.edu/davidr/or678/r2winbugs/index.html

https://sites.google.com/site/doingbayesiandataanalysis/software-installation

WinBUGS局限性很明显，如果刚开始学习bayesian分析，不如直接从Stan入手。

soccy 发表于 2014-12-31 00:35
WinBUGS局限性很明显，如果刚开始学习bayesian分析，不如直接从Stan入手。

大侠，我这个有点着急用，你看看这些代码之后怎么样才能绘出后验概率图？
请求您的帮忙，代码如下：

1. library(R2WinBUGS)
   
2. library(lattice)
   
3. library(coda)
   
4. dat<-read.csv("c:/data.csv")
   
5. attach(dat)
   
6. setwd("c:/")
   
7. sink("lme.model2.txt")
   
8. cat("
   
9. model {
   
10. 
    
11. # Priors
    
12. for (i in 1:6){               
    
13.     alpha[i] ~ dnorm(mu.int, tau.int)        # Random intercepts
    
14.     beta[i] ~ dnorm(mu.slope, tau.slope)# Random slopes
    
15. }
    
16. 
    
17. mu.int ~ dnorm(0, 0.001)                # Mean hyperparameter for random intercepts
    
18. tau.int <- 1 / (sigma.int * sigma.int)
    
19. sigma.int ~ dunif(0, 100)                # SD hyperparameter for random intercepts
    
20. 
    
21. mu.slope ~ dnorm(0, 0.001)                # Mean hyperparameter for random slopes
    
22. tau.slope <- 1 / (sigma.slope * sigma.slope)
    
23. sigma.slope ~ dunif(0, 100)                # SD hyperparameter for slopes
    
24. 
    
25. tau <- 1 / ( sigma * sigma)                # Residual precision
    
26. sigma ~ dunif(0, 100)                        # Residual standard deviation
    
27. 
    
28. # Likelihood
    
29. for (i in 1:n) {
    
30.     mass[i] ~ dnorm(mu[i], tau)
    
31.     mu[i] <- alpha[pop[i]] + beta[pop[i]]* length[i]
    
32. }
    
33. }
    
34. ",fill=TRUE)
    
35. sink()
    
36. 
    
37. # Bundle data
    
38. win.data <- list(mass = as.numeric(mass), pop = pop, length = length, n =length(mass))
    
39. 
    
40. # Inits function
    
41. inits <- function(){ list(alpha = rnorm(6, 0, 2), beta = rnorm(6, 10, 2),
    
42. mu.int = rnorm(1, 0, 1), sigma.int = rlnorm(1), mu.slope = rnorm(1, 0, 1),
    
43. sigma.slope = rlnorm(1), sigma = rlnorm(1))}
    
44. 
    
45. # Parameters to estimate
    
46. parameters <- c("alpha", "beta", "mu.int", "sigma.int", "mu.slope", "sigma.slope", "sigma")
    
47. 
    
48. ni <-10000
    
49. nb <- 1000
    
50. nt <- 2
    
51. nc <- 3
    
52. 
    
53. # Start Gibbs sampling
    
54. out1 <- bugs(win.data, inits, parameters, "lme.model2.txt", n.thin=nt,
    
55. n.chains=nc, n.burnin=nb, n.iter=ni, debug = TRUE)

复制代码

plot(out1)如何？

soccy 发表于 2014-12-31 19:17
plot(out1)如何？

不对的，我试过的，不知还有其他方法吗？数据也上传了，代码在上面，麻烦您试试怎么绘出？谢谢您了？