Ravi, Thanks a lot for that clarification. Now I see that the code works.
Deepankar On Thu, 2007-04-19 at 14:01 -0400, Ravi Varadhan wrote: > Hi Deepankar, > > Dimitris' code works just fine. Your problem is that the output of optim > does not have a corresponding "summary" method. Instead you should simply > type the name of the object returned by "optim" to look at the results. > > > out <- optim(mu.start, mlogl, method = "CG", y = women$J, X = cbind(1, > women$M, women$S)) > > out > $par > [1] -3.0612277 -1.4567141 0.3659251 > > $value > [1] 13.32251 > > $counts > function gradient > 357 101 > > $convergence > [1] 1 > > $message > NULL > > Hope this helps, > Ravi. > > ---------------------------------------------------------------------------- > ------- > > Ravi Varadhan, Ph.D. > > Assistant Professor, The Center on Aging and Health > > Division of Geriatric Medicine and Gerontology > > Johns Hopkins University > > Ph: (410) 502-2619 > > Fax: (410) 614-9625 > > Email: [EMAIL PROTECTED] > > Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html > > > > ---------------------------------------------------------------------------- > -------- > > -----Original Message----- > From: [EMAIL PROTECTED] > [mailto:[EMAIL PROTECTED] On Behalf Of Deepankar Basu > Sent: Thursday, April 19, 2007 12:42 PM > To: Dimitris Rizopoulos > Cc: r-help@stat.math.ethz.ch > Subject: Re: [R] Problems in programming a simple likelihood > > 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. ______________________________________________ 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.