Dimitris, Thanks a lot for your suggestion and also for suggestions that others have provided. I am learning fast and with the help of the R community will be able to get this going pretty soon. Of course, right now I am just trying to learn the language; so I am trying to program a standard probit model for which I know the answers. I could easily get the estimates with "glm". But I want to program the probit model to make sure I understand the subtleties of R.
I tried running the alternate script that you provided, but it is still not working. Am I making some mistake? Here is what I get when I run your script (which shows that the maximum number of iterations was reached without convergence): > source("probit1.R") > summary(out) Length Class Mode par 3 -none- numeric value 1 -none- numeric counts 2 -none- numeric convergence 1 -none- numeric message 0 -none- NULL Here is the script (exactly what you had suggested): 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(-3, -1.5, 0.5) 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 Thanks. Deepankar On Thu, 2007-04-19 at 09:26 +0200, Dimitris Rizopoulos wrote: > 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.