# [R] Problems in programming a simple likelihood

Dimitris Rizopoulos dimitris.rizopoulos at med.kuleuven.be
Thu Apr 19 09:26:09 CEST 2007

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

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 =

I hope it helps.

Best,
Dimitris

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

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)
> }
>
>
> # THE DATA SET CAN BE ACCESSED HERE
> # women <-
> # 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
>
> Deepankar
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help