全部版块 我的主页
论坛 计量经济学与统计论坛 五区 计量经济学与统计软件 winbugs及其他软件专版
2792 1
2013-12-17
如题,利用模拟数据构建probit回归模型,但总是出错。在winbug论坛看见如下描述,大意为probit模型超边,请问有大神遇到过同样问题么?如何解决?
a) 'undefined real result' indicates numerical overflow. Possible reasons include:
- initial values generated from a 'vague' prior distribution may be numerically extreme - specify appropriate initial values;
- numerically impossible values such as log of a non-positive number - check, for example, that no zero expectations have been given when Poisson modelling;
- numerical difficulties in sampling. Possible solutions include:
- better initial values;
- more informative priors - uniform priors might still be used but with their range restricted to plausible values;
- better parameterisation to improve orthogonality;
- standardisation of covariates to have mean 0 and standard deviation 1.
- can happen if all initial values are equal.

Probit models are particularly susceptible to this problem, i.e. generating undefined real results. If a probit is a stochastic node, it may help to put reasonable bounds on its distribution, e.g.

   probit(p) <- delta
   delta ~ dnorm(mu, tau)I(-5, 5)

This trap can sometimes be escaped from by simply clicking on the update button. The equivalent construction

   p <- phi(delta)

may be more forgiving.



程序源代码如下:library(MASS)
library(R2WinBUGS)
n.site <- 1500
g1<- rbinom(n = n.site, size=1, prob=0.4)   
c1<- rnorm(n = n.site, mean=0, sd=1)
c2<- rbinom(n = n.site, size=3, prob=0.3)
u<- rnorm(n = n.site, mean=1, sd=2)
xb <- 0+ 0.4*g1+0.6*c1-0.8*c2+u
probx <- 1/(1+exp(-xb))
x <- rbinom(n = n.site, size = 1, prob = probx)

sink("model.txt")                    #模型1-简单两阶段probit回归法
cat("
    model {
    # Priors
    alpha0 ~ dnorm(0,0.01)
    alpha1 ~ dnorm(0,0.01)
    alpha2 ~ dnorm(0,0.01)
    alpha3 ~ dnorm(0,0.01)
    alpha4 ~ dnorm(0,0.01)
    # Likelihood
    for (i in 1:n) {

    x ~ dbern(px)        # Note p before N
    probit(px) <- alpha0+alpha1*g1+alpha2*c1+alpha3*c2+alpha4*u

    }
    }
    ",fill=TRUE)
sink()

# Bundle data
win.data <- list(g1=g1,c2=c2,u=u,c1=c1,n = length(g1), x=x)

# Inits function
inits <- function(){ list(alpha0=0,alpha1=0,alpha2=0,alpha3=0,alpha4=0) }

# Parameters to estimate
params <- c("alpha0", "alpha2", "alpha1", "alpha3", "alpha4")

# MCMC settings
nc <- 3 #Number of Chains
ni <- 1200 #Number of draws from posterior
nb <- 200 #Number of draws to discard as burn-in
nt <- 2 #Thinning rate

# Start Gibbs sampling
out <- bugs(data=win.data, inits=inits, parameters.to.save=params,
model.file="model.txt", n.thin=nt, n.chains=nc, n.burnin=nb,
n.iter=ni, debug = TRUE,bugs.directory="D:/Program Files/WinBUGS14/")
print(out,dig=3)




二维码

扫码加我 拉你入群

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

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

全部回复
2013-12-30 20:19:03
This trap can sometimes be escaped from by simply clicking on the update button
二维码

扫码加我 拉你入群

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

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

相关推荐
栏目导航
热门文章
推荐文章

说点什么

分享

扫码加好友,拉您进群
各岗位、行业、专业交流群