rmnlIndepMetrop            package:bayesm            R Documentation

_M_C_M_C _A_l_g_o_r_i_t_h_m _f_o_r _M_u_l_t_i_n_o_m_i_a_l _L_o_g_i_t _M_o_d_e_l

_D_e_s_c_r_i_p_t_i_o_n:

     'rmnIndepMetrop' implements Independence Metropolis for the MNL.

_U_s_a_g_e:

     rmnlIndepMetrop(Data, Prior, Mcmc)

_A_r_g_u_m_e_n_t_s:

    Data: list(p,y,X)

   Prior: list(A,betabar)  optional

    Mcmc: list(R,keep,nu) 

_D_e_t_a_i_l_s:

     Model:   y ~ MNL(X,beta). Pr(y=j) =
     exp(x_j'beta)/sum_k{e^{x_k'beta}}. 

     Prior:   beta ~ N(betabar,A^{-1}) 

     list arguments contain:

   '_p' number of alternatives

   '_y' nobs vector of multinomial outcomes (1,..., p)

   '_X' nobs*p x nvar matrix

   '_A' nvar x nvar pds prior prec matrix (def: .01I)

   '_b_e_t_a_b_a_r' nvar x 1 prior mean (def: 0)

   '_R' number of MCMC draws

   '_k_e_e_p' MCMC thinning parm: keep every keepth draw (def: 1)

   '_n_u' degrees of freedom parameter for independence t density (def:
        6) 

_V_a_l_u_e:

     a list containing:  

betadraw: R/keep x nvar array of beta draws

 loglike: R/keep vector of loglike values for each draw

 acceptr: acceptance rate of Metropolis draws

_A_u_t_h_o_r(_s):

     Peter Rossi, Graduate School of Business, University of Chicago,
     Peter.Rossi@ChicagoGsb.edu.

_R_e_f_e_r_e_n_c_e_s:

     For further discussion, see _Bayesian Statistics and Marketing_ by
     Rossi, Allenby and McCulloch, Chapter 5. 
      <URL:
     http://faculty.chicagogsb.edu/peter.rossi/research/bsm.html>

_S_e_e _A_l_s_o:

     'rhierMnlRwMixture'

_E_x_a_m_p_l_e_s:

     ##

     if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=2000} else {R=10}

     set.seed(66)
     n=200; p=3; beta=c(1,-1,1.5,.5)

     simmnl= function(p,n,beta) {
       #   note: create X array with 2 alt.spec vars
       k=length(beta)
       X1=matrix(runif(n*p,min=-1,max=1),ncol=p)
       X2=matrix(runif(n*p,min=-1,max=1),ncol=p)
       X=createX(p,na=2,nd=NULL,Xd=NULL,Xa=cbind(X1,X2),base=1)
       Xbeta=X%*%beta # now do probs
       p=nrow(Xbeta)/n
       Xbeta=matrix(Xbeta,byrow=TRUE,ncol=p)
       Prob=exp(Xbeta)
       iota=c(rep(1,p))
       denom=Prob%*%iota
       Prob=Prob/as.vector(denom)
       # draw y
       y=vector("double",n)
       ind=1:p
       for (i in 1:n) 
             { yvec=rmultinom(1,1,Prob[i,]); y[i]=ind%*%yvec }
        return(list(y=y,X=X,beta=beta,prob=Prob))
     }

     simout=simmnl(p,n,beta)

     Data1=list(y=simout$y,X=simout$X,p=p); Mcmc1=list(R=R,keep=1)
     out=rmnlIndepMetrop(Data=Data1,Mcmc=Mcmc1)

     cat("Summary of beta draws",fill=TRUE)
     summary(out$betadraw,tvalues=beta)

     if(0){
     ## plotting examples
     plot(out$betadraw)
     }

