• 综合主站
  • WinBUGS在统计分析中的应用(第二部分)

谢谢楼主的讲座和书!
to priss111:谢老大和楼主都很牛--这么说就对啦!
4 个月 后
很感谢楼主的文章,让我受益很多!
我是winbugs的新手,我在做时遇到了一个问题,向您请教一下。模型是这样的: y=e+J*P, 期中 e ~ dnorm(u0,tau0), J~dnorm(u1,tau1), P~dpois(lanta), y的数据是已知的,我想要求出参数u0,u1,tau0,tau1,and lanta. 假设先验分布: u0,u1~dnorm, tau0,tau1~dgamma, lanta~ dbeta, and y~ dnorm,这些分布是有文献支持的. 但是我发现这样做出来的结果不好:首先我先假设参数的值,产生y的数据集,然后用y来求出这些参数,对比一下,发现差距比较大。我想是否我代码编错了?
model
{
for(i in 1:N) {
yy ~ dnorm (y,tau)
cq ~ dnorm (u1,tau1)
bq ~ dnorm(u0,tau0)
dq ~ dpois(lanta)
y <- bq+cq*dq
}
u1 ~dnorm (0.01, 0.0002)
u0 ~dnorm (1, 0.0002)
tau0~dgamma (4, 2)
tau1 ~dgamma (4, 2)
lanta~dbeta (2,1)
tau<- 1/(1/tau0+1/tau1*pow(1-exp(-lanta),2))
}
2 个月 后
[未知用户] 问题不是出在你的model上,而是出在你的实验设计上。你先假设参数的值,产生y的数据集,然后用y来求出这些参数,理论上来说,只有y无限多,且产生这些y的参数遍历其所有空间时,才能返回来求出这些参数的真值。举一个简单的例子:在normal distribution P(y|theta, tau)中,我已知theta=0,tau=1,用这样的参数来产生y;假如我产生了5次data, 非常非常巧,我5次data的值都是0,因此你的数据集y={0, 0, 0, 0, 0}.那么你能用这样的数据集来估计出真正的theta和tau么?答案是否定的,你只能估计出正确的theta,而tau却非常的不一样。所以说,你产生的data要能正确地反映出你model,你才能正确的估计它的参数。这是一个parameter identifiability 的问题。希望这样的解答对你有用。
[未知用户] icwei您好,也很想请教您一个问题:

您说的"问题不是出在你的model上,而是出在你的实验设计上”
这句话中的‘出在实验设计上’怎么理解,能不能再通俗的解释一下。

从您这句话中也受到启发。 因为我也遇到了一个这样的问题,
当然,不排除还有其他问题的可能性。
[未知用户] 我的意思是原作者有一个model,他为了检验这个model的正确性,便设计了一个实验,这个实验就是先假设所有参数已知,用这些参数产生data,然后再用这些data来估计这些参数,对比求得的参数和原参数的值来证明model是否正确。这样的实验没有大问题,但实施的时候要小心,你需要的检验你所产生的这些data是否能反映model的所有特性。

实际上,原作者提出的验证model的方法并不正确,即使他能够得到和原参数一样的值也并不能说明这个model对于实际的data就是正确的。他犯了一个基本的错误,既他用frequentist的方法来验证bayesian的正确性。在frenquentist的概念中,model是存在真值的,所以你可以提出一个hypothesis去检验估计值是否等于真值;而在bayesian中,我们是用credit interval 来表示所求的参数,既真值并不存在,而是与我们的prior knowledge相关的一个范围。

实际问题中,我们更多的是只得到一些数据,然后用这些数据去估计我们的参数。正确的判断参数是否能够反映系统正确性的方法是看posterior predictive distribution, 既P(y_pred|y). 如果你的posterior predictive distribution能够正确的模拟你的data,那说明这个model所估计的参数是真正能够反映你系统的特性的;反之不能。
[未知用户] 不好意思,我学统计学的时候用英文学的,所以可能用中文表达起来有些难懂。举个不完全恰当但简单的例子吧,像我上面说的,我有dataset={0, 0, 0, 0, 0},他们是从一个有噪音的系统产生的,噪音符合normal distribution,但均值和方差均不为0。那么我用一个normal distribution的model去fit这些data,得到的结果为mean=0, variance=0.这个model可以完美的fit我当前的data.可是当我用这个model求posterior predictive distribution的时候,我所得到的是P(y_pred|y)=N(0, 0), 其所相对应得predictive dataset为y_pred={0, 0, 0, 0, 0, .....}. 可是当我们让这个系统继续产生更多的data时,我们会发现我们model产生的值并不能很好的反映系统的特性(尽管估计的参数能够很好的fit我们原始的dataset)。

