Hi,
I tried to use rmnlIndepMetrop (bayesm package) for my MNL model with 4
choice alternatives, 5 independent variables, 69 observations,
dim(X) [1] 276 5, nu=6. So I run such code:
if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=2000} else {R=10}
set.seed(66)
df=read.table("X_metrop.dat",header=TRUE)
inp=as.matrix(df)
y=as.numeric(inp[,1])
n=length(y)
p=4
nvar=ncol(X)
X=cbind(1,inp[,2:5])
XTXI=chol2inv(chol(crossprod(X)))
bhat=XTXI%*%t(X)%*%y
#create X matrix with dim nobs*p x nvar
index=1:n
index2=rep(index,p)
index3=sort(index2)
X=X[index3,]
beta=array(bhat,dim=c(nvar))
A=diag(c(rep(.01,length(beta))));
betabar=rep(0,length(beta))
Data=list(y=y,X=X,p=p);
Mcmc=list(R=R,keep=1) ;
Prior=list(A=A,betabar=betabar)
out=rmnlIndepMetrop(Data=Data,Prior=Prior,Mcmc=Mcmc)
Ouput:
Starting Independence Metropolis Sampler for Multinomial Logit Model
with 4 alternatives
Prior Parms:
betabar
[1] 0 0 0 0 0
A
[,1] [,2] [,3] [,4] [,5]
[1,] 0.01 0.00 0.00 0.00 0.00
[2,] 0.00 0.01 0.00 0.00 0.00
[3,] 0.00 0.00 0.01 0.00 0.00
[4,] 0.00 0.00 0.00 0.01 0.00
[5,] 0.00 0.00 0.00 0.00 0.01
MCMC parms:
R= 10 keep= 1 nu (df for st candidates) = 6
Mistake in chol(mhess) : the leading Minor of the order 1 is not positive
defined
By debugging: debug(rmnlIndepMetrop)
rmnlIndepMetrop(Data=Data,Prior=Prior,Mcmc=Mcmc)
beta
[1] 1.46867390 -0.52889086 0.23710427 0.00188178 0.06021696
mhess=mnlHess(y,X,beta)
> mhess
X0 X0.1 X0.2 X0.3
0 0 0 0 0
X0 0 0 0 0 0
X0.1 0 0 0 0 0
X0.2 0 0 0 0 0
X0.3 0 0 0 0 0
>chol(mhess)
Mistake in chol(mhess) : the leading Minor of the order 1 is not positive
defined
Is something wrong with my data? But I tried another dataset and got the
same mistake.
Whether everybody knows how I could manage it?
I have no ideas anymore(.
Thank you in advance,
Tatjana
-----------------------------------
Tetyana Stepanchuk
Department of Electronic Commerce
Johann Wolfgang Goethe-University
Mertonstrasse 17
D-60054, Frankfurt/Main
Germany
phone: +49 069 798 22380
http://www.ecommerce.wiwi.uni-frankfurt.de
[[alternative HTML version deleted]]
______________________________________________
[email protected] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
