[R] Problems in programming a simple likelihood
Deepankar Basu
basu.15 at osu.edu
Thu Apr 19 20:29:32 CEST 2007
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: rvaradhan at jhmi.edu
>
> Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
>
>
>
> ----------------------------------------------------------------------------
> --------
>
> -----Original Message-----
> From: r-help-bounces at stat.math.ethz.ch
> [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Deepankar Basu
> Sent: Thursday, April 19, 2007 12:42 PM
> To: Dimitris Rizopoulos
> Cc: r-help at 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" <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)
> > > }
> > >
> > > 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 at 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 at 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.
More information about the R-help
mailing list