那么,我们的疑问是:我们的model错了么???但在这个特殊的例子里,我们可以说,model没有错,只是他估计的参数错了。那么参数不正确是由什么引起的呢?那就是我们的实验设计。这个实验设计为只得到5个数据,而这5个数据不足以反映系统的全部特性,因此才估计出错误的参数。

因此,我们说,参数是否能代表系统的特性,不但model要能反映真实系统的特性,数据也要能够反映真实系统的特性,这涉及到实验设计的问题。通常情况下,只要考虑posterior predictive distribution是否能反映系统产生的未来数据的特性就可以了。
谢谢icwen这么认真的回复,从中又受到一些启发。 哦,我基本上明白: 对于bayesian中的参数估计,需要用后验分布来检验所估计的参数是否能真实的反应实际情形;而其他类型参数估计的模拟实验,则直接用所估计到的参数与模拟时设定的参数进行比较,越吻合说明该估计方法对所估计的参数即有效也稳健。

不知道这样理解对不对。 另外,还想请教您一个问题:不只您是不是学数学的背景?
[未知用户] 你理解的基本正确,但注意不要混淆posterior distribution(后验分布)和posterior predictive distribution(后验估计分布?不知道这个中文名是否正确)的概念。他们的区别是:假设我们的参数为a,数据为y,那么后验分布是指P(a|y),而后验估计分布是指P(y_pred|y),它等于P(y_pred|y)=integral(P(y_pred|a)*P(a|y))da,既你对未来数据的预测是基于当前model的likelihood,以后验分布作为piror,进行加权而得到的。他是对你估计参数的所有值得一个综合性考量。

我不是纯数学背景出身,我大本学的是自控和信号处理,从博士开始专攻统计信号处理,博士毕业后开始转入统计学在生物医学方面的数学建模,现在主要兴趣在Bayesian statistics, mcmc, hierarchical graphical model, differential equation 以及 stochastic system。有兴趣的话可以和我联系,我们可以共同探讨。我的email是 chen.wei@mrc-bsu.cam.ac.uk
2 年 后
感谢楼主的文章和各位同学的讨论,在下受益匪浅。
在我最近的一篇文章中,我打算小试一下贝叶斯模型,在经过很多次尝试后,都已失败告终,感觉很深奥,所以想和各位探讨一二。
我的问题是找出各林分因子对树木死亡率的影响,即stand basal area(SBA),the ratio of focal species basal area to stand basal area (FSBA), stand age (SA)对 tree mortality rate 影响。假设我总共有PN个sample plots,每个plots,我有N个树是起初的活的株数,有S棵树是在一段时间(time)后的活树数量(不包括后来长大的),一般情况,每个plot的树木死亡率通过(ln(N)-ln(S))/time来计算获得,可是我的有些plot(一般以上的plots),树木没有死亡,或者是全部死亡,即N=S,or S=0,所以我想参照Condit R.(2006)和 Kraft J. N.(2010)的hieracichal bayies算死亡率(mp)的方法,把以上两种情况的plots都用在实际分析中。
他们的idea是这样的,对th species (plots, in my cases),树木的存活率(sp),N,S满足binomial distribution,即 S~dbin(sp,N),and, sp=exp(-mp*time),mp为死亡率,这样就实现了mp=(ln(N)-ln(S))/time这个运算过程,然后他们又假设among species,死亡率是满足lognormal distribution,并且得到数据上的验证,即有strong right skews,所以 mp~dlnorm(logmu,logsd),那么,我现在想和大家探讨的第一个问题是,我这样写出他们的winbugs model是不是正确的?
model
{ for(i in 1:speciesnumber)
{ S[i]~dbin(sp[i],N[i])
sp[i]<-exp(-mp[i]*time[i])
mp[i]~dlnorm(logmu,tau)
}
logmu~dnorm(1,0.00000001)
tau~dgamma(1,1)
}
接下来,就到我的case中了,我检查的一下我的among plots mortality rate distribution,我发现也stong right skews,所以我也想用了上面的方法得mp,然后我通过mp=a+b*SBA+c*FSBA+d*SA,得出各系数a,b,c,d。我的model是这样写的:
model
{ for (i in 1:PN)
{ S[i]~dbin(sp[i],N[i])
sp[i]<-exp(-mp[i]*time)
mp[i]~dlnorm(logmu,tau)
mp[i]<-a+b*SBA+c*FSBA+d*SA
}
logmu~dnorm(1,0.0000001)
tau~dgamma(1,1)
a~dnorm(1,0.0000001)
b~dnorm(1.0.0000001)
c~dnorm(1.0.0000001)
d~dnorm(1.0.0000001)
}
可是我这样写,却被提示“multiple definisitons of node mp[1]”,我想请问一下,我的model到底出了什么错误?
一万个感谢!

