9 Mar
2009
9 Mar
'09
8 p.m.
hi!
Quick (but fundamental question)
in exercise 1 gamma needs to be positive, so I re-parameterized...In
order to get the right MLE from optim() i didn't
forget to "re-parameterize" again...and i get the same result as the
one i found analytically. However, things change when i look for the
SE. Indeed, I get different answers analytically and with R (using the
hessian). I know this comes from the re-parameterization because I
don't have this issue when I change the whole function and do not
re-parameterize gamma. So my question is:
- if I re-parametterize, how do I apply the transformation to the
hessian to get the right result
how does that fit with the section notes that follows, why do we take
"pnorm" (this is the first transformation that is applied in the
ll.binom function) of "opt.1000 - 1.96*se" and not of "se" for
instance???
#binomial log-likelihood (N = # of trials for each observation)
ll.binom <- function(par, y, N){
# reparameterize pi; only search over [0,1]
p <- pnorm(par)
# log-likelihood
out <- sum(y*log(p) + (N-y)*log(1-p))
return(out)
}
# compare to wald ci
se <- sqrt(solve(-optim(par=2, fn=ll.binom, y=samp.1000, N=10,
method="BFGS",
control=list(fnscale=-1), hessian=T)$hessian)) #$
wald.ci <- c( pnorm(opt.1000 - 1.96*se), pnorm(opt.1000 + 1.96*se))
wald.ci # 0.7364839 0.7535663
- if i do not re-parameterize in order to be done with it and have
both my analytical
and my R result fit, how can i justify i am not re-parameterizing gamma?!
thanks!
charlotte
9 Mar
9 Mar
8:46 p.m.
I've got the same problem with re-parametrization and then I found out that instead
re-parametrization in this particular case we are advised to use L-BFGS-B method in
optim(), so we can set a lower and an upper limits for out point of interest. So, the
point of interest will never get out of this interval, which is sort of
reparameterization, right?
The section notes says that log-likelihood (I guess it's in general for any
LL-function?) is invariant to re-parametrization (any?), but I was not much successful
applying this technique for exponential log-likelihood function. With re-parametrization
done exactly the way it's in section note I get totally different result from what I
derived analytically.
As for pnorm(opt.1000 - 1.96*se) I think that it could be done separately for se and
opt.1000 - they both are found for reparameterized function. So, we basically looking for
opt and se in different coordinate system and then re-parameterized to get them back to
our original coordinate system. Am I right?
Sincerely,
Olena Ageyeva
Date: Mon, 9 Mar 2009 20:00:14 -0400
From: charlotte.cavaille at gmail.com
To: gov2001-l at lists.fas.harvard.edu
Subject: [gov2001-l] EXO 1
hi!
Quick (but fundamental question)
in exercise 1 gamma needs to be positive, so I re-parameterized...In
order to get the right MLE from optim() i didn't
forget to "re-parameterize" again...and i get the same result as the
one i found analytically. However, things change when i look for the
SE. Indeed, I get different answers analytically and with R (using the
hessian). I know this comes from the re-parameterization because I
don't have this issue when I change the whole function and do not
re-parameterize gamma. So my question is:
- if I re-parametterize, how do I apply the transformation to the
hessian to get the right result
how does that fit with the section notes that follows, why do we take
"pnorm" (this is the first transformation that is applied in the
ll.binom function) of "opt.1000 - 1.96*se" and not of "se" for
instance???
#binomial log-likelihood (N = # of trials for each observation)
ll.binom <- function(par, y, N){
# reparameterize pi; only search over [0,1]
p <- pnorm(par)
# log-likelihood
out <- sum(y*log(p) + (N-y)*log(1-p))
return(out)
}
# compare to wald ci
se <- sqrt(solve(-optim(par=2, fn=ll.binom, y=samp.1000, N=10,
method="BFGS",
control=list(fnscale=-1), hessian=T)$hessian)) #$
wald.ci <- c( pnorm(opt.1000 - 1.96*se), pnorm(opt.1000 + 1.96*se))
wald.ci # 0.7364839 0.7535663
- if i do not re-parameterize in order to be done with it and have
both my analytical
and my R result fit, how can i justify i am not re-parameterizing gamma?!
thanks!
charlotte
_______________________________________________
gov2001-l mailing list
gov2001-l at lists.fas.harvard.edu
http://lists.fas.harvard.edu/mailman/listinfo/gov2001-l
_________________________________________________________________
Windows Live? Groups: Create an online spot for your favorite groups to meet.
http://windowslive.com/online/groups?ocid=TXT_TAGLM_WL_groups_032009
10:44 p.m.
I agree with Olena in that you need to use L-BFGS-B to get over that. My question is
more of a "good-standard" practice:
I notice that the problem to the reparametrization trick is that the SE's gets messed
up when we reparamterize - which seems to be a bit kludgy since usually we care about the
mean and the SE's. If reparamterization is a standard technique that most
statistician use, then is there a standard technique with BFGS that helps us find the SE?
(e.g. without specifying the bounds using L-BFGS-B)
Any insights?
Thanks,
Clarence
From: gov2001-l-bounces at lists.fas.harvard.edu [mailto:gov2001-l-bounces at
lists.fas.harvard.edu] On Behalf Of Olena Ageyeva
Sent: Monday, March 09, 2009 8:46 PM
To: gov2001-l at lists.fas.harvard.edu
Subject: Re: [gov2001-l] EXO 1
I've got the same problem with re-parametrization and then I found out that instead
re-parametrization in this particular case we are advised to use L-BFGS-B method in
optim(), so we can set a lower and an upper limits for out point of interest. So, the
point of interest will never get out of this interval, which is sort of
reparameterization, right?
The section notes says that log-likelihood (I guess it's in general for any
LL-function?) is invariant to re-parametrization (any?), but I was not much successful
applying this technique for exponential log-likelihood function. With re-parametrization
done exactly the way it's in section note I get totally different result from what I
derived analytically.
As for pnorm(opt.1000 - 1.96*se) I think that it could be done separately for se and
opt.1000 - they both are found for reparameterized function. So, we basically looking for
opt and se in different coordinate system and then re-parameterized to get them back to
our original coordinate system. Am I right?
Sincerely,
Olena Ageyeva
Date: Mon, 9 Mar 2009 20:00:14 -0400
From: charlotte.cavaille at gmail.com
To: gov2001-l at lists.fas.harvard.edu
Subject: [gov2001-l] EXO 1
hi!
Quick (but fundamental question)
in exercise 1 gamma needs to be positive, so I re-parameterized...In
order to get the right MLE from optim() i didn't
forget to "re-parameterize" again...and i get the same result as the
one i found analytically. However, things change when i look for the
SE. Indeed, I get different answers analytically and with R (using the
hessian). I know this comes from the re-parameterization because I
don't have this issue when I change the whole function and do not
re-parameterize gamma. So my question is:
- if I re-parametterize, how do I apply the transformation to the
hessian to get the right result
how does that fit with the section notes that follows, why do we take
"pnorm" (this is the first transformation that is applied in the
ll.binom function) of "opt.1000 - 1.96*se" and not of "se" for
instance???
#binomial log-likelihood (N = # of trials for each observation)
ll.binom <- function(par, y, N){
# reparameterize pi; only search over [0,1]
p <- pnorm(par)
# log-likelihood
out <- sum(y*log(p) + (N-y)*log(1-p))
return(out)
}
# compare to wald ci
se <- sqrt(solve(-optim(par=2, fn=ll.binom, y=samp.1000, N=10,
method="BFGS",
control=list(fnscale=-1), hessian=T)$hessian)) #$
wald.ci <- c( pnorm(opt.1000 - 1.96*se), pnorm(opt.1000 + 1.96*se))
wald.ci # 0.7364839 0.7535663
- if i do not re-parameterize in order to be done with it and have
both my analytical
and my R result fit, how can i justify i am not re-parameterizing gamma?!
thanks!
charlotte
_______________________________________________
gov2001-l mailing list
gov2001-l at lists.fas.harvard.edu
http://lists.fas.harvard.edu/mailman/listinfo/gov2001-l
________________________________
Windows Live(tm) Groups: Create an online spot for your favorite groups to meet. Check it
out.<http://windowslive.com/online/groups?ocid=TXT_TAGLM_WL_groups_03200…
10 Mar
10 Mar
1:06 a.m.
Some clarification:
(1) Yes, MLEs are invariant to re-parameterization (i.e. going from an
bounded to an unbounded and back to an bounded scale when estimating your
parameters). But, measures of uncertainty are not invariant to
re-parameterization. There are two main ways to get your standard errors
(and confidence intervals) back to the original scale after you have
estimated your model by re-parameterizing your parameters. (1) The delta
method -- this is an analytical method to get the variance for our
parameters back to the original scale. Basically, it involves some math to
convert what we know about the distribution of our estimates on the
unbounded scale -- that is that they are distributed N(theta, Var(theta)) --
to get to what we want to know about the distribution of our estimates on
the bounded scale -- that is that they are distributed N(g(theta),
g(var(theta))). (2) The simpler way is to just simulate our parameters. We
will take repeated draws of our parameters from their distribution on the
unbounded scale (i.e. N(theta, sigma.sq) ) and convert each of these draws
back to the bounded scale. We can then simply calculate our standard errors
and confidence intervals by using this simulated sampling distribution.
(2) In PS 5 question 1.2.b, we recommend you use the "L-BFGS-B" method to
get your estimates on the original bounded scale. However, if you are
interested in the delta method, we encourage you to learn it and use it to
answer this question.
Best,
Miya
On Mon, Mar 9, 2009 at 9:44 PM, Lee, Clarence <clee at hbs.edu> wrote:
I agree with Olena in that you need to use L-BFGS-B
to get over that.
My question is more of a ?good-standard? practice:
I notice that the problem to the reparametrization trick is that the SE?s
gets messed up when we reparamterize ? which seems to be a bit kludgy since
usually we care about the mean and the SE?s. If reparamterization is a
standard technique that most statistician use, then is there a standard
technique with BFGS that helps us find the SE? (e.g. without specifying the
bounds using L-BFGS-B)
Any insights?
Thanks,
Clarence
*From:* gov2001-l-bounces at lists.fas.harvard.edu [mailto:
gov2001-l-bounces at lists.fas.harvard.edu] *On Behalf Of *Olena Ageyeva
*Sent:* Monday, March 09, 2009 8:46 PM
*To:* gov2001-l at lists.fas.harvard.edu
*Subject:* Re: [gov2001-l] EXO 1
I've got the same problem with re-parametrization and then I found out that
instead re-parametrization in this particular case we are advised to use
L-BFGS-B method in optim(), so we can set a lower and an upper limits for
out point of interest. So, the point of interest will never get out of this
interval, which is sort of reparameterization, right?
The section notes says that log-likelihood (I guess it's in general for any
LL-function?) is invariant to re-parametrization (any?), but I was not much
successful applying this technique for exponential log-likelihood function.
With re-parametrization done exactly the way it's in section note I get
totally different result from what I derived analytically.
As for pnorm(opt.1000 - 1.96*se) I think that it could be done separately
for se and opt.1000 - they both are found for reparameterized function. So,
we basically looking for opt and se in different coordinate system and then
re-parameterized to get them back to our original coordinate system. Am I
right?
Sincerely,
Olena Ageyeva
--
Miya Woolfalk
Ph.D. Student
Harvard University
Government and Social Policy
Date: Mon, 9 Mar 2009 20:00:14 -0400
From: charlotte.cavaille at gmail.com
To: gov2001-l at lists.fas.harvard.edu
Subject: [gov2001-l] EXO 1
hi!
Quick (but fundamental question)
in exercise 1 gamma needs to be positive, so I re-parameterized...In
order to get the right MLE from optim() i didn't
forget to "re-parameterize" again...and i get the same result as the
one i found analytically. However, things change when i look for the
SE. Indeed, I get different answers analytically and with R (using the
hessian). I know this comes from the re-parameterization because I
don't have this issue when I change the whole function and do not
re-parameterize gamma. So my question is:
- if I re-parametterize, how do I apply the transformation to the
hessian to get the right result
how does that fit with the section notes that follows, why do we take
"pnorm" (this is the first transformation that is applied in the
ll.binom function) of "opt.1000 - 1.96*se" and not of "se" for
instance???
#binomial log-likelihood (N = # of trials for each observation)
ll.binom <- function(par, y, N){
# reparameterize pi; only search over [0,1]
p <- pnorm(par)
# log-likelihood
out <- sum(y*log(p) + (N-y)*log(1-p))
return(out)
}
# compare to wald ci
se <- sqrt(solve(-optim(par=2, fn=ll.binom, y=samp.1000, N=10,
method="BFGS",
control=list(fnscale=-1), hessian=T)$hessian))
#$
wald.ci <- c( pnorm(opt.1000 - 1.96*se), pnorm(opt.1000 + 1.96*se))
wald.ci # 0.7364839 0.7535663
- if i do not re-parameterize in order to be done with it and have
both my analytical
and my R result fit, how can i justify i am not re-parameterizing gamma?!
thanks!
charlotte
_______________________________________________
gov2001-l mailing list
gov2001-l at lists.fas.harvard.edu
http://lists.fas.harvard.edu/mailman/listinfo/gov2001-l
------------------------------
Windows Live? Groups: Create an online spot for your favorite groups to
meet. Check it
out.<http://windowslive.com/online/groups?ocid=TXT_TAGLM_WL_groups_03200…
_______________________________________________
gov2001-l mailing list
gov2001-l at lists.fas.harvard.edu
http://lists.fas.harvard.edu/mailman/listinfo/gov2001-l
2:26 p.m.
New subject: [gov2001-l] Syntax to specify bounds in using L-BFGS-B?
What is the syntax to specify the bounds using L-BFGS-B. I've looked through the
R-help file but can't figure out how I tell optim what bounds to use....
From: gov2001-l-bounces at lists.fas.harvard.edu [mailto:gov2001-l-bounces at
lists.fas.harvard.edu] On Behalf Of Lee, Clarence
Sent: Monday, March 09, 2009 10:44 PM
To: gov2001-l at lists.fas.harvard.edu
Subject: Re: [gov2001-l] EXO 1
I agree with Olena in that you need to use L-BFGS-B to get over that. My question is
more of a "good-standard" practice:
I notice that the problem to the reparametrization trick is that the SE's gets messed
up when we reparamterize - which seems to be a bit kludgy since usually we care about the
mean and the SE's. If reparamterization is a standard technique that most
statistician use, then is there a standard technique with BFGS that helps us find the SE?
(e.g. without specifying the bounds using L-BFGS-B)
Any insights?
Thanks,
Clarence
From: gov2001-l-bounces at lists.fas.harvard.edu [mailto:gov2001-l-bounces at
lists.fas.harvard.edu] On Behalf Of Olena Ageyeva
Sent: Monday, March 09, 2009 8:46 PM
To: gov2001-l at lists.fas.harvard.edu
Subject: Re: [gov2001-l] EXO 1
I've got the same problem with re-parametrization and then I found out that instead
re-parametrization in this particular case we are advised to use L-BFGS-B method in
optim(), so we can set a lower and an upper limits for out point of interest. So, the
point of interest will never get out of this interval, which is sort of
reparameterization, right?
The section notes says that log-likelihood (I guess it's in general for any
LL-function?) is invariant to re-parametrization (any?), but I was not much successful
applying this technique for exponential log-likelihood function. With re-parametrization
done exactly the way it's in section note I get totally different result from what I
derived analytically.
As for pnorm(opt.1000 - 1.96*se) I think that it could be done separately for se and
opt.1000 - they both are found for reparameterized function. So, we basically looking for
opt and se in different coordinate system and then re-parameterized to get them back to
our original coordinate system. Am I right?
Sincerely,
Olena Ageyeva
Date: Mon, 9 Mar 2009 20:00:14 -0400
From: charlotte.cavaille at gmail.com
To: gov2001-l at lists.fas.harvard.edu
Subject: [gov2001-l] EXO 1
hi!
Quick (but fundamental question)
in exercise 1 gamma needs to be positive, so I re-parameterized...In
order to get the right MLE from optim() i didn't
forget to "re-parameterize" again...and i get the same result as the
one i found analytically. However, things change when i look for the
SE. Indeed, I get different answers analytically and with R (using the
hessian). I know this comes from the re-parameterization because I
don't have this issue when I change the whole function and do not
re-parameterize gamma. So my question is:
- if I re-parametterize, how do I apply the transformation to the
hessian to get the right result
how does that fit with the section notes that follows, why do we take
"pnorm" (this is the first transformation that is applied in the
ll.binom function) of "opt.1000 - 1.96*se" and not of "se" for
instance???
#binomial log-likelihood (N = # of trials for each observation)
ll.binom <- function(par, y, N){
# reparameterize pi; only search over [0,1]
p <- pnorm(par)
# log-likelihood
out <- sum(y*log(p) + (N-y)*log(1-p))
return(out)
}
# compare to wald ci
se <- sqrt(solve(-optim(par=2, fn=ll.binom, y=samp.1000, N=10,
method="BFGS",
control=list(fnscale=-1), hessian=T)$hessian)) #$
wald.ci <- c( pnorm(opt.1000 - 1.96*se), pnorm(opt.1000 + 1.96*se))
wald.ci # 0.7364839 0.7535663
- if i do not re-parameterize in order to be done with it and have
both my analytical
and my R result fit, how can i justify i am not re-parameterizing gamma?!
thanks!
charlotte
_______________________________________________
gov2001-l mailing list
gov2001-l at lists.fas.harvard.edu
http://lists.fas.harvard.edu/mailman/listinfo/gov2001-l
________________________________
Windows Live(tm) Groups: Create an online spot for your favorite groups to meet. Check it
out.<http://windowslive.com/online/groups?ocid=TXT_TAGLM_WL_groups_03200…
2:49 p.m.
New subject: [gov2001-l] Syntax to specify bounds in using L-BFGS-B?
in the optim help file under "Arguments", there are two arguments that you
are looking for.
On Tue, Mar 10, 2009 at 2:26 PM, Allen, Abigail <aallen at hbs.edu> wrote:
What is the syntax to specify the bounds using
L-BFGS-B. I?ve looked
through the R-help file but can?t figure out how I tell optim what bounds to
use?.
*From:* gov2001-l-bounces at lists.fas.harvard.edu [mailto:
gov2001-l-bounces at lists.fas.harvard.edu] *On Behalf Of *Lee, Clarence
*Sent:* Monday, March 09, 2009 10:44 PM
*To:* gov2001-l at lists.fas.harvard.edu
*Subject:* Re: [gov2001-l] EXO 1
I agree with Olena in that you need to use L-BFGS-B to get over that. My
question is more of a ?good-standard? practice:
I notice that the problem to the reparametrization trick is that the SE?s
gets messed up when we reparamterize ? which seems to be a bit kludgy since
usually we care about the mean and the SE?s. If reparamterization is a
standard technique that most statistician use, then is there a standard
technique with BFGS that helps us find the SE? (e.g. without specifying the
bounds using L-BFGS-B)
Any insights?
Thanks,
Clarence
*From:* gov2001-l-bounces at lists.fas.harvard.edu [mailto:
gov2001-l-bounces at lists.fas.harvard.edu] *On Behalf Of *Olena Ageyeva
*Sent:* Monday, March 09, 2009 8:46 PM
*To:* gov2001-l at lists.fas.harvard.edu
*Subject:* Re: [gov2001-l] EXO 1
I've got the same problem with re-parametrization and then I found out that
instead re-parametrization in this particular case we are advised to use
L-BFGS-B method in optim(), so we can set a lower and an upper limits for
out point of interest. So, the point of interest will never get out of this
interval, which is sort of reparameterization, right?
The section notes says that log-likelihood (I guess it's in general for any
LL-function?) is invariant to re-parametrization (any?), but I was not much
successful applying this technique for exponential log-likelihood function.
With re-parametrization done exactly the way it's in section note I get
totally different result from what I derived analytically.
As for pnorm(opt.1000 - 1.96*se) I think that it could be done separately
for se and opt.1000 - they both are found for reparameterized function. So,
we basically looking for opt and se in different coordinate system and then
re-parameterized to get them back to our original coordinate system. Am I
right?
Sincerely,
Olena Ageyeva
--
Patrick Lam
Department of Government and Institute for Quantitative Social Science,
Harvard University
http://www.people.fas.harvard.edu/~plam
Date: Mon, 9 Mar 2009 20:00:14 -0400
From: charlotte.cavaille at gmail.com
To: gov2001-l at lists.fas.harvard.edu
Subject: [gov2001-l] EXO 1
hi!
Quick (but fundamental question)
in exercise 1 gamma needs to be positive, so I re-parameterized...In
order to get the right MLE from optim() i didn't
forget to "re-parameterize" again...and i get the same result as the
one i found analytically. However, things change when i look for the
SE. Indeed, I get different answers analytically and with R (using the
hessian). I know this comes from the re-parameterization because I
don't have this issue when I change the whole function and do not
re-parameterize gamma. So my question is:
- if I re-parametterize, how do I apply the transformation to the
hessian to get the right result
how does that fit with the section notes that follows, why do we take
"pnorm" (this is the first transformation that is applied in the
ll.binom function) of "opt.1000 - 1.96*se" and not of "se" for
instance???
#binomial log-likelihood (N = # of trials for each observation)
ll.binom <- function(par, y, N){
# reparameterize pi; only search over [0,1]
p <- pnorm(par)
# log-likelihood
out <- sum(y*log(p) + (N-y)*log(1-p))
return(out)
}
# compare to wald ci
se <- sqrt(solve(-optim(par=2, fn=ll.binom, y=samp.1000, N=10,
method="BFGS",
control=list(fnscale=-1), hessian=T)$hessian))
#$
wald.ci <- c( pnorm(opt.1000 - 1.96*se), pnorm(opt.1000 + 1.96*se))
wald.ci # 0.7364839 0.7535663
- if i do not re-parameterize in order to be done with it and have
both my analytical
and my R result fit, how can i justify i am not re-parameterizing gamma?!
thanks!
charlotte
_______________________________________________
gov2001-l mailing list
gov2001-l at lists.fas.harvard.edu
http://lists.fas.harvard.edu/mailman/listinfo/gov2001-l
------------------------------
Windows Live? Groups: Create an online spot for your favorite groups to
meet. Check it
out.<http://windowslive.com/online/groups?ocid=TXT_TAGLM_WL_groups_03200…
_______________________________________________
gov2001-l mailing list
gov2001-l at lists.fas.harvard.edu
http://lists.fas.harvard.edu/mailman/listinfo/gov2001-l
3:10 p.m.
New subject: [gov2001-l] Syntax to specify bounds in using L-BFGS-B?
Ah thanks... can't believe I missed those :)
From: gov2001-l-bounces at lists.fas.harvard.edu [mailto:gov2001-l-bounces at
lists.fas.harvard.edu] On Behalf Of Patrick Lam
Sent: Tuesday, March 10, 2009 2:50 PM
To: gov2001-l at lists.fas.harvard.edu
Subject: Re: [gov2001-l] Syntax to specify bounds in using L-BFGS-B?
in the optim help file under "Arguments", there are two arguments that you are
looking for.
On Tue, Mar 10, 2009 at 2:26 PM, Allen, Abigail <aallen at hbs.edu<mailto:aallen at
hbs.edu>> wrote:
What is the syntax to specify the bounds using L-BFGS-B. I've looked through the
R-help file but can't figure out how I tell optim what bounds to use....
From: gov2001-l-bounces at lists.fas.harvard.edu<mailto:gov2001-l-bounces at
lists.fas.harvard.edu> [mailto:gov2001-l-bounces at
lists.fas.harvard.edu<mailto:gov2001-l-bounces at lists.fas.harvard.edu>] On Behalf
Of Lee, Clarence
Sent: Monday, March 09, 2009 10:44 PM
To: gov2001-l at lists.fas.harvard.edu<mailto:gov2001-l at lists.fas.harvard.edu>
Subject: Re: [gov2001-l] EXO 1
I agree with Olena in that you need to use L-BFGS-B to get over that. My question is
more of a "good-standard" practice:
I notice that the problem to the reparametrization trick is that the SE's gets messed
up when we reparamterize - which seems to be a bit kludgy since usually we care about the
mean and the SE's. If reparamterization is a standard technique that most
statistician use, then is there a standard technique with BFGS that helps us find the SE?
(e.g. without specifying the bounds using L-BFGS-B)
Any insights?
Thanks,
Clarence
From: gov2001-l-bounces at lists.fas.harvard.edu<mailto:gov2001-l-bounces at
lists.fas.harvard.edu> [mailto:gov2001-l-bounces at
lists.fas.harvard.edu<mailto:gov2001-l-bounces at lists.fas.harvard.edu>] On Behalf
Of Olena Ageyeva
Sent: Monday, March 09, 2009 8:46 PM
To: gov2001-l at lists.fas.harvard.edu<mailto:gov2001-l at lists.fas.harvard.edu>
Subject: Re: [gov2001-l] EXO 1
I've got the same problem with re-parametrization and then I found out that instead
re-parametrization in this particular case we are advised to use L-BFGS-B method in
optim(), so we can set a lower and an upper limits for out point of interest. So, the
point of interest will never get out of this interval, which is sort of
reparameterization, right?
The section notes says that log-likelihood (I guess it's in general for any
LL-function?) is invariant to re-parametrization (any?), but I was not much successful
applying this technique for exponential log-likelihood function. With re-parametrization
done exactly the way it's in section note I get totally different result from what I
derived analytically.
As for pnorm(opt.1000 - 1.96*se) I think that it could be done separately for se and
opt.1000 - they both are found for reparameterized function. So, we basically looking for
opt and se in different coordinate system and then re-parameterized to get them back to
our original coordinate system. Am I right?
Sincerely,
Olena Ageyeva
Date: Mon, 9 Mar 2009 20:00:14 -0400
From: charlotte.cavaille at gmail.com<mailto:charlotte.cavaille at gmail.com>
To: gov2001-l at lists.fas.harvard.edu<mailto:gov2001-l at lists.fas.harvard.edu>
Subject: [gov2001-l] EXO 1
hi!
Quick (but fundamental question)
in exercise 1 gamma needs to be positive, so I re-parameterized...In
order to get the right MLE from optim() i didn't
forget to "re-parameterize" again...and i get the same result as the
one i found analytically. However, things change when i look for the
SE. Indeed, I get different answers analytically and with R (using the
hessian). I know this comes from the re-parameterization because I
don't have this issue when I change the whole function and do not
re-parameterize gamma. So my question is:
- if I re-parametterize, how do I apply the transformation to the
hessian to get the right result
how does that fit with the section notes that follows, why do we take
"pnorm" (this is the first transformation that is applied in the
ll.binom function) of "opt.1000 - 1.96*se" and not of "se" for
instance???
#binomial log-likelihood (N = # of trials for each observation)
ll.binom <- function(par, y, N){
# reparameterize pi; only search over [0,1]
p <- pnorm(par)
# log-likelihood
out <- sum(y*log(p) + (N-y)*log(1-p))
return(out)
}
# compare to wald ci
se <- sqrt(solve(-optim(par=2, fn=ll.binom, y=samp.1000, N=10,
method="BFGS",
control=list(fnscale=-1), hessian=T)$hessian)) #$
wald.ci<http://wald.ci> <- c( pnorm(opt.1000 - 1.96*se), pnorm(opt.1000 +
1.96*se))
wald.ci<http://wald.ci> # 0.7364839 0.7535663
- if i do not re-parameterize in order to be done with it and have
both my analytical
and my R result fit, how can i justify i am not re-parameterizing gamma?!
thanks!
charlotte
_______________________________________________
gov2001-l mailing list
gov2001-l at lists.fas.harvard.edu<mailto:gov2001-l at lists.fas.harvard.edu>
http://lists.fas.harvard.edu/mailman/listinfo/gov2001-l
________________________________
Windows Live(tm) Groups: Create an online spot for your favorite groups to meet. Check it
out.<http://windowslive.com/online/groups?ocid=TXT_TAGLM_WL_groups_03200…
_______________________________________________
gov2001-l mailing list
gov2001-l at lists.fas.harvard.edu<mailto:gov2001-l at lists.fas.harvard.edu>
http://lists.fas.harvard.edu/mailman/listinfo/gov2001-l
--
Patrick Lam
Department of Government and Institute for Quantitative Social Science, Harvard
University
http://www.people.fas.harvard.edu/~plam
11 Mar
11 Mar
12:07 a.m.
New subject: [gov2001-l] problems with log(exp().... Sequential processing in R causing issues
In attempting to debug my function I came across the below anomaly which I cannot figure
out.
Whenever I try to evaluate log(exp(A)) for any A >= 710 R evaluates the result as
"Inf" which is clearly not the case. See R output below:
log(exp(708))
[1] 708
log(exp(709))
[1] 709
log(exp(710))
[1] Inf
log(exp(832))
[1] Inf
I assume that this happens because R evaluates the function inside out: i.e. first tries
to take the exp(710) which produces more digits than R will hold: hence R evaluates
exp(710) as Inf, and then applies the log fxn to this result--> log(Inf) = Inf.
Ok logical. However, incredibly inconvenient for purposes of this homework set, and for
optim() in general, as many negative values of B result in ln(1+e^(-xiB)) > ln(1+e^710)
--> creating an "-Inf" evaluation of the log likelihood. This sequential
calculation jams my function when I try to optimize since R cannot sum (Inf+510+712) for
example --> Error in -sum(formula in here :)) : invalid 'type' (character) of
argument.
Has anyone else encountered this issue?
Any suggestions?
Thanks,
Abigail
From: gov2001-l-bounces at lists.fas.harvard.edu [mailto:gov2001-l-bounces at
lists.fas.harvard.edu] On Behalf Of Patrick Lam
Sent: Tuesday, March 10, 2009 2:50 PM
To: gov2001-l at lists.fas.harvard.edu
Subject: Re: [gov2001-l] Syntax to specify bounds in using L-BFGS-B?
in the optim help file under "Arguments", there are two arguments that you are
looking for.
On Tue, Mar 10, 2009 at 2:26 PM, Allen, Abigail <aallen at hbs.edu<mailto:aallen at
hbs.edu>> wrote:
What is the syntax to specify the bounds using L-BFGS-B. I've looked through the
R-help file but can't figure out how I tell optim what bounds to use....
From: gov2001-l-bounces at lists.fas.harvard.edu<mailto:gov2001-l-bounces at
lists.fas.harvard.edu> [mailto:gov2001-l-bounces at
lists.fas.harvard.edu<mailto:gov2001-l-bounces at lists.fas.harvard.edu>] On Behalf
Of Lee, Clarence
Sent: Monday, March 09, 2009 10:44 PM
To: gov2001-l at lists.fas.harvard.edu<mailto:gov2001-l at lists.fas.harvard.edu>
Subject: Re: [gov2001-l] EXO 1
I agree with Olena in that you need to use L-BFGS-B to get over that. My question is
more of a "good-standard" practice:
I notice that the problem to the reparametrization trick is that the SE's gets messed
up when we reparamterize - which seems to be a bit kludgy since usually we care about the
mean and the SE's. If reparamterization is a standard technique that most
statistician use, then is there a standard technique with BFGS that helps us find the SE?
(e.g. without specifying the bounds using L-BFGS-B)
Any insights?
Thanks,
Clarence
From: gov2001-l-bounces at lists.fas.harvard.edu<mailto:gov2001-l-bounces at
lists.fas.harvard.edu> [mailto:gov2001-l-bounces at
lists.fas.harvard.edu<mailto:gov2001-l-bounces at lists.fas.harvard.edu>] On Behalf
Of Olena Ageyeva
Sent: Monday, March 09, 2009 8:46 PM
To: gov2001-l at lists.fas.harvard.edu<mailto:gov2001-l at lists.fas.harvard.edu>
Subject: Re: [gov2001-l] EXO 1
I've got the same problem with re-parametrization and then I found out that instead
re-parametrization in this particular case we are advised to use L-BFGS-B method in
optim(), so we can set a lower and an upper limits for out point of interest. So, the
point of interest will never get out of this interval, which is sort of
reparameterization, right?
The section notes says that log-likelihood (I guess it's in general for any
LL-function?) is invariant to re-parametrization (any?), but I was not much successful
applying this technique for exponential log-likelihood function. With re-parametrization
done exactly the way it's in section note I get totally different result from what I
derived analytically.
As for pnorm(opt.1000 - 1.96*se) I think that it could be done separately for se and
opt.1000 - they both are found for reparameterized function. So, we basically looking for
opt and se in different coordinate system and then re-parameterized to get them back to
our original coordinate system. Am I right?
Sincerely,
Olena Ageyeva
Date: Mon, 9 Mar 2009 20:00:14 -0400
From: charlotte.cavaille at gmail.com<mailto:charlotte.cavaille at gmail.com>
To: gov2001-l at lists.fas.harvard.edu<mailto:gov2001-l at lists.fas.harvard.edu>
Subject: [gov2001-l] EXO 1
hi!
Quick (but fundamental question)
in exercise 1 gamma needs to be positive, so I re-parameterized...In
order to get the right MLE from optim() i didn't
forget to "re-parameterize" again...and i get the same result as the
one i found analytically. However, things change when i look for the
SE. Indeed, I get different answers analytically and with R (using the
hessian). I know this comes from the re-parameterization because I
don't have this issue when I change the whole function and do not
re-parameterize gamma. So my question is:
- if I re-parametterize, how do I apply the transformation to the
hessian to get the right result
how does that fit with the section notes that follows, why do we take
"pnorm" (this is the first transformation that is applied in the
ll.binom function) of "opt.1000 - 1.96*se" and not of "se" for
instance???
#binomial log-likelihood (N = # of trials for each observation)
ll.binom <- function(par, y, N){
# reparameterize pi; only search over [0,1]
p <- pnorm(par)
# log-likelihood
out <- sum(y*log(p) + (N-y)*log(1-p))
return(out)
}
# compare to wald ci
se <- sqrt(solve(-optim(par=2, fn=ll.binom, y=samp.1000, N=10,
method="BFGS",
control=list(fnscale=-1), hessian=T)$hessian)) #$
wald.ci<http://wald.ci> <- c( pnorm(opt.1000 - 1.96*se), pnorm(opt.1000 +
1.96*se))
wald.ci<http://wald.ci> # 0.7364839 0.7535663
- if i do not re-parameterize in order to be done with it and have
both my analytical
and my R result fit, how can i justify i am not re-parameterizing gamma?!
thanks!
charlotte
_______________________________________________
gov2001-l mailing list
gov2001-l at lists.fas.harvard.edu<mailto:gov2001-l at lists.fas.harvard.edu>
http://lists.fas.harvard.edu/mailman/listinfo/gov2001-l
________________________________
Windows Live(tm) Groups: Create an online spot for your favorite groups to meet. Check it
out.<http://windowslive.com/online/groups?ocid=TXT_TAGLM_WL_groups_03200…
_______________________________________________
gov2001-l mailing list
gov2001-l at lists.fas.harvard.edu<mailto:gov2001-l at lists.fas.harvard.edu>
http://lists.fas.harvard.edu/mailman/listinfo/gov2001-l
--
Patrick Lam
Department of Government and Institute for Quantitative Social Science, Harvard
University
http://www.people.fas.harvard.edu/~plam
1:02 p.m.
New subject: [gov2001-l] problems with log(exp().... Sequential processing in R causing issues
I would try to check that your variable doesn't get grater then 710.
This simple modification makes the algorithm work:
m=605
n=110
log(exp(m)*exp(n))
[1] Inf
log(exp(m))+log(exp(n))
[1] 715
I would just create a function for log(exp(x)), where x is a large number.
log_exp_for_large_numbers<-function(x)
{
x_array<-c()
while (x>0)
{
ifelse (x>710, x_array<-cbind(x_array,709), x_array<-cbind(x_array,x))
x<-x-709
}
out<-0
for (i in 1:length(x_array))
{
out<-out+log(exp(x_array[i]))
}
return(out)
}
log_exp_for_large_numbers(31568)
[1] 31568
Sincerely,
Olena Ageyeva
From: aallen at hbs.edu
To: gov2001-l at lists.fas.harvard.edu
Date: Wed, 11 Mar 2009 00:07:38 -0400
Subject: [gov2001-l] problems with log(exp().... Sequential processing in R causing
issues
In attempting to debug my function I came across the below anomaly
which I cannot figure out.
Whenever I try to evaluate log(exp(A)) for any A >= 710 R
evaluates the result as ?Inf? which is clearly not the case. See R
output below:
log(exp(708))
[1] 708
log(exp(709))
[1] 709
log(exp(710))
[1] Inf
log(exp(832))
[1] Inf
I assume that this happens because R evaluates the function
inside out: i.e. first tries to take the exp(710) which produces more digits
than R will hold: hence R evaluates exp(710) as Inf, and then applies the log fxn
to this result? log(Inf) = Inf.
Ok logical. However, incredibly inconvenient for purposes of
this homework set, and for optim() in general, as many negative values of B
result in ln(1+e^(-xiB)) > ln(1+e^710) ? creating an ?-Inf?
evaluation of the log likelihood. This sequential calculation jams my function
when I try to optimize since R cannot sum (Inf+510+712) for example ? Error
in -sum(formula in here J) : invalid 'type' (character) of
argument.
Has anyone else encountered this issue?
Any suggestions?
Thanks,
Abigail
From:
gov2001-l-bounces at lists.fas.harvard.edu [mailto:gov2001-l-bounces at
lists.fas.harvard.edu]
On Behalf Of Patrick Lam
Sent: Tuesday, March 10, 2009 2:50 PM
To: gov2001-l at lists.fas.harvard.edu
Subject: Re: [gov2001-l] Syntax to specify bounds in using L-BFGS-B?
in the optim help file under
"Arguments", there are two arguments that you are looking for.
On Tue, Mar 10, 2009 at 2:26 PM, Allen, Abigail <aallen at hbs.edu> wrote:
What is the syntax to specify
the bounds using L-BFGS-B. I?ve looked through the R-help file but
can?t figure out how I tell optim what bounds to use?.
From: gov2001-l-bounces at lists.fas.harvard.edu [mailto:gov2001-l-bounces at
lists.fas.harvard.edu]
On Behalf Of Lee, Clarence
Sent: Monday, March 09, 2009 10:44 PM
To: gov2001-l at lists.fas.harvard.edu
Subject: Re: [gov2001-l] EXO 1
I agree with Olena in that you
need to use L-BFGS-B to get over that. My question is more of a
?good-standard? practice:
I notice that the problem to
the reparametrization trick is that the SE?s gets messed up when we
reparamterize ? which seems to be a bit kludgy since usually we care
about the mean and the SE?s. If reparamterization is a
standard technique that most statistician use, then is there a standard
technique with BFGS that helps us find the SE? (e.g. without specifying the
bounds using L-BFGS-B)
Any insights?
Thanks,
Clarence
From: gov2001-l-bounces at lists.fas.harvard.edu [mailto:gov2001-l-bounces at
lists.fas.harvard.edu]
On Behalf Of Olena Ageyeva
Sent: Monday, March 09, 2009 8:46 PM
To: gov2001-l at lists.fas.harvard.edu
Subject: Re: [gov2001-l] EXO 1
I've got the
same problem with re-parametrization and then I found out that instead
re-parametrization in this particular case we are advised to use L-BFGS-B
method in optim(), so we can set a lower and an upper limits for out point of
interest. So, the point of interest will never get out of this interval, which
is sort of reparameterization, right?
The section notes says that log-likelihood (I guess it's in general for any
LL-function?) is invariant to re-parametrization (any?), but I was not much
successful applying this technique for exponential log-likelihood function.
With re-parametrization done exactly the way it's in section note I get totally
different result from what I derived analytically.
As for pnorm(opt.1000 - 1.96*se) I think that it could be done separately for
se and opt.1000 - they both are found for reparameterized function. So, we
basically looking for opt and se in different coordinate system and then
re-parameterized to get them back to our original coordinate system. Am I
right?
Sincerely,
Olena Ageyeva
Date: Mon, 9 Mar 2009 20:00:14 -0400
From: charlotte.cavaille at gmail.com
To: gov2001-l at lists.fas.harvard.edu
Subject: [gov2001-l] EXO 1
hi!
Quick (but fundamental question)
in exercise 1 gamma needs to be positive, so I
re-parameterized...In
order to get the right MLE from optim() i didn't
forget to "re-parameterize" again...and i
get the same result as
the
one i found analytically. However, things change when
i look for the
SE. Indeed, I get different answers analytically and
with R (using the
hessian). I know this comes from the
re-parameterization because I
don't have this issue when I change the whole
function and do not
re-parameterize gamma. So my question is:
- if I re-parametterize, how do I apply the
transformation to the
hessian to get the right result
how does that fit with the section notes that follows,
why do we take
"pnorm" (this is the first transformation
that is applied in the
ll.binom function) of "opt.1000 - 1.96*se"
and not of
"se" for
instance???
#binomial log-likelihood (N = # of trials for each
observation)
ll.binom <- function(par, y, N){
# reparameterize pi; only search over [0,1]
p <- pnorm(par)
# log-likelihood
out <- sum(y*log(p) + (N-y)*log(1-p))
return(out)
}
# compare to wald ci
se <- sqrt(solve(-optim(par=2, fn=ll.binom,
y=samp.1000, N=10,
method="BFGS",
control=list(fnscale=-1), hessian=T)$hessian)) #$
wald.ci <- c(
pnorm(opt.1000 - 1.96*se),
pnorm(opt.1000 + 1.96*se))
wald.ci # 0.7364839 0.7535663
- if i do not re-parameterize in order to be done with
it and have
both my analytical
and my R result fit, how can i justify i am not
re-parameterizing gamma?!
thanks!
charlotte
_______________________________________________
gov2001-l mailing list
gov2001-l at lists.fas.harvard.edu
Windows Live? Groups: Create an online
spot for your favorite groups to meet. Check it out.
_______________________________________________
gov2001-l mailing list
gov2001-l at lists.fas.harvard.edu
http://lists.fas.harvard.edu/mailman/listinfo/gov2001-l
--
Patrick Lam
Department of Government and Institute for Quantitative Social Science, Harvard
University
http://www.people.fas.harvard.edu/~plam
_________________________________________________________________
Windows Live? Contacts: Organize your contact list.
http://windowslive.com/connect/post/marcusatmicrosoft.spaces.live.com-Blog-…
1:06 p.m.
New subject: [gov2001-l] problems with log(exp().... Sequential processing in R causing issues
Abigail,
Before trying to implement a work around, you should make sure that your log
likelihood function is programed correctly. I have not heard of others
encountering this issue so my suspicion is that you have an error in your
function.
Best,
Miya
On Wed, Mar 11, 2009 at 1:02 PM, Olena Ageyeva <eageeva at hotmail.com> wrote:
I would try to check that your variable doesn't
get grater then 710.
This simple modification makes the algorithm work:
--
Miya Woolfalk
Ph.D. Student
Harvard University
Government and Social Policy
m=605
n=110
log(exp(m)*exp(n))
[1] Inf
log(exp(m))+log(exp(n))
[1] 715
I would just create a function for log(exp(x)), where x is a large number.
log_exp_for_large_numbers<-function(x)
{
x_array<-c()
while (x>0)
{
ifelse (x>710, x_array<-cbind(x_array,709), x_array<-cbind(x_array,x))
x<-x-709
}
out<-0
for (i in 1:length(x_array))
{
out<-out+log(exp(x_array[i]))
}
return(out)
}
log_exp_for_large_numbers(31568)
[1]
31568
Sincerely,
Olena Ageyeva
------------------------------
From: aallen at hbs.edu
To: gov2001-l at lists.fas.harvard.edu
Date: Wed, 11 Mar 2009 00:07:38 -0400
Subject: [gov2001-l] problems with log(exp().... Sequential processing in R
causing issues
In attempting to debug my function I came across the below anomaly which
I cannot figure out.
Whenever I try to evaluate log(exp(A)) for any A >= 710 R evaluates the
result as ?Inf? which is clearly not the case. See R output below:
log(exp(708))
[1] 708
log(exp(709))
[1] 709
*> log(exp(710))*
*[1] Inf*
*> log(exp(832))*
*[1] Inf*
I assume that this happens because R evaluates the function inside out:
i.e. first tries to take the exp(710) which produces more digits than R will
hold: hence R evaluates exp(710) as Inf, and then applies the log fxn to
this result? log(Inf) = Inf.
Ok logical. However, incredibly inconvenient for purposes of this homework
set, and for optim() in general, as many negative values of B result in
ln(1+e^(-xiB)) > ln(1+e^710) ? creating an ?-Inf? evaluation of the log
likelihood. This sequential calculation jams my function when I try to
optimize since R cannot sum (Inf+510+712) for example ? Error in -sum(formula
in here J) : invalid 'type' (character) of argument.
Has anyone else encountered this issue?
Any suggestions?
Thanks,
Abigail
*From:* gov2001-l-bounces at lists.fas.harvard.edu [mailto:
gov2001-l-bounces at lists.fas.harvard.edu] *On Behalf Of *Patrick Lam
*Sent:* Tuesday, March 10, 2009 2:50 PM
*To:* gov2001-l at lists.fas.harvard.edu
*Subject:* Re: [gov2001-l] Syntax to specify bounds in using L-BFGS-B?
in the optim help file under "Arguments", there are two arguments that you
are looking for.
On Tue, Mar 10, 2009 at 2:26 PM, Allen, Abigail <aallen at hbs.edu> wrote:
What is the syntax to specify the bounds using L-BFGS-B. I?ve looked
through the R-help file but can?t figure out how I tell optim what bounds to
use?.
*From:* gov2001-l-bounces at lists.fas.harvard.edu [mailto:
gov2001-l-bounces at lists.fas.harvard.edu] *On Behalf Of *Lee, Clarence
*Sent:* Monday, March 09, 2009 10:44 PM
*To:* gov2001-l at lists.fas.harvard.edu
*Subject:* Re: [gov2001-l] EXO 1
I agree with Olena in that you need to use L-BFGS-B to get over that. My
question is more of a ?good-standard? practice:
I notice that the problem to the reparametrization trick is that the SE?s
gets messed up when we reparamterize ? which seems to be a bit kludgy since
usually we care about the mean and the SE?s. If reparamterization is a
standard technique that most statistician use, then is there a standard
technique with BFGS that helps us find the SE? (e.g. without specifying the
bounds using L-BFGS-B)
Any insights?
Thanks,
Clarence
*From:* gov2001-l-bounces at lists.fas.harvard.edu [mailto:
gov2001-l-bounces at lists.fas.harvard.edu] *On Behalf Of *Olena Ageyeva
*Sent:* Monday, March 09, 2009 8:46 PM
*To:* gov2001-l at lists.fas.harvard.edu
*Subject:* Re: [gov2001-l] EXO 1
I've got the same problem with re-parametrization and then I found out that
instead re-parametrization in this particular case we are advised to use
L-BFGS-B method in optim(), so we can set a lower and an upper limits for
out point of interest. So, the point of interest will never get out of this
interval, which is sort of reparameterization, right?
The section notes says that log-likelihood (I guess it's in general for any
LL-function?) is invariant to re-parametrization (any?), but I was not much
successful applying this technique for exponential log-likelihood function.
With re-parametrization done exactly the way it's in section note I get
totally different result from what I derived analytically.
As for pnorm(opt.1000 - 1.96*se) I think that it could be done separately
for se and opt.1000 - they both are found for reparameterized function. So,
we basically looking for opt and se in different coordinate system and then
re-parameterized to get them back to our original coordinate system. Am I
right?
Sincerely,
Olena Ageyeva
Date: Mon, 9 Mar 2009 20:00:14 -0400
From: charlotte.cavaille at gmail.com
To: gov2001-l at lists.fas.harvard.edu
Subject: [gov2001-l] EXO 1
hi!
Quick (but fundamental question)
in exercise 1 gamma needs to be positive, so I re-parameterized...In
order to get the right MLE from optim() i didn't
forget to "re-parameterize" again...and i get the same result as the
one i found analytically. However, things change when i look for the
SE. Indeed, I get different answers analytically and with R (using the
hessian). I know this comes from the re-parameterization because I
don't have this issue when I change the whole function and do not
re-parameterize gamma. So my question is:
- if I re-parametterize, how do I apply the transformation to the
hessian to get the right result
how does that fit with the section notes that follows, why do we take
"pnorm" (this is the first transformation that is applied in the
ll.binom function) of "opt.1000 - 1.96*se" and not of "se" for
instance???
#binomial log-likelihood (N = # of trials for each observation)
ll.binom <- function(par, y, N){
# reparameterize pi; only search over [0,1]
p <- pnorm(par)
# log-likelihood
out <- sum(y*log(p) + (N-y)*log(1-p))
return(out)
}
# compare to wald ci
se <- sqrt(solve(-optim(par=2, fn=ll.binom, y=samp.1000, N=10,
method="BFGS",
control=list(fnscale=-1), hessian=T)$hessian))
#$
wald.ci <- c( pnorm(opt.1000 - 1.96*se), pnorm(opt.1000 + 1.96*se))
wald.ci # 0.7364839 0.7535663
- if i do not re-parameterize in order to be done with it and have
both my analytical
and my R result fit, how can i justify i am not re-parameterizing gamma?!
thanks!
charlotte
_______________________________________________
gov2001-l mailing list
gov2001-l at lists.fas.harvard.edu
http://lists.fas.harvard.edu/mailman/listinfo/gov2001-l
------------------------------
Windows Live? Groups: Create an online spot for your favorite groups to
meet. Check it out.
_______________________________________________
gov2001-l mailing list
gov2001-l at lists.fas.harvard.edu
http://lists.fas.harvard.edu/mailman/listinfo/gov2001-l
--
Patrick Lam
Department of Government and Institute for Quantitative Social Science,
Harvard University
http://www.people.fas.harvard.edu/~plam
------------------------------
Windows Live? Contacts: Organize your contact list. Check it
out.<http://windowslive.com/connect/post/marcusatmicrosoft.spaces.live.c…
_______________________________________________
gov2001-l mailing list
gov2001-l at lists.fas.harvard.edu
http://lists.fas.harvard.edu/mailman/listinfo/gov2001-l
6:57 p.m.
New subject: [gov2001-l] problems with log(exp().... Sequential processing in R causing issues
Hi guys,
I am running into the same issues as Abigail, and I'm taking the log likelihood
exactly as is, derived from slide 3 of Gary's lecture titled "Single Equation
Models." Are we not suppose to use this likelihood?
I've checked the inputs for dimensionality and errors as well, and I do not think
that's causing the problem. This seems to be a numerical overflow issue when dealing
with Logit Models. Do you guys know if there is a general workaround for this aside from
manually coding in approximations ourselves?
Clarence
From: gov2001-l-bounces at lists.fas.harvard.edu [mailto:gov2001-l-bounces at
lists.fas.harvard.edu] On Behalf Of Miya Woolfalk
Sent: Wednesday, March 11, 2009 1:06 PM
To: gov2001-l at lists.fas.harvard.edu
Subject: Re: [gov2001-l] problems with log(exp().... Sequential processing in R causing
issues
Abigail,
Before trying to implement a work around, you should make sure that your log likelihood
function is programed correctly. I have not heard of others encountering this issue so my
suspicion is that you have an error in your function.
Best,
Miya
On Wed, Mar 11, 2009 at 1:02 PM, Olena Ageyeva <eageeva at
hotmail.com<mailto:eageeva at hotmail.com>> wrote:
I would try to check that your variable doesn't get grater then 710.
This simple modification makes the algorithm work:
m=605
n=110
log(exp(m)*exp(n))
[1] Inf
log(exp(m))+log(exp(n))
[1] 715
I would just create a function for log(exp(x)), where x is a large number.
log_exp_for_large_numbers<-function(x)
{
x_array<-c()
while (x>0)
{
ifelse (x>710, x_array<-cbind(x_array,709), x_array<-cbind(x_array,x))
x<-x-709
}
out<-0
for (i in 1:length(x_array))
{
out<-out+log(exp(x_array[i]))
}
return(out)
}
log_exp_for_large_numbers(31568)
[1] 31568
Sincerely,
Olena Ageyeva
________________________________
From: aallen at hbs.edu<mailto:aallen at hbs.edu>
To: gov2001-l at lists.fas.harvard.edu<mailto:gov2001-l at lists.fas.harvard.edu>
Date: Wed, 11 Mar 2009 00:07:38 -0400
Subject: [gov2001-l] problems with log(exp().... Sequential processing in R causing
issues
In attempting to debug my function I came across the below anomaly which I cannot figure
out.
Whenever I try to evaluate log(exp(A)) for any A >= 710 R evaluates the result as
"Inf" which is clearly not the case. See R output below:
log(exp(708))
[1] 708
log(exp(709))
[1] 709
log(exp(710))
[1] Inf
log(exp(832))
[1] Inf
I assume that this happens because R evaluates the function inside out: i.e. first tries
to take the exp(710) which produces more digits than R will hold: hence R evaluates
exp(710) as Inf, and then applies the log fxn to this result--> log(Inf) = Inf.
Ok logical. However, incredibly inconvenient for purposes of this homework set, and for
optim() in general, as many negative values of B result in ln(1+e^(-xiB)) > ln(1+e^710)
--> creating an "-Inf" evaluation of the log likelihood. This sequential
calculation jams my function when I try to optimize since R cannot sum (Inf+510+712) for
example --> Error in -sum(formula in here :)) : invalid 'type' (character) of
argument.
Has anyone else encountered this issue?
Any suggestions?
Thanks,
Abigail
From: gov2001-l-bounces at lists.fas.harvard.edu<mailto:gov2001-l-bounces at
lists.fas.harvard.edu> [mailto:gov2001-l-bounces at
lists.fas.harvard.edu<mailto:gov2001-l-bounces at lists.fas.harvard.edu>] On Behalf
Of Patrick Lam
Sent: Tuesday, March 10, 2009 2:50 PM
To: gov2001-l at lists.fas.harvard.edu<mailto:gov2001-l at lists.fas.harvard.edu>
Subject: Re: [gov2001-l] Syntax to specify bounds in using L-BFGS-B?
in the optim help file under "Arguments", there are two arguments that you are
looking for.
On Tue, Mar 10, 2009 at 2:26 PM, Allen, Abigail <aallen at hbs.edu> wrote:
What is the syntax to specify the bounds using L-BFGS-B. I've looked through the
R-help file but can't figure out how I tell optim what bounds to use....
From: gov2001-l-bounces at lists.fas.harvard.edu [mailto:gov2001-l-bounces at
lists.fas.harvard.edu] On Behalf Of Lee, Clarence
Sent: Monday, March 09, 2009 10:44 PM
To: gov2001-l at lists.fas.harvard.edu
Subject: Re: [gov2001-l] EXO 1
I agree with Olena in that you need to use L-BFGS-B to get over that. My question is
more of a "good-standard" practice:
I notice that the problem to the reparametrization trick is that the SE's gets messed
up when we reparamterize - which seems to be a bit kludgy since usually we care about the
mean and the SE's. If reparamterization is a standard technique that most
statistician use, then is there a standard technique with BFGS that helps us find the SE?
(e.g. without specifying the bounds using L-BFGS-B)
Any insights?
Thanks,
Clarence
From: gov2001-l-bounces at lists.fas.harvard.edu [mailto:gov2001-l-bounces at
lists.fas.harvard.edu] On Behalf Of Olena Ageyeva
Sent: Monday, March 09, 2009 8:46 PM
To: gov2001-l at lists.fas.harvard.edu
Subject: Re: [gov2001-l] EXO 1
I've got the same problem with re-parametrization and then I found out that instead
re-parametrization in this particular case we are advised to use L-BFGS-B method in
optim(), so we can set a lower and an upper limits for out point of interest. So, the
point of interest will never get out of this interval, which is sort of
reparameterization, right?
The section notes says that log-likelihood (I guess it's in general for any
LL-function?) is invariant to re-parametrization (any?), but I was not much successful
applying this technique for exponential log-likelihood function. With re-parametrization
done exactly the way it's in section note I get totally different result from what I
derived analytically.
As for pnorm(opt.1000 - 1.96*se) I think that it could be done separately for se and
opt.1000 - they both are found for reparameterized function. So, we basically looking for
opt and se in different coordinate system and then re-parameterized to get them back to
our original coordinate system. Am I right?
Sincerely,
Olena Ageyeva
Date: Mon, 9 Mar 2009 20:00:14 -0400
From: charlotte.cavaille at gmail.com
To: gov2001-l at lists.fas.harvard.edu
Subject: [gov2001-l] EXO 1
hi!
Quick (but fundamental question)
in exercise 1 gamma needs to be positive, so I re-parameterized...In
order to get the right MLE from optim() i didn't
forget to "re-parameterize" again...and i get the same result as the
one i found analytically. However, things change when i look for the
SE. Indeed, I get different answers analytically and with R (using the
hessian). I know this comes from the re-parameterization because I
don't have this issue when I change the whole function and do not
re-parameterize gamma. So my question is:
- if I re-parametterize, how do I apply the transformation to the
hessian to get the right result
how does that fit with the section notes that follows, why do we take
"pnorm" (this is the first transformation that is applied in the
ll.binom function) of "opt.1000 - 1.96*se" and not of "se" for
instance???
#binomial log-likelihood (N = # of trials for each observation)
ll.binom <- function(par, y, N){
# reparameterize pi; only search over [0,1]
p <- pnorm(par)
# log-likelihood
out <- sum(y*log(p) + (N-y)*log(1-p))
return(out)
}
# compare to wald ci
se <- sqrt(solve(-optim(par=2, fn=ll.binom, y=samp.1000, N=10,
method="BFGS",
control=list(fnscale=-1), hessian=T)$hessian)) #$
wald.ci <- c( pnorm(opt.1000 - 1.96*se), pnorm(opt.1000 + 1.96*se))
wald.ci # 0.7364839 0.7535663
- if i do not re-parameterize in order to be done with it and have
both my analytical
and my R result fit, how can i justify i am not re-parameterizing gamma?!
thanks!
charlotte
_______________________________________________
gov2001-l mailing list
gov2001-l at lists.fas.harvard.edu
http://lists.fas.harvard.edu/mailman/listinfo/gov2001-l
________________________________
Windows Live(tm) Groups: Create an online spot for your favorite groups to meet. Check it
out.
_______________________________________________
gov2001-l mailing list
gov2001-l at lists.fas.harvard.edu
http://lists.fas.harvard.edu/mailman/listinfo/gov2001-l
--
Patrick Lam
Department of Government and Institute for Quantitative Social Science, Harvard
University
http://www.people.fas.harvard.edu/~plam
________________________________
Windows Live(tm) Contacts: Organize your contact list. Check it
out.<http://windowslive.com/connect/post/marcusatmicrosoft.spaces.live.c…
_______________________________________________
gov2001-l mailing list
gov2001-l at lists.fas.harvard.edu<mailto:gov2001-l at lists.fas.harvard.edu>
http://lists.fas.harvard.edu/mailman/listinfo/gov2001-l
--
Miya Woolfalk
Ph.D. Student
Harvard University
Government and Social Policy
7:03 p.m.
New subject: [gov2001-l] problems with log(exp().... Sequential processing in R causing issues
I hate to pile on, but count me as another voice that is seeing the same problem. It feels
like there should be an elegant work around, but elegance evades me. Olena offered a
solution that does work but it is disappointing to think R would require that kind of
trickeration.
From: gov2001-l-bounces at lists.fas.harvard.edu [mailto:gov2001-l-bounces at
lists.fas.harvard.edu] On Behalf Of Lee, Clarence
Sent: Wednesday, March 11, 2009 6:57 PM
To: gov2001-l at lists.fas.harvard.edu
Subject: Re: [gov2001-l] problems with log(exp().... Sequential processing in R causing
issues
Hi guys,
I am running into the same issues as Abigail, and I'm taking the log likelihood
exactly as is, derived from slide 3 of Gary's lecture titled "Single Equation
Models." Are we not suppose to use this likelihood?
I've checked the inputs for dimensionality and errors as well, and I do not think
that's causing the problem. This seems to be a numerical overflow issue when dealing
with Logit Models. Do you guys know if there is a general workaround for this aside from
manually coding in approximations ourselves?
Clarence
From: gov2001-l-bounces at lists.fas.harvard.edu [mailto:gov2001-l-bounces at
lists.fas.harvard.edu] On Behalf Of Miya Woolfalk
Sent: Wednesday, March 11, 2009 1:06 PM
To: gov2001-l at lists.fas.harvard.edu
Subject: Re: [gov2001-l] problems with log(exp().... Sequential processing in R causing
issues
Abigail,
Before trying to implement a work around, you should make sure that your log likelihood
function is programed correctly. I have not heard of others encountering this issue so my
suspicion is that you have an error in your function.
Best,
Miya
On Wed, Mar 11, 2009 at 1:02 PM, Olena Ageyeva <eageeva at
hotmail.com<mailto:eageeva at hotmail.com>> wrote:
I would try to check that your variable doesn't get grater then 710.
This simple modification makes the algorithm work:
m=605
n=110
log(exp(m)*exp(n))
[1] Inf
log(exp(m))+log(exp(n))
[1] 715
I would just create a function for log(exp(x)), where x is a large number.
log_exp_for_large_numbers<-function(x)
{
x_array<-c()
while (x>0)
{
ifelse (x>710, x_array<-cbind(x_array,709), x_array<-cbind(x_array,x))
x<-x-709
}
out<-0
for (i in 1:length(x_array))
{
out<-out+log(exp(x_array[i]))
}
return(out)
}
log_exp_for_large_numbers(31568)
[1] 31568
Sincerely,
Olena Ageyeva
________________________________
From: aallen at hbs.edu<mailto:aallen at hbs.edu>
To: gov2001-l at lists.fas.harvard.edu<mailto:gov2001-l at lists.fas.harvard.edu>
Date: Wed, 11 Mar 2009 00:07:38 -0400
Subject: [gov2001-l] problems with log(exp().... Sequential processing in R causing
issues
In attempting to debug my function I came across the below anomaly which I cannot figure
out.
Whenever I try to evaluate log(exp(A)) for any A >= 710 R evaluates the result as
"Inf" which is clearly not the case. See R output below:
log(exp(708))
[1] 708
log(exp(709))
[1] 709
log(exp(710))
[1] Inf
log(exp(832))
[1] Inf
I assume that this happens because R evaluates the function inside out: i.e. first tries
to take the exp(710) which produces more digits than R will hold: hence R evaluates
exp(710) as Inf, and then applies the log fxn to this result--> log(Inf) = Inf.
Ok logical. However, incredibly inconvenient for purposes of this homework set, and for
optim() in general, as many negative values of B result in ln(1+e^(-xiB)) > ln(1+e^710)
--> creating an "-Inf" evaluation of the log likelihood. This sequential
calculation jams my function when I try to optimize since R cannot sum (Inf+510+712) for
example --> Error in -sum(formula in here :)) : invalid 'type' (character) of
argument.
Has anyone else encountered this issue?
Any suggestions?
Thanks,
Abigail
From: gov2001-l-bounces at lists.fas.harvard.edu<mailto:gov2001-l-bounces at
lists.fas.harvard.edu> [mailto:gov2001-l-bounces at
lists.fas.harvard.edu<mailto:gov2001-l-bounces at lists.fas.harvard.edu>] On Behalf
Of Patrick Lam
Sent: Tuesday, March 10, 2009 2:50 PM
To: gov2001-l at lists.fas.harvard.edu<mailto:gov2001-l at lists.fas.harvard.edu>
Subject: Re: [gov2001-l] Syntax to specify bounds in using L-BFGS-B?
in the optim help file under "Arguments", there are two arguments that you are
looking for.
On Tue, Mar 10, 2009 at 2:26 PM, Allen, Abigail <aallen at hbs.edu> wrote:
What is the syntax to specify the bounds using L-BFGS-B. I've looked through the
R-help file but can't figure out how I tell optim what bounds to use....
From: gov2001-l-bounces at lists.fas.harvard.edu [mailto:gov2001-l-bounces at
lists.fas.harvard.edu] On Behalf Of Lee, Clarence
Sent: Monday, March 09, 2009 10:44 PM
To: gov2001-l at lists.fas.harvard.edu
Subject: Re: [gov2001-l] EXO 1
I agree with Olena in that you need to use L-BFGS-B to get over that. My question is
more of a "good-standard" practice:
I notice that the problem to the reparametrization trick is that the SE's gets messed
up when we reparamterize - which seems to be a bit kludgy since usually we care about the
mean and the SE's. If reparamterization is a standard technique that most
statistician use, then is there a standard technique with BFGS that helps us find the SE?
(e.g. without specifying the bounds using L-BFGS-B)
Any insights?
Thanks,
Clarence
From: gov2001-l-bounces at lists.fas.harvard.edu [mailto:gov2001-l-bounces at
lists.fas.harvard.edu] On Behalf Of Olena Ageyeva
Sent: Monday, March 09, 2009 8:46 PM
To: gov2001-l at lists.fas.harvard.edu
Subject: Re: [gov2001-l] EXO 1
I've got the same problem with re-parametrization and then I found out that instead
re-parametrization in this particular case we are advised to use L-BFGS-B method in
optim(), so we can set a lower and an upper limits for out point of interest. So, the
point of interest will never get out of this interval, which is sort of
reparameterization, right?
The section notes says that log-likelihood (I guess it's in general for any
LL-function?) is invariant to re-parametrization (any?), but I was not much successful
applying this technique for exponential log-likelihood function. With re-parametrization
done exactly the way it's in section note I get totally different result from what I
derived analytically.
As for pnorm(opt.1000 - 1.96*se) I think that it could be done separately for se and
opt.1000 - they both are found for reparameterized function. So, we basically looking for
opt and se in different coordinate system and then re-parameterized to get them back to
our original coordinate system. Am I right?
Sincerely,
Olena Ageyeva
Date: Mon, 9 Mar 2009 20:00:14 -0400
From: charlotte.cavaille at gmail.com
To: gov2001-l at lists.fas.harvard.edu
Subject: [gov2001-l] EXO 1
hi!
Quick (but fundamental question)
in exercise 1 gamma needs to be positive, so I re-parameterized...In
order to get the right MLE from optim() i didn't
forget to "re-parameterize" again...and i get the same result as the
one i found analytically. However, things change when i look for the
SE. Indeed, I get different answers analytically and with R (using the
hessian). I know this comes from the re-parameterization because I
don't have this issue when I change the whole function and do not
re-parameterize gamma. So my question is:
- if I re-parametterize, how do I apply the transformation to the
hessian to get the right result
how does that fit with the section notes that follows, why do we take
"pnorm" (this is the first transformation that is applied in the
ll.binom function) of "opt.1000 - 1.96*se" and not of "se" for
instance???
#binomial log-likelihood (N = # of trials for each observation)
ll.binom <- function(par, y, N){
# reparameterize pi; only search over [0,1]
p <- pnorm(par)
# log-likelihood
out <- sum(y*log(p) + (N-y)*log(1-p))
return(out)
}
# compare to wald ci
se <- sqrt(solve(-optim(par=2, fn=ll.binom, y=samp.1000, N=10,
method="BFGS",
control=list(fnscale=-1), hessian=T)$hessian)) #$
wald.ci <- c( pnorm(opt.1000 - 1.96*se), pnorm(opt.1000 + 1.96*se))
wald.ci # 0.7364839 0.7535663
- if i do not re-parameterize in order to be done with it and have
both my analytical
and my R result fit, how can i justify i am not re-parameterizing gamma?!
thanks!
charlotte
_______________________________________________
gov2001-l mailing list
gov2001-l at lists.fas.harvard.edu
http://lists.fas.harvard.edu/mailman/listinfo/gov2001-l
________________________________
Windows Live(tm) Groups: Create an online spot for your favorite groups to meet. Check it
out.
_______________________________________________
gov2001-l mailing list
gov2001-l at lists.fas.harvard.edu
http://lists.fas.harvard.edu/mailman/listinfo/gov2001-l
--
Patrick Lam
Department of Government and Institute for Quantitative Social Science, Harvard
University
http://www.people.fas.harvard.edu/~plam
________________________________
Windows Live(tm) Contacts: Organize your contact list. Check it
out.<http://windowslive.com/connect/post/marcusatmicrosoft.spaces.live.c…
_______________________________________________
gov2001-l mailing list
gov2001-l at lists.fas.harvard.edu<mailto:gov2001-l at lists.fas.harvard.edu>
http://lists.fas.harvard.edu/mailman/listinfo/gov2001-l
--
Miya Woolfalk
Ph.D. Student
Harvard University
Government and Social Policy
7:53 p.m.
New subject: [gov2001-l] problems with log(exp().... Sequential processing in R causing issues
Just a curiosity, which problem did you get this issue with?
Sincerely,
Olena Ageyeva
From: wschmidt at hbs.edu
To: gov2001-l at lists.fas.harvard.edu
Date: Wed, 11 Mar 2009 19:03:52 -0400
Subject: Re: [gov2001-l] problems with log(exp().... Sequential processing in R causing
issues
I hate to pile on, but count me as another voice that is seeing
the same problem. It feels like there should be an elegant work around, but
elegance evades me. Olena offered a solution that does work but it is disappointing
to think R would require that kind of trickeration.
From:
gov2001-l-bounces at lists.fas.harvard.edu
[mailto:gov2001-l-bounces at lists.fas.harvard.edu] On Behalf Of Lee,
Clarence
Sent: Wednesday, March 11, 2009 6:57 PM
To: gov2001-l at lists.fas.harvard.edu
Subject: Re: [gov2001-l] problems with log(exp().... Sequential
processing in R causing issues
Hi guys,
I am running into the same issues as Abigail, and I?m
taking the log likelihood exactly as is, derived from slide 3 of Gary?s
lecture titled ?Single Equation Models.? Are we not suppose
to use this likelihood?
I?ve checked the inputs for dimensionality and errors as
well, and I do not think that?s causing the problem. This seems to
be a numerical overflow issue when dealing with Logit Models. Do
you guys know if there is a general workaround for this aside from manually
coding in approximations ourselves?
Clarence
From:
gov2001-l-bounces at lists.fas.harvard.edu
[mailto:gov2001-l-bounces at lists.fas.harvard.edu] On Behalf Of Miya
Woolfalk
Sent: Wednesday, March 11, 2009 1:06 PM
To: gov2001-l at lists.fas.harvard.edu
Subject: Re: [gov2001-l] problems with log(exp().... Sequential
processing in R causing issues
Abigail,
Before trying to implement a work around, you should make sure that your log
likelihood function is programed correctly. I have not heard of others
encountering this issue so my suspicion is that you have an error in your
function.
Best,
Miya
On Wed, Mar 11, 2009 at 1:02 PM, Olena Ageyeva <eageeva at hotmail.com> wrote:
I would try to check that your variable doesn't get grater
then 710.
This simple modification makes the algorithm work:
m=605
n=110
log(exp(m)*exp(n))
[1] Inf
log(exp(m))+log(exp(n))
[1] 715
I would just create a function for log(exp(x)), where x is a large number.
log_exp_for_large_numbers<-function(x)
{
x_array<-c()
while (x>0)
{
ifelse (x>710, x_array<-cbind(x_array,709), x_array<-cbind(x_array,x))
x<-x-709
}
out<-0
for (i in 1:length(x_array))
{
out<-out+log(exp(x_array[i]))
}
return(out)
}
log_exp_for_large_numbers(31568)
[1] 31568
Sincerely,
Olena Ageyeva
From: aallen at hbs.edu
To: gov2001-l at lists.fas.harvard.edu
Date: Wed, 11 Mar 2009 00:07:38 -0400
Subject: [gov2001-l] problems with log(exp().... Sequential processing in R
causing issues
In attempting to debug my function
I came across the below anomaly which I cannot figure out.
Whenever I try to evaluate
log(exp(A)) for any A >= 710 R evaluates the result as
?Inf? which is clearly not the case. See R output below:
log(exp(708))
[1] 708
log(exp(709))
[1] 709
log(exp(710))
[1] Inf
log(exp(832))
[1] Inf
I assume that this happens
because R evaluates the function inside out: i.e. first tries to take the
exp(710) which produces more digits than R will hold: hence R evaluates
exp(710) as Inf, and then applies the log fxn to this result? log(Inf) = Inf.
Ok logical. However, incredibly
inconvenient for purposes of this homework set, and for optim() in general, as
many negative values of B result in ln(1+e^(-xiB)) > ln(1+e^710) ? creating an ?-Inf?
evaluation of the log likelihood. This sequential calculation jams my
function when I try to optimize since R cannot sum (Inf+510+712) for example ? Error in
-sum(formula in here J) : invalid 'type' (character) of
argument.
Has anyone else encountered
this issue?
Any suggestions?
Thanks,
Abigail
From: gov2001-l-bounces at lists.fas.harvard.edu [mailto:gov2001-l-bounces at
lists.fas.harvard.edu]
On Behalf Of Patrick Lam
Sent: Tuesday, March 10, 2009 2:50 PM
To: gov2001-l at lists.fas.harvard.edu
Subject: Re: [gov2001-l] Syntax to specify bounds in using L-BFGS-B?
in the optim help file under
"Arguments", there are two arguments that you are looking for.
On Tue, Mar 10, 2009 at 2:26 PM, Allen, Abigail <aallen at hbs.edu>
wrote:
What is the
syntax to specify the bounds using L-BFGS-B. I?ve looked through
the R-help file but can?t figure out how I tell optim what bounds to
use?.
From: gov2001-l-bounces at lists.fas.harvard.edu
[mailto:gov2001-l-bounces at lists.fas.harvard.edu] On Behalf Of Lee,
Clarence
Sent: Monday, March 09, 2009 10:44 PM
To: gov2001-l at lists.fas.harvard.edu
Subject: Re: [gov2001-l] EXO 1
I agree with Olena in that you
need to use L-BFGS-B to get over that. My question is more of a
?good-standard? practice:
I notice that the problem to the
reparametrization trick is that the SE?s gets messed up when we
reparamterize ? which seems to be a bit kludgy since usually we care
about the mean and the SE?s. If reparamterization is a
standard technique that most statistician use, then is there a standard
technique with BFGS that helps us find the SE? (e.g. without specifying the
bounds using L-BFGS-B)
Any insights?
Thanks,
Clarence
From: gov2001-l-bounces at lists.fas.harvard.edu
[mailto:gov2001-l-bounces at lists.fas.harvard.edu] On Behalf Of Olena
Ageyeva
Sent: Monday, March 09, 2009 8:46 PM
To: gov2001-l at lists.fas.harvard.edu
Subject: Re: [gov2001-l] EXO 1
I've got the
same problem with re-parametrization and then I found out that instead
re-parametrization in this particular case we are advised to use L-BFGS-B
method in optim(), so we can set a lower and an upper limits for out point of
interest. So, the point of interest will never get out of this interval, which
is sort of reparameterization, right?
The section notes says that log-likelihood (I guess it's in general for any
LL-function?) is invariant to re-parametrization (any?), but I was not much
successful applying this technique for exponential log-likelihood function.
With re-parametrization done exactly the way it's in section note I get totally
different result from what I derived analytically.
As for pnorm(opt.1000 - 1.96*se) I think that it could be done separately for
se and opt.1000 - they both are found for reparameterized function. So, we
basically looking for opt and se in different coordinate system and then
re-parameterized to get them back to our original coordinate system. Am I
right?
Sincerely,
Olena Ageyeva
Date: Mon, 9 Mar 2009 20:00:14 -0400
From: charlotte.cavaille at gmail.com
To: gov2001-l at lists.fas.harvard.edu
Subject: [gov2001-l] EXO 1
hi!
Quick (but fundamental question)
in exercise 1 gamma needs to be positive, so I
re-parameterized...In
order to get the right MLE from optim() i didn't
forget to "re-parameterize" again...and i
get the same result as
the
one i found analytically. However, things change when
i look for the
SE. Indeed, I get different answers analytically and
with R (using the
hessian). I know this comes from the
re-parameterization because I
don't have this issue when I change the whole
function and do not
re-parameterize gamma. So my question is:
- if I re-parametterize, how do I apply the
transformation to the
hessian to get the right result
how does that fit with the section notes that follows,
why do we take
"pnorm" (this is the first transformation
that is applied in the
ll.binom function) of "opt.1000 - 1.96*se"
and not of
"se" for
instance???
#binomial log-likelihood (N = # of trials for each
observation)
ll.binom <- function(par, y, N){
# reparameterize pi; only search over [0,1]
p <- pnorm(par)
# log-likelihood
out <- sum(y*log(p) + (N-y)*log(1-p))
return(out)
}
# compare to wald ci
se <- sqrt(solve(-optim(par=2, fn=ll.binom,
y=samp.1000, N=10,
method="BFGS",
control=list(fnscale=-1), hessian=T)$hessian)) #$
wald.ci <- c( pnorm(opt.1000 - 1.96*se),
pnorm(opt.1000 + 1.96*se))
wald.ci # 0.7364839 0.7535663
- if i do not re-parameterize in order to be done with
it and have
both my analytical
and my R result fit, how can i justify i am not
re-parameterizing gamma?!
thanks!
charlotte
_______________________________________________
gov2001-l mailing list
gov2001-l at lists.fas.harvard.edu
Windows Live? Groups:
Create an online spot for your favorite groups to meet. Check it out.
_______________________________________________
gov2001-l mailing list
gov2001-l at lists.fas.harvard.edu
http://lists.fas.harvard.edu/mailman/listinfo/gov2001-l
--
Patrick Lam
Department of Government and Institute for Quantitative Social Science, Harvard
University
http://www.people.fas.harvard.edu/~plam
Windows Live? Contacts: Organize your contact list. Check it out.
_______________________________________________
gov2001-l mailing list
gov2001-l at lists.fas.harvard.edu
http://lists.fas.harvard.edu/mailman/listinfo/gov2001-l
--
Miya Woolfalk
Ph.D. Student
Harvard University
Government and Social Policy
_________________________________________________________________
Windows Live?: Life without walls.
http://windowslive.com/explore?ocid=TXT_TAGLM_WL_allup_1a_explore_032009
9:11 p.m.
New subject: [gov2001-l] problems with log(exp().... Sequential processing in R causing issues
Clarence figured out the fix for his code and sent me and Abigail a note since we seemed
to be the only ones affected. I made the same error as him. Miya was correct, it was a
problem with my function. That said, R does hit a wall when taking an exp() of any value
greater than 710. This could be problematic with other data and I am not sure there is a
clean solution when using the log-likelihood of a logit with a large argument for the
exp(). Thanks for your help.
From: gov2001-l-bounces at lists.fas.harvard.edu [mailto:gov2001-l-bounces at
lists.fas.harvard.edu] On Behalf Of Olena Ageyeva
Sent: Wednesday, March 11, 2009 7:53 PM
To: gov2001-l at lists.fas.harvard.edu
Subject: Re: [gov2001-l] problems with log(exp().... Sequential processing in R causing
issues
Just a curiosity, which problem did you get this issue with?
Sincerely,
Olena Ageyeva
________________________________
From: wschmidt at hbs.edu
To: gov2001-l at lists.fas.harvard.edu
Date: Wed, 11 Mar 2009 19:03:52 -0400
Subject: Re: [gov2001-l] problems with log(exp().... Sequential processing in R causing
issues
I hate to pile on, but count me as another voice that is seeing the same problem. It feels
like there should be an elegant work around, but elegance evades me. Olena offered a
solution that does work but it is disappointing to think R would require that kind of
trickeration.
From: gov2001-l-bounces at lists.fas.harvard.edu [mailto:gov2001-l-bounces at
lists.fas.harvard.edu] On Behalf Of Lee, Clarence
Sent: Wednesday, March 11, 2009 6:57 PM
To: gov2001-l at lists.fas.harvard.edu
Subject: Re: [gov2001-l] problems with log(exp().... Sequential processing in R causing
issues
Hi guys,
I am running into the same issues as Abigail, and I'm taking the log likelihood
exactly as is, derived from slide 3 of Gary's lecture titled "Single Equation
Models." Are we not suppose to use this likelihood?
I've checked the inputs for dimensionality and errors as well, and I do not think
that's causing the problem. This seems to be a numerical overflow issue when dealing
with Logit Models. Do you guys know if there is a general workaround for this aside from
manually coding in approximations ourselves?
Clarence
From: gov2001-l-bounces at lists.fas.harvard.edu [mailto:gov2001-l-bounces at
lists.fas.harvard.edu] On Behalf Of Miya Woolfalk
Sent: Wednesday, March 11, 2009 1:06 PM
To: gov2001-l at lists.fas.harvard.edu
Subject: Re: [gov2001-l] problems with log(exp().... Sequential processing in R causing
issues
Abigail,
Before trying to implement a work around, you should make sure that your log likelihood
function is programed correctly. I have not heard of others encountering this issue so my
suspicion is that you have an error in your function.
Best,
Miya
On Wed, Mar 11, 2009 at 1:02 PM, Olena Ageyeva <eageeva at hotmail.com> wrote:
I would try to check that your variable doesn't get grater then 710.
This simple modification makes the algorithm work:
m=605
n=110
log(exp(m)*exp(n))
[1] Inf
log(exp(m))+log(exp(n))
[1] 715
I would just create a function for log(exp(x)), where x is a large number.
log_exp_for_large_numbers<-function(x)
{
x_array<-c()
while (x>0)
{
ifelse (x>710, x_array<-cbind(x_array,709), x_array<-cbind(x_array,x))
x<-x-709
}
out<-0
for (i in 1:length(x_array))
{
out<-out+log(exp(x_array[i]))
}
return(out)
}
log_exp_for_large_numbers(31568)
[1] 31568
Sincerely,
Olena Ageyeva
________________________________
From: aallen at hbs.edu
To: gov2001-l at lists.fas.harvard.edu
Date: Wed, 11 Mar 2009 00:07:38 -0400
Subject: [gov2001-l] problems with log(exp().... Sequential processing in R causing
issues
In attempting to debug my function I came across the below anomaly which I cannot figure
out.
Whenever I try to evaluate log(exp(A)) for any A >= 710 R evaluates the result as
"Inf" which is clearly not the case. See R output below:
log(exp(708))
[1] 708
log(exp(709))
[1] 709
log(exp(710))
[1] Inf
log(exp(832))
[1] Inf
I assume that this happens because R evaluates the function inside out: i.e. first tries
to take the exp(710) which produces more digits than R will hold: hence R evaluates
exp(710) as Inf, and then applies the log fxn to this result--> log(Inf) = Inf.
Ok logical. However, incredibly inconvenient for purposes of this homework set, and for
optim() in general, as many negative values of B result in ln(1+e^(-xiB)) > ln(1+e^710)
--> creating an "-Inf" evaluation of the log likelihood. This sequential
calculation jams my function when I try to optimize since R cannot sum (Inf+510+712) for
example --> Error in -sum(formula in here :)) : invalid 'type' (character) of
argument.
Has anyone else encountered this issue?
Any suggestions?
Thanks,
Abigail
From: gov2001-l-bounces at lists.fas.harvard.edu [mailto:gov2001-l-bounces at
lists.fas.harvard.edu] On Behalf Of Patrick Lam
Sent: Tuesday, March 10, 2009 2:50 PM
To: gov2001-l at lists.fas.harvard.edu
Subject: Re: [gov2001-l] Syntax to specify bounds in using L-BFGS-B?
in the optim help file under "Arguments", there are two arguments that you are
looking for.
On Tue, Mar 10, 2009 at 2:26 PM, Allen, Abigail <aallen at hbs.edu> wrote:
What is the syntax to specify the bounds using L-BFGS-B. I've looked through the
R-help file but can't figure out how I tell optim what bounds to use....
From: gov2001-l-bounces at lists.fas.harvard.edu [mailto:gov2001-l-bounces at
lists.fas.harvard.edu] On Behalf Of Lee, Clarence
Sent: Monday, March 09, 2009 10:44 PM
To: gov2001-l at lists.fas.harvard.edu
Subject: Re: [gov2001-l] EXO 1
I agree with Olena in that you need to use L-BFGS-B to get over that. My question is
more of a "good-standard" practice:
I notice that the problem to the reparametrization trick is that the SE's gets messed
up when we reparamterize - which seems to be a bit kludgy since usually we care about the
mean and the SE's. If reparamterization is a standard technique that most
statistician use, then is there a standard technique with BFGS that helps us find the SE?
(e.g. without specifying the bounds using L-BFGS-B)
Any insights?
Thanks,
Clarence
From: gov2001-l-bounces at lists.fas.harvard.edu [mailto:gov2001-l-bounces at
lists.fas.harvard.edu] On Behalf Of Olena Ageyeva
Sent: Monday, March 09, 2009 8:46 PM
To: gov2001-l at lists.fas.harvard.edu
Subject: Re: [gov2001-l] EXO 1
I've got the same problem with re-parametrization and then I found out that instead
re-parametrization in this particular case we are advised to use L-BFGS-B method in
optim(), so we can set a lower and an upper limits for out point of interest. So, the
point of interest will never get out of this interval, which is sort of
reparameterization, right?
The section notes says that log-likelihood (I guess it's in general for any
LL-function?) is invariant to re-parametrization (any?), but I was not much successful
applying this technique for exponential log-likelihood function. With re-parametrization
done exactly the way it's in section note I get totally different result from what I
derived analytically.
As for pnorm(opt.1000 - 1.96*se) I think that it could be done separately for se and
opt.1000 - they both are found for reparameterized function. So, we basically looking for
opt and se in different coordinate system and then re-parameterized to get them back to
our original coordinate system. Am I right?
Sincerely,
Olena Ageyeva
Date: Mon, 9 Mar 2009 20:00:14 -0400
From: charlotte.cavaille at gmail.com
To: gov2001-l at lists.fas.harvard.edu
Subject: [gov2001-l] EXO 1
hi!
Quick (but fundamental question)
in exercise 1 gamma needs to be positive, so I re-parameterized...In
order to get the right MLE from optim() i didn't
forget to "re-parameterize" again...and i get the same result as the
one i found analytically. However, things change when i look for the
SE. Indeed, I get different answers analytically and with R (using the
hessian). I know this comes from the re-parameterization because I
don't have this issue when I change the whole function and do not
re-parameterize gamma. So my question is:
- if I re-parametterize, how do I apply the transformation to the
hessian to get the right result
how does that fit with the section notes that follows, why do we take
"pnorm" (this is the first transformation that is applied in the
ll.binom function) of "opt.1000 - 1.96*se" and not of "se" for
instance???
#binomial log-likelihood (N = # of trials for each observation)
ll.binom <- function(par, y, N){
# reparameterize pi; only search over [0,1]
p <- pnorm(par)
# log-likelihood
out <- sum(y*log(p) + (N-y)*log(1-p))
return(out)
}
# compare to wald ci
se <- sqrt(solve(-optim(par=2, fn=ll.binom, y=samp.1000, N=10,
method="BFGS",
control=list(fnscale=-1), hessian=T)$hessian)) #$
wald.ci <- c( pnorm(opt.1000 - 1.96*se), pnorm(opt.1000 + 1.96*se))
wald.ci # 0.7364839 0.7535663
- if i do not re-parameterize in order to be done with it and have
both my analytical
and my R result fit, how can i justify i am not re-parameterizing gamma?!
thanks!
charlotte
_______________________________________________
gov2001-l mailing list
gov2001-l at lists.fas.harvard.edu
http://lists.fas.harvard.edu/mailman/listinfo/gov2001-l
________________________________
Windows Live(tm) Groups: Create an online spot for your favorite groups to meet. Check it
out.
_______________________________________________
gov2001-l mailing list
gov2001-l at lists.fas.harvard.edu
http://lists.fas.harvard.edu/mailman/listinfo/gov2001-l
--
Patrick Lam
Department of Government and Institute for Quantitative Social Science, Harvard
University
http://www.people.fas.harvard.edu/~plam
________________________________
Windows Live(tm) Contacts: Organize your contact list. Check it out.
_______________________________________________
gov2001-l mailing list
gov2001-l at lists.fas.harvard.edu
http://lists.fas.harvard.edu/mailman/listinfo/gov2001-l
--
Miya Woolfalk
Ph.D. Student
Harvard University
Government and Social Policy
________________________________
Windows Live(tm): Keep your life in sync. Check it
out.<http://windowslive.com/explore?ocid=TXT_TAGLM_WL_allup_1b_explore_0…
9:26 p.m.
New subject: [gov2001-l] problems with log(exp().... Sequential processing in R causing issues
Glad to hear it was resolved. If something like this happens, one quick fix
is to rescale some of your X's. Also, this is somewhat related to why we
like to work with the log scale (as in log-likelihood). It's really easy to
hit numerical limits when we multiply stuff rather than adding stuff.
On Wed, Mar 11, 2009 at 9:11 PM, Schmidt, William <wschmidt at hbs.edu> wrote:
Clarence figured out the fix for his code and sent me
and Abigail a note
since we seemed to be the only ones affected. I made the same error as him.
Miya was correct, it was a problem with my function. That said, R does hit a
wall when taking an exp() of any value greater than 710. This could be
problematic with other data and I am not sure there is a clean solution when
using the log-likelihood of a logit with a large argument for the exp().
Thanks for your help.
*From:* gov2001-l-bounces at lists.fas.harvard.edu [mailto:
gov2001-l-bounces at lists.fas.harvard.edu] *On Behalf Of *Olena Ageyeva
*Sent:* Wednesday, March 11, 2009 7:53 PM
*To:* gov2001-l at lists.fas.harvard.edu
*Subject:* Re: [gov2001-l] problems with log(exp().... Sequential
processing in R causing issues
Just a curiosity, which problem did you get this issue with?
Sincerely,
Olena Ageyeva
------------------------------
From: wschmidt at hbs.edu
To: gov2001-l at lists.fas.harvard.edu
Date: Wed, 11 Mar 2009 19:03:52 -0400
Subject: Re: [gov2001-l] problems with log(exp().... Sequential processing
in R causing issues
I hate to pile on, but count me as another voice that is seeing the same
problem. It feels like there should be an elegant work around, but elegance
evades me. Olena offered a solution that does work but it is disappointing
to think R would require that kind of trickeration.
*From:* gov2001-l-bounces at lists.fas.harvard.edu [mailto:
gov2001-l-bounces at lists.fas.harvard.edu] *On Behalf Of *Lee, Clarence
*Sent:* Wednesday, March 11, 2009 6:57 PM
*To:* gov2001-l at lists.fas.harvard.edu
*Subject:* Re: [gov2001-l] problems with log(exp().... Sequential
processing in R causing issues
Hi guys,
I am running into the same issues as Abigail, and I?m taking the log
likelihood exactly as is, derived from slide 3 of Gary?s lecture titled
?Single Equation Models.? Are we not suppose to use this likelihood?
I?ve checked the inputs for dimensionality and errors as well, and I do not
think that?s causing the problem. This seems to be a numerical overflow
issue when dealing with Logit Models. Do you guys know if there is a
general workaround for this aside from manually coding in approximations
ourselves?
Clarence
*From:* gov2001-l-bounces at lists.fas.harvard.edu [mailto:
gov2001-l-bounces at lists.fas.harvard.edu] *On Behalf Of *Miya Woolfalk
*Sent:* Wednesday, March 11, 2009 1:06 PM
*To:* gov2001-l at lists.fas.harvard.edu
*Subject:* Re: [gov2001-l] problems with log(exp().... Sequential
processing in R causing issues
Abigail,
Before trying to implement a work around, you should make sure that your
log likelihood function is programed correctly. I have not heard of others
encountering this issue so my suspicion is that you have an error in your
function.
Best,
Miya
On Wed, Mar 11, 2009 at 1:02 PM, Olena Ageyeva <eageeva at hotmail.com>
wrote:
I would try to check that your variable doesn't get grater then 710.
This simple modification makes the algorithm work:
--
Patrick Lam
Department of Government and Institute for Quantitative Social Science,
Harvard University
http://www.people.fas.harvard.edu/~plam
m=605
n=110
log(exp(m)*exp(n))
[1] Inf
log(exp(m))+log(exp(n))
[1] 715
I would just create a function for log(exp(x)), where x is a large number.
log_exp_for_large_numbers<-function(x)
{
x_array<-c()
while (x>0)
{
ifelse (x>710, x_array<-cbind(x_array,709), x_array<-cbind(x_array,x))
x<-x-709
}
out<-0
for (i in 1:length(x_array))
{
out<-out+log(exp(x_array[i]))
}
return(out)
}
log_exp_for_large_numbers(31568)
[1]
31568
Sincerely,
Olena Ageyeva
------------------------------
From: aallen at hbs.edu
To: gov2001-l at lists.fas.harvard.edu
Date: Wed, 11 Mar 2009 00:07:38 -0400
Subject: [gov2001-l] problems with log(exp().... Sequential processing in R
causing issues
In attempting to debug my function I came across the below anomaly which I
cannot figure out.
Whenever I try to evaluate log(exp(A)) for any A >= 710 R evaluates the
result as ?Inf? which is clearly not the case. See R output below:
log(exp(708))
[1] 708
log(exp(709))
[1] 709
*> log(exp(710))*
*[1] Inf*
*> log(exp(832))*
*[1] Inf*
I assume that this happens because R evaluates the function inside out:
i.e. first tries to take the exp(710) which produces more digits than R will
hold: hence R evaluates exp(710) as Inf, and then applies the log fxn to
this result? log(Inf) = Inf.
Ok logical. However, incredibly inconvenient for purposes of this homework
set, and for optim() in general, as many negative values of B result in
ln(1+e^(-xiB)) > ln(1+e^710) ? creating an ?-Inf? evaluation of the log
likelihood. This sequential calculation jams my function when I try to
optimize since R cannot sum (Inf+510+712) for example ? Error in -sum(formula
in here J) : invalid 'type' (character) of argument.
Has anyone else encountered this issue?
Any suggestions?
Thanks,
Abigail
*From:* gov2001-l-bounces at lists.fas.harvard.edu [mailto:
gov2001-l-bounces at lists.fas.harvard.edu] *On Behalf Of *Patrick Lam
*Sent:* Tuesday, March 10, 2009 2:50 PM
*To:* gov2001-l at lists.fas.harvard.edu
*Subject:* Re: [gov2001-l] Syntax to specify bounds in using L-BFGS-B?
in the optim help file under "Arguments", there are two arguments that you
are looking for.
On Tue, Mar 10, 2009 at 2:26 PM, Allen, Abigail <aallen at hbs.edu> wrote:
What is the syntax to specify the bounds using L-BFGS-B. I?ve looked
through the R-help file but can?t figure out how I tell optim what bounds to
use?.
*From:* gov2001-l-bounces at lists.fas.harvard.edu [mailto:
gov2001-l-bounces at lists.fas.harvard.edu] *On Behalf Of *Lee, Clarence
*Sent:* Monday, March 09, 2009 10:44 PM
*To:* gov2001-l at lists.fas.harvard.edu
*Subject:* Re: [gov2001-l] EXO 1
I agree with Olena in that you need to use L-BFGS-B to get over that. My
question is more of a ?good-standard? practice:
I notice that the problem to the reparametrization trick is that the SE?s
gets messed up when we reparamterize ? which seems to be a bit kludgy since
usually we care about the mean and the SE?s. If reparamterization is a
standard technique that most statistician use, then is there a standard
technique with BFGS that helps us find the SE? (e.g. without specifying the
bounds using L-BFGS-B)
Any insights?
Thanks,
Clarence
*From:* gov2001-l-bounces at lists.fas.harvard.edu [mailto:
gov2001-l-bounces at lists.fas.harvard.edu] *On Behalf Of *Olena Ageyeva
*Sent:* Monday, March 09, 2009 8:46 PM
*To:* gov2001-l at lists.fas.harvard.edu
*Subject:* Re: [gov2001-l] EXO 1
I've got the same problem with re-parametrization and then I found out that
instead re-parametrization in this particular case we are advised to use
L-BFGS-B method in optim(), so we can set a lower and an upper limits for
out point of interest. So, the point of interest will never get out of this
interval, which is sort of reparameterization, right?
The section notes says that log-likelihood (I guess it's in general for any
LL-function?) is invariant to re-parametrization (any?), but I was not much
successful applying this technique for exponential log-likelihood function.
With re-parametrization done exactly the way it's in section note I get
totally different result from what I derived analytically.
As for pnorm(opt.1000 - 1.96*se) I think that it could be done separately
for se and opt.1000 - they both are found for reparameterized function. So,
we basically looking for opt and se in different coordinate system and then
re-parameterized to get them back to our original coordinate system. Am I
right?
Sincerely,
Olena Ageyeva
Date: Mon, 9 Mar 2009 20:00:14 -0400
From: charlotte.cavaille at gmail.com
To: gov2001-l at lists.fas.harvard.edu
Subject: [gov2001-l] EXO 1
hi!
Quick (but fundamental question)
in exercise 1 gamma needs to be positive, so I re-parameterized...In
order to get the right MLE from optim() i didn't
forget to "re-parameterize" again...and i get the same result as the
one i found analytically. However, things change when i look for the
SE. Indeed, I get different answers analytically and with R (using the
hessian). I know this comes from the re-parameterization because I
don't have this issue when I change the whole function and do not
re-parameterize gamma. So my question is:
- if I re-parametterize, how do I apply the transformation to the
hessian to get the right result
how does that fit with the section notes that follows, why do we take
"pnorm" (this is the first transformation that is applied in the
ll.binom function) of "opt.1000 - 1.96*se" and not of "se" for
instance???
#binomial log-likelihood (N = # of trials for each observation)
ll.binom <- function(par, y, N){
# reparameterize pi; only search over [0,1]
p <- pnorm(par)
# log-likelihood
out <- sum(y*log(p) + (N-y)*log(1-p))
return(out)
}
# compare to wald ci
se <- sqrt(solve(-optim(par=2, fn=ll.binom, y=samp.1000, N=10,
method="BFGS",
control=list(fnscale=-1), hessian=T)$hessian))
#$
wald.ci <- c( pnorm(opt.1000 - 1.96*se), pnorm(opt.1000 + 1.96*se))
wald.ci # 0.7364839 0.7535663
- if i do not re-parameterize in order to be done with it and have
both my analytical
and my R result fit, how can i justify i am not re-parameterizing gamma?!
thanks!
charlotte
_______________________________________________
gov2001-l mailing list
gov2001-l at lists.fas.harvard.edu
http://lists.fas.harvard.edu/mailman/listinfo/gov2001-l
------------------------------
Windows Live? Groups: Create an online spot for your favorite groups to
meet. Check it out.
_______________________________________________
gov2001-l mailing list
gov2001-l at lists.fas.harvard.edu
http://lists.fas.harvard.edu/mailman/listinfo/gov2001-l
--
Patrick Lam
Department of Government and Institute for Quantitative Social Science,
Harvard University
http://www.people.fas.harvard.edu/~plam
------------------------------
Windows Live? Contacts: Organize your contact list. Check it out.
_______________________________________________
gov2001-l mailing list
gov2001-l at lists.fas.harvard.edu
http://lists.fas.harvard.edu/mailman/listinfo/gov2001-l
--
Miya Woolfalk
Ph.D. Student
Harvard University
Government and Social Policy
------------------------------
Windows Live?: Keep your life in sync. Check it
out.<http://windowslive.com/explore?ocid=TXT_TAGLM_WL_allup_1b_explore_0…
_______________________________________________
gov2001-l mailing list
gov2001-l at lists.fas.harvard.edu
http://lists.fas.harvard.edu/mailman/listinfo/gov2001-l
12 Mar
12 Mar
12:33 a.m.
New subject: [gov2001-l] Using Zelig to check our results in 2.3-2.5
Does anyone know the commands in Zelig to quickly check our results for 2.3-2.5?
Thanks,
Clarence
2:14 a.m.
New subject: [gov2001-l] Using Zelig to check our results in 2.3-2.5
for 2.3, for example, something like
data <- as.data.frame(cbind(y,X)) %which I defined earlier to be the vector
and matrix of interest
names(data) <- c("incident", "intercept", "temperature",
"pressure")
output <- zelig(incident ~ temperature +
pressure, data = data, model = "logit")
Should confirm your results.
On Thu, Mar 12, 2009 at 12:33 AM, Lee, Clarence <clee at hbs.edu> wrote:
Does anyone know the commands in Zelig to quickly
check our results for
2.3-2.5?
Thanks,
Clarence
_______________________________________________
gov2001-l mailing list
gov2001-l at lists.fas.harvard.edu
http://lists.fas.harvard.edu/mailman/listinfo/gov2001-l
2:23 a.m.
New subject: [gov2001-l] Using Zelig to check our results in 2.3-2.5
Hi,
Maybe I?m missing something in my understanding, but how do you evaluated the difference
in the expected probability in the data when checking for Zelig? Do you mind describing
how you specified the inputs y and X?
Thanks!
Clarence
From: gov2001-l-bounces at lists.fas.harvard.edu [mailto:gov2001-l-bounces at
lists.fas.harvard.edu] On Behalf Of Z. A. Townsend
Sent: Thursday, March 12, 2009 2:14 AM
To: gov2001-l at lists.fas.harvard.edu
Subject: Re: [gov2001-l] Using Zelig to check our results in 2.3-2.5
for 2.3, for example, something like
data <- as.data.frame(cbind(y,X)) %which I defined earlier to be the vector and matrix
of interest
names(data) <- c("incident", "intercept", "temperature",
"pressure")
output <- zelig(incident ~ temperature + pressure, data = data, model =
"logit")
Should confirm your results.
On Thu, Mar 12, 2009 at 12:33 AM, Lee, Clarence <clee at hbs.edu<mailto:clee at
hbs.edu>> wrote:
Does anyone know the commands in Zelig to quickly check our results for 2.3-2.5?
Thanks,
Clarence
_______________________________________________
gov2001-l mailing list
gov2001-l at lists.fas.harvard.edu<mailto:gov2001-l at lists.fas.harvard.edu>
http://lists.fas.harvard.edu/mailman/listinfo/gov2001-l
2:25 a.m.
New subject: [gov2001-l] Using Zelig to check our results in 2.3-2.5
it might help to look at the Zelig manual, particularly the setx and sim
functions:
http://gking.harvard.edu/zelig/docs/index.html
On Thu, Mar 12, 2009 at 2:23 AM, Lee, Clarence <clee at hbs.edu> wrote:
Hi,
Maybe I?m missing something in my understanding, but how do you evaluated
the difference in the expected probability in the data when checking for
Zelig? Do you mind describing how you specified the inputs y and X?
Thanks!
Clarence
*From:* gov2001-l-bounces at lists.fas.harvard.edu [mailto:
gov2001-l-bounces at lists.fas.harvard.edu] *On Behalf Of *Z. A. Townsend
*Sent:* Thursday, March 12, 2009 2:14 AM
*To:* gov2001-l at lists.fas.harvard.edu
*Subject:* Re: [gov2001-l] Using Zelig to check our results in 2.3-2.5
for 2.3, for example, something like
data <- as.data.frame(cbind(y,X)) %which I defined earlier to be the vector
and matrix of interest
names(data) <- c("incident", "intercept", "temperature",
"pressure")
output <- zelig(incident ~ temperature +
pressure, data = data, model = "logit")
Should confirm your results.
On Thu, Mar 12, 2009 at 12:33 AM, Lee, Clarence <clee at hbs.edu> wrote:
Does anyone know the commands in Zelig to quickly check our results for
2.3-2.5?
Thanks,
Clarence
_______________________________________________
gov2001-l mailing list
gov2001-l at lists.fas.harvard.edu
http://lists.fas.harvard.edu/mailman/listinfo/gov2001-l
_______________________________________________
gov2001-l mailing list
gov2001-l at lists.fas.harvard.edu
http://lists.fas.harvard.edu/mailman/listinfo/gov2001-l
--
Patrick Lam
Department of Government and Institute for Quantitative Social Science,
Harvard University
http://www.people.fas.harvard.edu/~plam
11:48 a.m.
New subject: [gov2001-l] Using Zelig to check our results in 2.3-2.5
I tried this but for some reason got a "NULL" result.
Quoting Patrick Lam <plam at fas.harvard.edu>:
> it might help to look at the Zelig manual, particularly the setx and sim
> functions:
> http://gking.harvard.edu/zelig/docs/index.html
> On Thu, Mar 12, 2009 at 2:23 AM, Lee,
Clarence <clee at hbs.edu> wrote:
> > Hi,
> >
> > > Maybe I?m
missing something in my understanding, but how do you evaluated
> > the difference in the expected probability in the data when checking for
> > Zelig? Do you mind describing how you specified the inputs y and X?
> >
> > > Thanks!
> >
> > Clarence
> >
> > > *From:*
gov2001-l-bounces at lists.fas.harvard.edu [mailto:
> > gov2001-l-bounces at lists.fas.harvard.edu] *On Behalf Of *Z. A. Townsend
> > *Sent:* Thursday, March 12, 2009 2:14 AM
> > *To:* gov2001-l at lists.fas.harvard.edu
> > *Subject:* Re: [gov2001-l] Using Zelig to check our results in 2.3-2.5
> >
> > > for 2.3,
for example, something like
> > > data <- as.data.frame(cbind(y,X))
%which I defined earlier to be the vector
> > and matrix of interest
> > names(data) <- c("incident", "intercept",
"temperature", "pressure")
> > output <- zelig(incident ~ temperature +
> > pressure, data = data, model = "logit")
> >
> > > Should
confirm your results.
> >
> > > On Thu, Mar
12, 2009 at 12:33 AM, Lee, Clarence <clee at hbs.edu> wrote:
> > > Does anyone know the commands in Zelig
to quickly check our results for
> > 2.3-2.5?
> >
> > > Thanks,
> >
> > > Clarence
> >
> > _______________________________________________
> > gov2001-l mailing list
> > gov2001-l at lists.fas.harvard.edu
> > http://lists.fas.harvard.edu/mailman/listinfo/gov2001-l
> >
> > >
_______________________________________________
> > gov2001-l mailing list
> > gov2001-l at lists.fas.harvard.edu
> > http://lists.fas.harvard.edu/mailman/listinfo/gov2001-l
> >
> --
> Patrick Lam
> Department of Government and Institute for Quantitative Social Science,
> Harvard University
> http://www.people.fas.harvard.edu/~plam
check <- zelig(Incident ~ Temperature + Pressure,
model="logit",
+ data=data.frame)
How to cite this model in Zelig:
Kosuke Imai, Gary King, and Oliva Lau. 2008. "logit: Logistic Regression for
Dichotomous Dependent Variables" in Kosuke Imai, Gary King, and Olivia Lau,
"Zelig: Everyone's Statistical Software," http://gking.harvard.edu/zelig
check
Call: zelig(formula = Incident ~ Temperature + Pressure, model = "logit",
data = data.frame)
Coefficients:
(Intercept) Temperature Pressure
13.2924 -0.2287 0.0104
Degrees of Freedom: 22 Total (i.e. Null); 20 Residual
Null Deviance: 28.3
Residual Deviance: 18.8 AIC: 24.8
x.out.53 <- setx(object = check, Temperature = 53)
x.out.53
(Intercept) Temperature Pressure
1 1 53 152.2
x.out.31 <- setx(object = check, Temperature = 31)
x.out.31
(Intercept) Temperature Pressure
1 1 31 152.2
s.out <- sim(object = check, x = x.out.53, x1 =
x.out.31)
s.out$ev
NULL
1:36 p.m.
New subject: [gov2001-l] Using Zelig to check our results in 2.3-2.5
Assuming that everything up to the point where you extract your quantities
of interest from your sim object is correct, your problem seems to be how
you call the value from the sim.out object. You need to use sim.out$qi$...
where ... is the quantity of interest (e.g. s.out$qi$ev).
Hope this helps.
Miya
On Thu, Mar 12, 2009 at 11:48 AM, Christopher Michael Muller <
muller at fas.harvard.edu> wrote:
I tried this but for some reason got a
"NULL" result.
Quoting Patrick Lam <plam at fas.harvard.edu>:
> it might help to look at the Zelig manual, particularly the setx and sim
> functions:
> http://gking.harvard.edu/zelig/docs/index.html
> On Thu, Mar 12, 2009 at 2:23 AM, Lee,
Clarence <clee at hbs.edu> wrote:
> > Hi,
> > > > >
Maybe I?m missing something in my understanding, but how do you
evaluated
> > > >
Thanks!
> > > > Clarence
> > > > >
*From:* gov2001-l-bounces at lists.fas.harvard.edu [mailto:
> > gov2001-l-bounces at lists.fas.harvard.edu] *On Behalf Of *Z. A. Townsend
> > *Sent:* Thursday, March 12, 2009 2:14 AM
> > *To:* gov2001-l at lists.fas.harvard.edu
> > *Subject:* Re: [gov2001-l] Using Zelig to check our results in 2.3-2.5
> > > > > for
2.3, for example, something like
> > > data <-
as.data.frame(cbind(y,X)) %which I defined earlier to be the
vector
> > and matrix of interest
> > names(data) <- c("incident", "intercept",
"temperature", "pressure")
> > output <- zelig(incident ~ temperature +
> > pressure, data = data, model = "logit")
> > > > >
Should confirm your results.
> > > > > On
Thu, Mar 12, 2009 at 12:33 AM, Lee, Clarence <clee at hbs.edu> wrote:
> > > Does anyone know the
commands in Zelig to quickly check our results for
> > 2.3-2.5?
> > > > >
Thanks,
> > > > >
Clarence
> > > > _______________________________________________
> > gov2001-l mailing list
> > gov2001-l at lists.fas.harvard.edu
> > http://lists.fas.harvard.edu/mailman/listinfo/gov2001-l
> > > > >
_______________________________________________
> > gov2001-l mailing list
> > gov2001-l at lists.fas.harvard.edu
> > http://lists.fas.harvard.edu/mailman/listinfo/gov2001-l
> > > --
> Patrick Lam
> Department of Government and Institute for Quantitative Social Science,
> Harvard University
>
http://www.people.fas.harvard.edu/~plam<http://www.people.fas.harvard.ed…
_______________________________________________
gov2001-l mailing list
gov2001-l at lists.fas.harvard.edu
http://lists.fas.harvard.edu/mailman/listinfo/gov2001-l
--
Miya Woolfalk
Ph.D. Student
Harvard University
Government and Social Policy
check <- zelig(Incident ~ Temperature +
Pressure, model="logit",
+ data=data.frame)
How to cite this model in Zelig:
Kosuke Imai, Gary King, and Oliva Lau. 2008. "logit: Logistic Regression
for
Dichotomous Dependent Variables" in Kosuke Imai, Gary King, and Olivia Lau,
"Zelig: Everyone's Statistical Software," http://gking.harvard.edu/zelig
check
Call: zelig(formula = Incident ~ Temperature + Pressure, model = "logit",
data = data.frame)
Coefficients:
(Intercept) Temperature Pressure
13.2924 -0.2287 0.0104
Degrees of Freedom: 22 Total (i.e. Null); 20 Residual
Null Deviance: 28.3
Residual Deviance: 18.8 AIC: 24.8
x.out.53 <- setx(object = check, Temperature = 53)
x.out.53
(Intercept) Temperature Pressure
1 1 53 152.2
x.out.31 <- setx(object = check, Temperature =
31)
x.out.31
(Intercept) Temperature Pressure
1 1 31 152.2
s.out <- sim(object = check, x = x.out.53, x1
= x.out.31)
s.out$ev
NULL
> the difference in the expected probability
in the data when checking
for
> > Zelig? Do you mind describing how you specified the inputs y and X?
>
27 Apr
27 Apr
7:08 p.m.
New subject: [gov2001-l] CEM
If anyone is comfortable with CEM, I have a couple of questions:
1. After running CEM, how do I access the matched data for further analysis? Using
MatchIt, match.data calls the data. What is the equivalent in cem?
2. I am trying to refine my CEM by using k2k matching. I have tried this two ways:
Option 1:
CEMMatch <- cem(treatment="Treated", data=Test, drop=dropvar,
cutpoints=MyCuts, k2k=TRUE)
Option 2:
CEMMatch <- cem(treatment="Treated", data=Test, drop=dropvar,
cutpoints=MyCuts)
CEMMatch1 <- k2k(CEMMatch, Test, "euclidean", 1)
In both options, I am getting the results below. We get the same set of warnings when we
replicate the example in the R help file for cem. Do these warnings mean anything that
needs to be addressed?
There were 50 or more warnings (use warnings() to see the first 50)
warnings()
Warning messages:
1: In min(x, na.rm = TRUE) : no non-missing arguments to min; returning Inf
9:32 p.m.
New subject: [gov2001-l] CEM
CEM objects have a value strata which returns the vector of the stratum
number in which each observation belongs. This can be used to fairly easily
reconstruct the matched data. I don't know of any a more straightforward
way, but that should work.
As for the warnings- I don't know what the issue is. I tried it over here
and I get the same problem, which appears to be contingent on specifying
k2k=TRUE. Error notwithstanding the procedure seems to work just fine (it
does produce a 1:1 match and a full object), but I understand why you'd be
concerned. I'll leave the answer to better informed people...
Brandon
On Mon, Apr 27, 2009 at 7:08 PM, Schmidt, William <wschmidt at hbs.edu> wrote:
If anyone is comfortable with CEM, I have a couple of
questions:
1. After running CEM, how do I access the matched data for further
analysis? Using MatchIt, match.data calls the data. What is the equivalent
in cem?
2. I am trying to refine my CEM by using k2k matching. I have tried
this two ways:
Option 1:
CEMMatch <- cem(treatment="Treated", data=Test, drop=dropvar,
cutpoints=MyCuts, k2k=TRUE)
Option 2:
CEMMatch <- cem(treatment="Treated", data=Test, drop=dropvar,
cutpoints=MyCuts)
CEMMatch1 <- k2k(CEMMatch, Test, "euclidean", 1)
In both options, I am getting the results below. We get the same set of
warnings when we replicate the example in the R help file for cem. Do these
warnings mean anything that needs to be addressed?
There were 50 or more warnings (use warnings() to see the first 50)
warnings()
Warning messages:
1: In min(x, na.rm = TRUE) : no non-missing arguments to min; returning Inf
_______________________________________________
gov2001-l mailing list
gov2001-l at mailman.fas.harvard.edu
http://mailman.fas.harvard.edu/mailman/listinfo/gov2001-l
5483
days inactive
5532
days old
gov2001@lists.gking.harvard.edu
22 comments
10 participants
participants (10)
-
aallen＠hbs.edu -
bstewart＠fas.harvard.edu -
charlotte.cavaille＠gmail.com -
clee＠hbs.edu -
eageeva＠hotmail.com -
muller＠fas.harvard.edu -
plam＠fas.harvard.edu -
woolfalk＠fas.harvard.edu -
wschmidt＠hbs.edu -
ztownsend＠gmail.com