Deepankar, Some general advice from a non-expert:
Write your likelihoods without a for loop. This is important because the likelihood is evaluated multiple times in the maximization process and you don't want to be looping looping looping ... Always try multiple starting values Sometimes it helps to try different optimization functions (e.g., optim) Make sure your likelihood is correct. Check it against existing software if possible If necessary, simplify your model, to a single parameter even, and build it up from there Generate data under the model and see if your estimates are getting close to the truth Good luck, Stephen Rochester, Minn. USA On 4/18/07, Deepankar Basu <[EMAIL PROTECTED]> wrote: > As part of carrying out a complicated maximum likelihood estimation, I > am trying to learn to program likelihoods in R. I started with a simple > probit model but am unable to get the code to work. Any help or > suggestions are most welcome. I give my code below: > > ************************************ > mlogl <- function(mu, y, X) { > n <- nrow(X) > zeta <- X%*%mu > llik <- 0 > for (i in 1:n) { > if (y[i]==1) > llik <- llik + log(pnorm(zeta[i,], mean=0, sd=1)) > else > llik <- llik + log(1-pnorm(zeta[i,], mean=0, sd=1)) > } > return(-llik) > } > > women <- read.table("~/R/Examples/Women13.txt", header=TRUE) # DATA > > # THE DATA SET CAN BE ACCESSED HERE > # women <- > read.table("http://wps.aw.com/wps/media/objects/2228/2281678/Data_Sets/ASCII/Women13.txt", > header=TRUE) > # I HAVE CHANGED THE NAMES OF THE VARIABLES > # J is changed to "work" > # M is changed to "mar" > # S is changed to "school" > > attach(women) > > # THE VARIABLES OF USE ARE > # work: binary dependent variable > # mar: whether married or not > # school: years of schooling > > mu.start <- c(3, -1.5, 10) > data <- cbind(1, mar, school) > out <- nlm(mlogl, mu.start, y=work, X=data) > cat("Results", "\n") > out$estimate > > detach(women) > > ************************************* > > When I try to run the code, this is what I get: > > > source("probit.R") > Results > Warning messages: > 1: NA/Inf replaced by maximum positive value > 2: NA/Inf replaced by maximum positive value > 3: NA/Inf replaced by maximum positive value > 4: NA/Inf replaced by maximum positive value > > Thanks in advance. > Deepankar > > ______________________________________________ > R-help@stat.math.ethz.ch 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. > ______________________________________________ R-help@stat.math.ethz.ch 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.