Condit, R., Ashton, P., Bunyavejchewin, S., Dattaraja, H.S., Davies, S., Esufali, S., Ewango, C., Foster, R., Gunatilleke, I.A.U.N., Gunatilleke, C.V.S., Hall, P., Harms, K.E., Hart, T., Hernandez, C., Hubbell, S., Itoh, A., Kiratiprayoon, S., LaFrankie, J., de Lao, S.L., Makana, J.R., Noor, M.N.S., Kassim, A.R., Russo, S., Sukumar, R., Samper, C., Suresh, H.S., Tan, S., Thomas, S., Valencia, R., Vallejo, M., Villa, G., and Zillio, T. 2006. The importance of demographic niches to tree diversity. Science 313(5783): 98-101.
Kraft, N.J.B., Metz, M.R., Condit, R.S., and Chave, J. 2010. The relationship between wood density and mortality in a global tropical forest data set. New Phytologist 188(4): 1124-1136. doi: 10.1111/j.1469-8137.2010.03444.x
8 个月 后
我想要用MCMC方法估计参数,看到您这篇文章以后做了个程序,可是原理不是很清楚,不知道对不对,cigam总是不出结果。
分布函数是F(x)=(1-F(x) )G(x-k)+F(k)=1-〖[1+ε((x-k)/σ)]〗^(-1⁄ε)
ε服从Pareto(a,c)

σ服从gamma(a,b)
不知道这样写程序是否正确呢?

model
{
for(i in 1:N)
{
p <-tao*pow(1+cigam*tao*(x-k),-1/cigam-1)
}
tao~ dgamma(0.01,0.01)
cigam~ dpar(0.1,0.1)
}
data
list(N=53,k=308,x=c(309,310,311,313,313,316,319,320,328,330,332,332,338,341,341,349,351,354,356,365,368,370,371,375,384,384,385,389,389,391,396,408,417,418,418,435,447,448,460,463,465,473,500,507,521,522,554,618,702,754,842,858,1152))
initial values
list(tao=1,cigam=1)

谢谢,麻烦你啦嘻嘻
1 年 后
[未知用户] icwei,您好!有个问题请教您:我有个问题想请教您:对于无信息先验分布(non - informative prior)是怎么定义的?我看了一些材料,多处都是按照以下方式定义的。但是对于无信息先验分布,方差不是应该足够大吗(In case of non-informative priors, the variances should have been specified as very large or instead uniform distributions covering the entire range of plausible values should have been used. )?能否帮我解释下?Thanks.
model
{
for (i in 1:N)
{
pd<-pow(d,c)
pred<-1.3+a*(1-exp(-b*pd))

h~dnorm(pred,prec)
}

# Priors
a ~ dnorm(0, 1.0E-6)
b ~ dnorm(0, 1.0E-6)
c ~ dnorm(0, 1.0E-6)
prec~dgamma(0.001,0.001)

}
4 个月 后
带入您的代码,发现如下提示:
this chain contains uninitialized variables
请更正...谢谢!