[R] Problems in programming a simple likelihood

[Dimitris Rizopoulos]

try the following:

mlogl <- function (mu, y, X) {
    zeta <- as.vector(X %*% mu)
    y.logic <- as.logical(y)
    lgLik <- numeric(length(y))
    lgLik[y.logic] <- pnorm(zeta[y.logic], log.p = TRUE)
    lgLik[!y.logic] <- pnorm(zeta[!y.logic], lower.tail = FALSE, log.p 
= TRUE)
    -sum(lgLik)
}

women <- 
read.table("http://wps.aw.com/wps/media/objects/2228/2281678/Data_Sets/ASCII/Women13.txt", 
header=TRUE)

mu.start <- c(0, -1.5, 0.01)
out <- optim(mu.start, mlogl, method = "BFGS", y = women$J, X = 
cbind(1, women$M, women$S))
out

glm.fit(x = cbind(1, women$M, women$S), y = women$J, family = 
binomial(link = "probit"))$coefficients


I hope it helps.

Best,
Dimitris

----
Dimitris Rizopoulos
Ph.D. Student
Biostatistical Centre
School of Public Health
Catholic University of Leuven

Address: Kapucijnenvoer 35, Leuven, Belgium
Tel: +32/(0)16/336899
Fax: +32/(0)16/337015
Web: http://med.kuleuven.be/biostat/
     http://www.student.kuleuven.be/~m0390867/dimitris.htm


----- Original Message ----- 
From: "Deepankar Basu" <basu.15 at osu.edu>
To: <r-help at stat.math.ethz.ch>
Sent: Thursday, April 19, 2007 12:38 AM
Subject: [R] Problems in programming a simple likelihood


> 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.
> 


Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm