[R] Problems in programming a simple likelihood

[Deepankar Basu]

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