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" <[EMAIL PROTECTED]> To: <r-help@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@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.