Re: [R] Problems in programming a simple likelihood 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"  To:  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@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 ______________________________________________ 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.