发帖
雷达卡
fit.xls
(15.5 KB)
我在利用R2WinBUGS计算时，编译通不过，以下是我的代码，附件是数据。请高手帮我看看，特别是数据导入那块，先谢谢了！
library(R2WinBUGS)
library(lattice)
library(coda)
#write model
sink("curve.txt")
cat("
model
{
for (i in 1:N)
{
pred<-1.3+a*pow(d,b)
h~dnorm(pred,prec)
}
# Priors
a ~ dnorm(2.0,1.0E-6)
b ~ dnorm(0.3,1.0E-6)
prec~dgamma(0.001,0.001) # prior for the precision
# Assess model fit using a sums-of-squares-type discrepancy
for (i in 1:N)
{
residual <- h-pred # Residuals
predicted <- pred # Predicted values
sq <- pow(residual, 2) # Squared residuals for observed data
# Generate replicate data and compute fit stats for them
h.new ~ dnorm(pred, prec) # one new data set at each MCMC iteration
sq.new <- pow(h.new-predicted, 2) # Squared residuals for new data
}
fit <- sum(sq[]) # Sum of squared residuals for actual data set
fit.new <- sum(sq.new[]) # Sum of squared residuals for new data set
bpvalue <- mean(test) # Bayesian p-value
}
",fill=TRUE)
sink()
hd <- read.csv("fit.csv", header=TRUE)
# Inits function
inits <- function(){ list(a=rnorm(1.5), b=rnorm(0.1), prec=rlnorm(1))}
# Parameters to estimate
params <- c("a","b", "prec", "p.decline", "fit", "fit.new", "bpvalue","residual", "predicted")
# MCMC settings
nc=1 ; ni=10000 ; nb=200 ; nt=1
# Start Gibbs sampler
out <- bugs(hd, inits=inits, parameters=params, model="curve.txt", n.thin=nt, n.chains=nc, n.burnin=nb, n.iter=ni, debug=TRUE)
print(out, dig=3)
京公网安备 11010802022788号
