Hi, all:
Can anybody check where is wrong with my code? I tried a lot of times, but
did not find an error. The parameters' estimator is not accurate. It's a
simple model about a multiple regression, with five covariates. rwmetrop is
supposed to give a much more accurate estimand. Thanks a lot.
rm(list=ls())
n=100; p=5;
xTrue=matrix(rnorm(n*p),nrow=n, ncol=p)
betaTrue=c(1,2,0,3,1)
yTrue=xTrue%*%betaTrue+rnorm(n)
d=list(y=yTrue, x=xTrue)
datapost=function(theta,data){
x=data$x
y=data$y
mu=rep(0,times=100)
for(j in 1:5){
mu=mu+x[,j]*theta[j]
}
logdensity=-(y-mu)^2/2-log(sqrt(2*pi))
sum(logdensity)
}
covariance=array(0,dim=c(p,p))
covariance[row(covariance)==col(covariance)]=1
proposal=list(var=covariance, scale=2)
start=c(1,1,1,1,1)
fit=rwmetrop(datapost, proposal, start, 100000, d)
colMeans(fit$par[50001:100000,])
[[alternative HTML version deleted]]