[R] Discrete choice model maximum likelihood estimation

infinitehorizon barisvardar at hotmail.com
Wed May 16 02:26:55 CEST 2012


Hello John,

Thanks for your response. Can you explain why it is better to use NM or nmk
instead of optim? By the way, my original model consists 24 parameters. In
this case is it still better to use NM or nmk?

Best,

Marc



John C Nash wrote
> 
> Your function will not evaluate as coded. i.e., llfn(start.par) doesn't
> "work", as there
> are unequal numbers of arguments. Also, while R allows you to use
> variables that are not
> explicitly defined for a function, I've almost always got into trouble if
> I don't pass
> them VERY carefully.
> 
> Finally, as author of the code used to put CG in optim, I'll advise
> against its use. One
> of my least successful pieces of code. Rcgmin is better, but you really do
> need analytic
> derivatives to make it sing. For 5 parameters, use NM, or better the nmk
> from dfoptim package.
> 
> Best, JN
> 
> 
> On 05/15/2012 06:00 AM, r-help-request@ wrote:
>> Message: 13
>> Date: Mon, 14 May 2012 04:21:57 -0700 (PDT)
>> From: infinitehorizon <barisvardar@>
>> To: r-help@
>> Subject: Re: [R] Discrete choice model maximum likelihood estimation
>> Message-ID: <1336994517063-4629910.post at .nabble>
>> Content-Type: text/plain; charset=us-ascii
>> 
>> Hello again,
>> 
>> I changed the name to tt. 
>> and for a and tt actually I was getting them from data, I didn't put them
>> here in the question. Now I restructured my code and below I copy the
>> full
>> code, I tried many things but still getting the same error, I don't
>> understand where is the mistake.
>> 
>> I also defined one more variable to increase comprehension of the
>> problem.
>> Instead of data, I define three representative vectors in the beginning.
>> If
>> you run this code, you will see the error message.
>> 
>> # Variables, a: discrete choice variable-dependent, x and tt are
>> independent
>> variables, tt is binary
>> a  <- c(1,1,2,3,2,3,1,2,2,3,1,1)
>> x  <- c(23,26,12,27,10,30,13,20,23,44,17,15)
>> tt <- c(1,0,0,1,1,1,0,0,1,0,1,1)
>> 
>> # First, just to see, the linear model
>> 
>> LM	<-lm(a ~ x+tt)
>> coefLM	<- coefficients(LM)
>> summary(LM)
>> 
>> # Probabilities for discrete choices for a=3, a=2 and a=1 respectively 
>> P3 <- function(bx,b3,b,tt) { 
>> P <- exp(bx*x+b3+b*(tt==1))/(1-exp(bx*x+b3+b*(tt==1))) 
>> return(P) 
>> } 
>> P2 <- function(bx,b2,b,tt) { 
>> P <- exp(bx*x+b2+b*(tt==1))/(1-exp(bx*x+b2+b*(tt==1))) 
>> return(P) 
>> } 
>> P1 <- function(bx,b1,b,tt) { 
>> P <- exp(bx*x+b1+b*(tt==1))/(1-exp(bx*x+b1+b*(tt==1))) 
>> return(P) 
>> }
>> 
>> # Likelihood functions for discrete choices for a=3, a=2 and a=1
>> respectively
>> 
>> L3 <- function(bx,b1,b2,b3,b,tt) { 
>> P11 <- P1(bx,b1,b,tt) 
>> P22 <- P2(bx,b2,b,tt) 
>> P33 <- P3(bx,b3,b,tt) 
>> 
>> L3l <- P11*P22*P33 
>> return(L3l) 
>> } 
>> 
>> L2 <- function(bx,b1,b2,b,tt) { 
>> P11 <- P1(bx,b1,b,tt) 
>> P22 <- P2(bx,b2,b,tt) 
>> 
>> L2l <- P11*P22 
>> return(L2l) 
>> } 
>> 
>> L1 <- function(bx,b1,b,tt) { 
>> P11 <- P1(bx,b1,b,tt) 
>> 
>> L1l <- P11 
>> return(L1l) 
>> }
>> 
>> # Log-likelihood function 
>> 
>> llfn <- function(param) { 
>> 
>> bx <- param[1] 
>> b1 <- param[2] 
>> b2 <- param[3] 
>> b3 <- param[4] 
>> b <- param[5] 
>> 
>> lL1 <- log(L1(bx,b1,b2,b,tt)) 
>> lL2 <- log(L2(bx,b1,b2,b3,b,tt)) 
>> lL3 <- log(L3(bx,b1,b2,b3,b,tt)) 
>> 
>> llfn <- (a==1)*lL1+(a==2)*lL2+(a==3)*lL3 
>> } 
>> start.par <- c(0,0,0,0,0) 
>> est <- optim(param=start.par,llfn,x=x, a=a, tt=tt, method =
>> c("CG"),control=list(trace=2,maxit=2000), hessian=TRUE)
>>
> 
> ______________________________________________
> R-help@ mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 


--
View this message in context: http://r.789695.n4.nabble.com/Discrete-choice-model-maximum-likelihood-estimation-tp4629877p4630199.html
Sent from the R help mailing list archive at Nabble.com.



More information about the R-help mailing list