发帖
雷达卡
如题，利用模拟数据构建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)
京公网安备 11010802022788号
