[R] Problems in programming a simple likelihood
Stephen Weigand
weigand.stephen at gmail.com
Thu Apr 19 04:26:24 CEST 2007
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 <basu.15 at osu.edu> 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 at 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.
>
More information about the R-help
mailing list