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.
