Re: [R] maximum-likelihood-estimation with mle()

[Ben Bolker]

peter dalgaard  gmail.com> writes:

> 
> You are being over-optimistic with your starting values, and/or 
> with constrains on the parameter space. 
> Your fit is diverging in sigma for some reason known
>  only to nonlinear-optimizer gurus...
> 
> For me, it works either to put in an explicit
>  constraint or to reparametrize with log(sigma).

  [snip]

  Another few strategies:

set.seed(101)
x <- 1:10
y <- 3*x - 1 + rnorm(length(x), mean=0, sd=0.5)

library("stats4")
nLL <- function(a, b, sigma) {
  -sum(dnorm(y, mean=a*x+b, sd=sigma, log=TRUE))
}

fit <- mle(nLL, start=list(a=0, b=0, sigma=1), nobs=length(y))


   Nelder-Mead is generally more robust than BFGS, although slower:


fit2 <- mle(nLL, start=list(a=0, b=0, sigma=1),
 method="Nelder-Mead",nobs=length(y))

   bbmle can be used to simplify things slightly (this repeats
Peter Dalgaard's transformation of sigma):

library("bbmle")
fit3 <- mle2(y~dnorm(mean=mu,sd=exp(logsigma)),
   parameters=list(mu~x),
   data=data.frame(x,y),
   start=list(mu=0, logsigma=0))

  You can also profile out the sigma:

## dnorm with sd profiled out
dnorm2 <- function(x,mean,log=FALSE) {
  ssq <- sum((x-mean)^2)
  dnorm(x,mean,sd=sqrt(ssq/length(x)),log=log)
}
fit4 <- mle2(y~dnorm2(mean=mu),
   parameters=list(mu~x),
   data=data.frame(x,y),
   start=list(mu=0))

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] maximum-likelihood-estimation with mle()

[Ronald Kölpin]

Hi everyone,

I have a problem with maximum-likelihood-estimation in the following
situation:

Assume a functional relation y = f(x) (the specific form of f should be
irrelevant). For my observations I assume (for simplicity) white noise,
such that hat(y_i) = f(x_i) + epsilon_i, with the epsilon_i iid N(0,
sigma^2). Y_i should then be N(f(x_i), sigma^2)-distributed and due to
the iid assumption the density of Y = (Y_1, ..., Y_n) is simply the
product of the individual densities, taking the log gives the the sum of
the log of individual densities.

I tried coding this in R with a simple example: f(x) = a*x + b (simple
linear regression). This way I wanted to compare the results from my
ml-estimation (specifying the log-likelihood manually and estimating
with mle()) with the results from using lm(y~x). In my example however
it doesn't work:

x <- 1:10
y <- 3*x - 1 + rnorm(length(x), mean=0, sd=0.5)

library("stats4")
nLL <- function(a, b, sigma)
{
  -sum(dnorm(y, mean=a*x+b, sd=sigma, log=TRUE))
}

fit <- mle(nLL1, start=list(a=0, b=0, sigma=1), nobs=length(y))

summary(lm(y~x))
summary(fit)

These should be the same but the aren't. I must have made some mistake
specifying the (negative) log-likehood (but I just don't see it). I also
actually don't care much (at the moment) for estimating sigma but I
don't know of a way to specify (and estimate) the (negative)
log-likelihood without estimating sigma.

Thanks a lot
and kind regards

Ronald Koelpin

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] maximum-likelihood-estimation with mle()

[peter dalgaard]

You are being over-optimistic with your starting values, and/or with constrains 
on the parameter space. 
Your fit is diverging in sigma for some reason known only to 
nonlinear-optimizer gurus...

For me, it works either to put in an explicit constraint or to reparametrize 
with log(sigma).

E.g.

> fit <- mle(nLL, start=list(a=0, b=0, sigma=1), method="L-BFGS-B", 
> lower=c(-Inf, -Inf, 0.001))
> summary(fit)
Maximum likelihood estimation

Call:
mle(minuslogl = nLL, start = list(a = 0, b = 0, sigma = 1), method = 
"L-BFGS-B", 
lower = c(-Inf, -Inf, 0.001))

Coefficients:
Estimate Std. Error
a  3.0293645 0.01904026
b -1.1516406 0.11814174
sigma  0.1729418 0.03866673

-2 log L: -6.717779 

or 

> nLL <- function(a, b, lsigma)
+  -sum(dnorm(y, mean=a*x+b, sd=exp(lsigma), log=TRUE))
> fit <- mle(nLL, start=list(a=0, b=0, lsigma=0))
> summary(fit)
Maximum likelihood estimation

Call:
mle(minuslogl = nLL, start = list(a = 0, b = 0, lsigma = 0))

Coefficients:
Estimate Std. Error
a   3.029365 0.01903975
b  -1.151642 0.11813856
lsigma -1.754827 0.22360675

-2 log L: -6.717779 

both of which reproduce lm() except for DF issues.

To fix sigma, use fixed=list(sigma=0.19) in mle(). This also stabilizes the 
convergence (which it blinking well should, since it is now a purely quadratic 
problem).

-pd


On 11 Sep 2015, at 14:20 , Ronald Kölpin  wrote:

> Hi everyone,
> 
> I have a problem with maximum-likelihood-estimation in the following
> situation:
> 
> Assume a functional relation y = f(x) (the specific form of f should be
> irrelevant). For my observations I assume (for simplicity) white noise,
> such that hat(y_i) = f(x_i) + epsilon_i, with the epsilon_i iid N(0,
> sigma^2). Y_i should then be N(f(x_i), sigma^2)-distributed and due to
> the iid assumption the density of Y = (Y_1, ..., Y_n) is simply the
> product of the individual densities, taking the log gives the the sum of
> the log of individual densities.
> 
> I tried coding this in R with a simple example: f(x) = a*x + b (simple
> linear regression). This way I wanted to compare the results from my
> ml-estimation (specifying the log-likelihood manually and estimating
> with mle()) with the results from using lm(y~x). In my example however
> it doesn't work:
> 
> x <- 1:10
> y <- 3*x - 1 + rnorm(length(x), mean=0, sd=0.5)
> 
> library("stats4")
> nLL <- function(a, b, sigma)
> {
>  -sum(dnorm(y, mean=a*x+b, sd=sigma, log=TRUE))
> }
> 
> fit <- mle(nLL1, start=list(a=0, b=0, sigma=1), nobs=length(y))
> 
> summary(lm(y~x))
> summary(fit)
> 
> These should be the same but the aren't. I must have made some mistake
> specifying the (negative) log-likehood (but I just don't see it). I also
> actually don't care much (at the moment) for estimating sigma but I
> don't know of a way to specify (and estimate) the (negative)
> log-likelihood without estimating sigma.
> 
> Thanks a lot
> and kind regards
> 
> Ronald Koelpin
> 
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

-- 
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Office: A 4.23
Email: pd@cbs.dk  Priv: pda...@gmail.com

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] maximum likelihood estimation

[pari hesabi]

Hello,As an example for Exponential distribution the MLE is got by this 
structure:t - rexp(100, 2)loglik - function(theta){ log(theta) - theta*t}a - 
maxLik(loglik, start=1)print(a)Exponential distribution has a simple 
loglikelihood function.  But if a new pdf has a more complicated form 
like:pr(N(2)=n)= integral( ((2-x)^n)*(exp(ax-2))) - integral (((5-ax)^n)), both 
integrals are defined over the interval(0,2) with respect to x.   I  also need 
to use the loglike of the form : [F log(pr(n))]=lnL  where F is the vector of 
observations and (n) is the vector of input for  the defined pmf.  Do you think 
how I can insert (define) the  pr(n) in the loop below?  loglik - function(a) 
sum(F*(log(pr(n)))?  n - c(0,1,2,3,4,5,6,7,8)  F-c(0,0,1,3,5,7,8,11,10) 
 loglik - function(a) sum(F*(log(pr(n)))??  re - maxLik (loglik, 
start=.5)  summary(re)I would be grateful if you let me know your idea.Best 
Regards,pari
[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] maximum likelihood estimation

[Arne Henningsen]

On 10 October 2014 08:04, pari hesabi statistic...@hotmail.com wrote:
 Hello,As an example for Exponential distribution the MLE is got by this 
 structure:t - rexp(100, 2)loglik - function(theta){ log(theta) - theta*t}a 
 - maxLik(loglik, start=1)print(a)Exponential distribution has a simple 
 loglikelihood function.  But if a new pdf has a more complicated form 
 like:pr(N(2)=n)= integral( ((2-x)^n)*(exp(ax-2))) - integral (((5-ax)^n)), 
 both integrals are defined over the interval(0,2) with respect to x.   I  
 also need to use the loglike of the form : [F log(pr(n))]=lnL  where F is the 
 vector of observations and (n) is the vector of input for  the defined pmf.  
 Do you think how I can insert (define) the  pr(n) in the loop below?  loglik 
 - function(a) sum(F*(log(pr(n)))?  n - c(0,1,2,3,4,5,6,7,8)  
 F-c(0,0,1,3,5,7,8,11,10)  loglik - function(a) sum(F*(log(pr(n)))?? 
  re - maxLik (loglik, start=.5)  summary(re)I would be grateful if you 
 let me know your idea.Best Regards,pari

 [...]
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

Your example is not reproducible:

 pr(N(2)=n)= integral( ((2-x)^n)*(exp(ax-2))) - integral (((5-ax)^n))
Error: unexpected '=' in pr(N(2)=

Please format your email in a way that one can easily copy the R code
to the R console.

Best regards,
Arne

-- 
Arne Henningsen
http://www.arne-henningsen.name

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] Maximum likelihood Estimation

[Parvin Dehghani]

maximum likelihood estimation
pari hesabi  
6:04 AM 
To: r-help@r-project.org
Hello,
As an example for Exponential distribution the MLE is got by this structure:
t - rexp(100, 2)
loglik - function(theta){ log(theta) - theta*t}
a - maxLik(loglik, start=1)
print(a)

Exponential distribution has a simple loglikelihood function.  But if a new pdf 
has a more complicated form like:
pr(N(2)=n)= integral( ((2-x)^n)*(exp(ax-2))) - integral (((5-ax)^n)), both 
integrals are defined over the interval(0,2) with respect to x.  
 I  also need to use the loglike of the form : [F log(pr(n))]=lnL  where F is 
the vector of observations and (n) is the vector of input for  the defined pmf. 
 Do you think how I can insert (define) the  pr(n) in the loop below? 
My question is related to:loglik - function(a) sum(F*(log(pr(n)))?



  n - c(0,1,2,3,4,5,6,7,8)
  F-c(0,0,1,3,5,7,8,11,10)
  loglik - function(a) sum(F*(log(pr(n)))??
  re - maxLik (loglik, start=.5)
  summary(re)

I would be grateful if you let me know your idea.
Best Regards,
pari

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] maximum likelihood estimation

[Arne Henningsen]

Dear Pari

On 7 October 2014 10:55, pari hesabi statistic...@hotmail.com wrote:
 HelloI am trying to estimate the parameter of a function by the Maximum 
 Likelihood Estimation method.If  the function is the difference between two 
 integrals: C-function(n){integrand3-function(x) {((2-x)^n)*(exp(ax-2))}cc- 
 integrate (integrand3,0,2)print(cc)}
 D-function(n){integrand4-function(x) {((5-ax)^n)}cc- integrate 
 (integrand4,0,2)print(cc)}
 f(n) = C(n) - D(n)
  I need to estimate parameter (a).   loglikelihood function is in the form of 
  sum[F(k) log(f(k))]=lnL. I am wondering how to introduce my logliklihood 
 function to the loop below.  Can anybody help me for correcting the following 
 loop? if there are some other packages better than this , please let me 
 know.Thank you,Diba
 n - c(0,1,2,3,4,5,6,7,8)
 F-c(0,0,1,3,5,7,8,11,10)
 loglik - function(a) sum(F*log(C(n)-D(n)))
 re - maxLik (loglik, start=.5)
 summary(re)

Unfortunately, I cannot reproduce your example, because I get an error
message about an unexpected symbol.

 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

Best regards,
Arne

-- 
Arne Henningsen
http://www.arne-henningsen.name

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] maximum likelihood estimation

[pari hesabi]

HelloI am trying to estimate the parameter of a function by the Maximum 
Likelihood Estimation method.If  the function is the difference between two 
integrals: C-function(n){integrand3-function(x) {((2-x)^n)*(exp(ax-2))}cc- 
integrate (integrand3,0,2)print(cc)}
D-function(n){integrand4-function(x) {((5-ax)^n)}cc- integrate 
(integrand4,0,2)print(cc)}
f(n) = C(n) - D(n)
 I need to estimate parameter (a).   loglikelihood function is in the form of  
sum[F(k) log(f(k))]=lnL. I am wondering how to introduce my logliklihood 
function to the loop below.  Can anybody help me for correcting the following 
loop? if there are some other packages better than this , please let me 
know.Thank you,Diba
n - c(0,1,2,3,4,5,6,7,8)
F-c(0,0,1,3,5,7,8,11,10)
loglik - function(a) sum(F*log(C(n)-D(n)))
re - maxLik (loglik, start=.5)
summary(re)   
[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Maximum likelihood estimation (stats4::mle)

[peter dalgaard]

On 22 Jul 2014, at 06:04 , David Winsemius dwinsem...@comcast.net wrote:

 
 On Jul 21, 2014, at 12:10 PM, Ronald Kölpin wrote:
 
 Dear R-Community,
 
 I'm trying to estimate the parameters of a probability distribution
 function by maximum likelihood estimation (using the stats4 function
 mle()) but can't seem to get it working.
 
 For each unit of observation I have a pair of observations (a, r)
 which I assume (both) to be log-normal distributed (iid). Taking the
 log of both values I have (iid) normally distributed random variables
 and the likelihood function to be estimated is:
 
 L = Product(F(x_i) - F(y_i), i=1..n)
 
 where F is the Normal PDF and (x,y) := (log(a), log(r)). Taking the
 log and multiplying by -1 gives the negative loglikelihood
 
 I don't see the need to multiply by -1. The log of any probability is of 
 necessity less than (or possibly equal to) 0 since probabilities are bounded 
 above by 1.


Well, mle() wants to minimize -log L by definition. That's just because 
optimizers tend to prefer to have an objective function to minimize rather than 
maximize.

However, what is keeping the terms of Product(F(x_i) - F(y_i), i=1..n) from 
going negative? As far as I can tell, the x_i are bigger than the corresponding 
y_i, and F is an increasing function so all terms are negative and the log of 
each term is undefined. Switching the order should help.

Also, what is the rationale for this form of likelihood? Surely not that the 
x,y pairs are independent lognormals. Looks more like what you'd get from 
interval censored data. If so, it should be possible to come up with better 
starting values than mu=0, sd=1.

-pd


 So sums of these number will be negative which then allows you to maximize 
 their sums.
 
 
 l = Sum(log( F(x_i) - F(y_i) ), i=1..n)
 
 However estimation by mle() produces the error vmmin is not finite
 
 
 As I would have predicted. If one maximizes numbers that get larger as 
 probabilities get small this is what would be expected.
 
 -- 
 David.
 
 and NaN have been created - even though put bound on the parameters
 mu and sigma (see code below).
 
 
 library(stats4)
 
 gaps - matrix(nrow=10, ncol=4, dimnames=list(c(1:10),c(r_i, a_i,
 log(r_i), log(a_i
 gaps[,1] - c(2.6, 1.4, 2.2, 2.9, 2.9, 1.7, 1.3, 1.7, 3.8, 4.5)
 gaps[,2] - c(9.8, 20.5, 8.7, 7.2, 10.3, 11, 4.5, 5.2, 6.7, 7.6)
 gaps[,3] - log(gaps[,1])
 gaps[,4] - log(gaps[,2])
 
 nll - function(mu, sigma)
 {
   if(sigma = 0  mu = 0)
   {
   -sum(log(pnorm(gaps[,3], mean=mu, sd=sigma) - pnorm(gaps[,4],
 mean=mu, sd=sigma)))
   }
   else
   {
   NA
   }
 }
 
 fit - mle(nll, start=list(mu=0, sigma=1), nobs=10)
 print(fit)
 
 
 To be honest, I'm stumped and don't quite know what the problem is...
 
 Regards and Thanks,
 
 Ronald Kölpin
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.
 
 David Winsemius
 Alameda, CA, USA
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

-- 
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Email: pd@cbs.dk  Priv: pda...@gmail.com

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Maximum likelihood estimation (stats4::mle)

[Ronald Kölpin]

Thanks, that was exactly it -- switching the values did the trick (and
was actually correct in terms of theory.) And of course, you are right
-- i changed the starting values to mean(x) - mean(y) for mu and
sqrt(var(x-y)) for sigma.

I also see your point about the theoretical justification for the
specific form of the likelihood function: F(x) - F(y) is supposed to
express the probability that a certain lies in the intervall (x,y),
that is x, y are takes as lower and upper bounds for the estimate
respectively. I think the theoretical justification is sound, but i'm
still trying to work out the details.

Thanks a lot anyways for solving my coding problem!

RK

On 22.07.2014 15:26, peter dalgaard wrote:
 
 On 22 Jul 2014, at 06:04 , David Winsemius dwinsem...@comcast.net
 wrote:
 
 
 On Jul 21, 2014, at 12:10 PM, Ronald Kölpin wrote:
 
 Dear R-Community,
 
 I'm trying to estimate the parameters of a probability
 distribution function by maximum likelihood estimation (using
 the stats4 function mle()) but can't seem to get it working.
 
 For each unit of observation I have a pair of observations (a,
 r) which I assume (both) to be log-normal distributed (iid).
 Taking the log of both values I have (iid) normally distributed
 random variables and the likelihood function to be estimated
 is:
 
 L = Product(F(x_i) - F(y_i), i=1..n)
 
 where F is the Normal PDF and (x,y) := (log(a), log(r)). Taking
 the log and multiplying by -1 gives the negative loglikelihood
 
 I don't see the need to multiply by -1. The log of any
 probability is of necessity less than (or possibly equal to) 0
 since probabilities are bounded above by 1.
 
 
 Well, mle() wants to minimize -log L by definition. That's just
 because optimizers tend to prefer to have an objective function to
 minimize rather than maximize.
 
 However, what is keeping the terms of Product(F(x_i) - F(y_i),
 i=1..n) from going negative? As far as I can tell, the x_i are
 bigger than the corresponding y_i, and F is an increasing function
 so all terms are negative and the log of each term is undefined.
 Switching the order should help.
 
 Also, what is the rationale for this form of likelihood? Surely not
 that the x,y pairs are independent lognormals. Looks more like what
 you'd get from interval censored data. If so, it should be possible
 to come up with better starting values than mu=0, sd=1.
 
 -pd
 
 
 So sums of these number will be negative which then allows you to
 maximize their sums.
 
 
 l = Sum(log( F(x_i) - F(y_i) ), i=1..n)
 
 However estimation by mle() produces the error vmmin is not
 finite
 
 
 As I would have predicted. If one maximizes numbers that get
 larger as probabilities get small this is what would be
 expected.
 
 -- David.
 
 and NaN have been created - even though put bound on the
 parameters mu and sigma (see code below).
 
 
 library(stats4)
 
 gaps - matrix(nrow=10, ncol=4, dimnames=list(c(1:10),c(r_i,
 a_i, log(r_i), log(a_i gaps[,1] - c(2.6, 1.4, 2.2,
 2.9, 2.9, 1.7, 1.3, 1.7, 3.8, 4.5) gaps[,2] - c(9.8, 20.5,
 8.7, 7.2, 10.3, 11, 4.5, 5.2, 6.7, 7.6) gaps[,3] -
 log(gaps[,1]) gaps[,4] - log(gaps[,2])
 
 nll - function(mu, sigma) { if(sigma = 0  mu = 0) { 
 -sum(log(pnorm(gaps[,3], mean=mu, sd=sigma) - pnorm(gaps[,4], 
 mean=mu, sd=sigma))) } else { NA } }
 
 fit - mle(nll, start=list(mu=0, sigma=1), nobs=10) print(fit)
 
 
 To be honest, I'm stumped and don't quite know what the problem
 is...
 
 Regards and Thanks,
 
 Ronald Kölpin
 
 __ 
 R-help@r-project.org mailing list 
 https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the
 posting guide http://www.R-project.org/posting-guide.html and
 provide commented, minimal, self-contained, reproducible code.
 
 David Winsemius Alameda, CA, USA
 
 __ 
 R-help@r-project.org mailing list 
 https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the
 posting guide http://www.R-project.org/posting-guide.html and
 provide commented, minimal, self-contained, reproducible code.


__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] Maximum likelihood estimation (stats4::mle)

[Ronald Kölpin]

Dear R-Community,

I'm trying to estimate the parameters of a probability distribution
function by maximum likelihood estimation (using the stats4 function
mle()) but can't seem to get it working.

For each unit of observation I have a pair of observations (a, r)
which I assume (both) to be log-normal distributed (iid). Taking the
log of both values I have (iid) normally distributed random variables
and the likelihood function to be estimated is:

L = Product(F(x_i) - F(y_i), i=1..n)

where F is the Normal PDF and (x,y) := (log(a), log(r)). Taking the
log and multiplying by -1 gives the negative loglikelihood

l = Sum(log( F(x_i) - F(y_i) ), i=1..n)

However estimation by mle() produces the error vmmin is not finite
and NaN have been created - even though put bound on the parameters
mu and sigma (see code below).


library(stats4)

gaps - matrix(nrow=10, ncol=4, dimnames=list(c(1:10),c(r_i, a_i,
log(r_i), log(a_i
gaps[,1] - c(2.6, 1.4, 2.2, 2.9, 2.9, 1.7, 1.3, 1.7, 3.8, 4.5)
gaps[,2] - c(9.8, 20.5, 8.7, 7.2, 10.3, 11, 4.5, 5.2, 6.7, 7.6)
gaps[,3] - log(gaps[,1])
gaps[,4] - log(gaps[,2])

nll - function(mu, sigma)
{
if(sigma = 0  mu = 0)
{
-sum(log(pnorm(gaps[,3], mean=mu, sd=sigma) - pnorm(gaps[,4],
mean=mu, sd=sigma)))
}
else
{
NA
}
}

fit - mle(nll, start=list(mu=0, sigma=1), nobs=10)
print(fit)


To be honest, I'm stumped and don't quite know what the problem is...

Regards and Thanks,

Ronald Kölpin

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Maximum likelihood estimation (stats4::mle)

[David Winsemius]

On Jul 21, 2014, at 12:10 PM, Ronald Kölpin wrote:

 Dear R-Community,
 
 I'm trying to estimate the parameters of a probability distribution
 function by maximum likelihood estimation (using the stats4 function
 mle()) but can't seem to get it working.
 
 For each unit of observation I have a pair of observations (a, r)
 which I assume (both) to be log-normal distributed (iid). Taking the
 log of both values I have (iid) normally distributed random variables
 and the likelihood function to be estimated is:
 
 L = Product(F(x_i) - F(y_i), i=1..n)
 
 where F is the Normal PDF and (x,y) := (log(a), log(r)). Taking the
 log and multiplying by -1 gives the negative loglikelihood

I don't see the need to multiply by -1. The log of any probability is of 
necessity less than (or possibly equal to) 0 since probabilities are bounded 
above by 1. So sums of these number will be negative which then allows you to 
maximize their sums.

 
 l = Sum(log( F(x_i) - F(y_i) ), i=1..n)
 
 However estimation by mle() produces the error vmmin is not finite


As I would have predicted. If one maximizes numbers that get larger as 
probabilities get small this is what would be expected.

-- 
David.

 and NaN have been created - even though put bound on the parameters
 mu and sigma (see code below).
 
 
 library(stats4)
 
 gaps - matrix(nrow=10, ncol=4, dimnames=list(c(1:10),c(r_i, a_i,
 log(r_i), log(a_i
 gaps[,1] - c(2.6, 1.4, 2.2, 2.9, 2.9, 1.7, 1.3, 1.7, 3.8, 4.5)
 gaps[,2] - c(9.8, 20.5, 8.7, 7.2, 10.3, 11, 4.5, 5.2, 6.7, 7.6)
 gaps[,3] - log(gaps[,1])
 gaps[,4] - log(gaps[,2])
 
 nll - function(mu, sigma)
 {
if(sigma = 0  mu = 0)
{
-sum(log(pnorm(gaps[,3], mean=mu, sd=sigma) - pnorm(gaps[,4],
 mean=mu, sd=sigma)))
}
else
{
NA
}
 }
 
 fit - mle(nll, start=list(mu=0, sigma=1), nobs=10)
 print(fit)
 
 
 To be honest, I'm stumped and don't quite know what the problem is...
 
 Regards and Thanks,
 
 Ronald Kölpin
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

David Winsemius
Alameda, CA, USA

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] Maximum likelihood estimation of ARMA(1,1)-GARCH(1,1)

[Andy Yeh]

Hello

Following some standard textbooks on ARMA(1,1)-GARCH(1,1) (e.g. Ruey
Tsay's Analysis of Financial Time Series), I try to write an R program
to estimate the key parameters of an ARMA(1,1)-GARCH(1,1) model for
Intel's stock returns. For some random reason, I cannot decipher what
is wrong with my R program. The R package fGarch already gives me the
answer, but my customized function does not seem to produce the same
result.

I would like to build an R program that helps estimate the baseline
ARMA(1,1)-GARCH(1,1) model. Then I would like to adapt this baseline
script to fit different GARCH variants (e.g. EGARCH, NGARCH, and
TGARCH). It would be much appreciated if you could provide some
guidance in this case. The code below is the R script for estimating
the 6 parameters of an ARMA(1,1)-GARCH(1,1) model for Intel's stock
returns. At any rate, I would be glad to know your thoughts and
insights. If you have a similar example, please feel free to share
your extant code in R. Many thanks in advance.

Emily



# This R script offers a suite of functions for estimating  the
volatility dynamics based on the standard ARMA(1,1)-GARCH(1,1) model
and its variants.
# The baseline ARMA(1,1) model characterizes the dynamic evolution of
the return generating process.
# The baseline GARCH(1,1) model depicts the the return volatility
dynamics over time.
# We can extend the GARCH(1,1) volatility model to a variety of
alternative specifications to capture the potential asymmetry for a
better comparison:
# GARCH(1,1), EGARCH(1,1),  NGARCH(1,1), and TGARCH(1,1).

options(scipen=10)

intel= read.csv(file=intel.csv)
summary(intel)

raw_data= as.matrix(intel$logret)

library(fGarch)
garchFit(~arma(1,1)+garch(1,1), data=raw_data, trace=FALSE)


negative_log_likelihood_arma11_garch11=
function(theta, data)
{mean =theta[1]
 delta=theta[2]
 gamma=theta[3]
 omega=theta[4]
 alpha=theta[5]
 beta= theta[6]

 r= ts(data)
 n= length(r)

 u= vector(length=n)
 u= ts(u)
 u[1]= r[1]- mean

 for (t in 2:n)
 {u[t]= r[t]- mean- delta*r[t-1]- gamma*u[t-1]}

 h= vector(length=n)
 h= ts(h)
 h[1]= omega/(1-alpha-beta)

 for (t in 2:n)
 {h[t]= omega+ alpha*(u[t-1]^2)+ beta*h[t-1]}

 #return(-sum(dnorm(u[2:n], mean=mean, sd=sqrt(h[2:n]), log=TRUE)))
 pi=3.141592653589793238462643383279502884197169399375105820974944592
 return(-sum(-0.5*log(2*pi) -0.5*log(h[2:n]) -0.5*(u[2:n]^2)/h[2:n]))
}


#theta0=c(0, +0.78, -0.79, +0.018, +0.06, +0.93, 0.01)
theta0=rep(0.01,6)
negative_log_likelihood_arma11_garch11(theta=theta0, data=raw_data)


alpha= proc.time()
maximum_likelihood_fit_arma11_garch11=
nlm(negative_log_likelihood_arma11_garch11,
p=theta0,
data=raw_data,
hessian=TRUE,
iterlim=500)
#optim(theta0,
#  negative_log_likelihood_arma11_garch11,
#  data=raw_data,
#  method=L-BFGS-B,
#  
upper=c(+0.,+0.,+0.,0.,0.,0.),
#  
lower=c(-0.,-0.,-0.,0.0001,0.0001,0.0001),
#  hessian=TRUE)

# We record the end time and calculate the total runtime for the above work.
omega= proc.time()
runtime= omega-alpha
zhours = floor(runtime/60/60)
zminutes=floor(runtime/60- zhours*60)
zseconds=floor(runtime- zhours*60*60- zminutes*60)
print(paste(It takes ,zhours,hour(s), zminutes, minute(s) ,and
, zseconds,second(s) to finish running this R program,sep=))


maximum_likelihood_fit_arma11_garch11

sqrt(diag(solve(maximum_likelihood_fit_arma11_garch11$hessian)))

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] maximum likelihood estimation in R

[Ben Bolker]

David Winsemius dwinsemius at comcast.net writes:

 
 
 On Nov 10, 2012, at 9:22 PM, mmosalman wrote:
 
  I want to find ML estimates of a model using mle2 in bbmle package. When I
  insert new parameters (for new covariates) in model the log-likelihood value
  does not change and the estimated value is exactly the initial value that I
  determined. What's the problem? This is the code and the result:
 
 Nope. The code is not here. It may be visible in the Nabble universe
  but in the Rhelp universe it is completely missing.
 
  As you see the estimated values for b2 , b3 and b4 are the initial values of
  them. The log-likelihood value did not change!

  Copying and pasting the relevant code from Nabble:

library(GB2) 
library(bbmle)

lgb1=function(a,b0,b2,b3,b4,p,q){
 xb=b0+b2*fsex1[,2]+b3*fvtype1[,2]+b4*fvuse1[,2]
ll=sum(log(dgb2(loss1,a,exp(xb),p,q)))
 return(-ll)
 }

start=list(a=3.1,b0=2.5,
b2=.2,b3=1,b4=-.5,
p=7.2,q=.3)

mle2(lgb1,start)-fit1

summary(fit1)
Maximum likelihood estimation

Call:
mle2(minuslogl = lgb1, start = start)

Coefficients:
  Estimate  Std. Error z value Pr(z)
a   3.0747e+00  6.4741e-01  4.7492e+00 2.042e-06 ***
b0  2.5327e+00  4.6887e-01  5.4016e+00 6.605e-08 ***
b2  2.e-01  3.9686e-11  5.0396e+09  2.2e-16 ***
b3  1.e+00  7.6565e-12  1.3061e+11  2.2e-16 ***
b4 -5.e-01  1.5312e-11 -3.2653e+10  2.2e-16 ***
p   7.1281e+00  8.0269e+00  8.8800e-010.3745
q   3.5098e-01  8.6902e-02  4.0388e+00 5.372e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

-2 log L: 10137.56


  I agree that it looks like the parameters got stuck at
their initial values, but without the data it's really hard to
know what went wrong.  I assume you started with reasonable
parameter values and did some sanity checking on the data ... ?

  You can use the 'trace' and 'browse_obj' arguments to 
get more detail about what's going on.

  For what it's worth, with a little more effort you could
use the formula interface to mle2: something like

## set up a d* function with a log argument
  dgb2B - function(...,log=FALSE) {
 r - dgb2(...)
 if (log) log(r) else r
  }

## set up a data frame to hold the data
dframe - data.frame(loss1,sex=fsex1[,2],type=fvtype1[,2],use=fvuse1[,2])

mle2(loss1~dgb2B(shape1=a,scale=exp(logscale),shape2=p,shape3=q),
   parameters=list(logscale~sex+type+use),
   data=dframe,
   start=...)

[see the help page for the details of how start is specified
in this case]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] maximum likelihood estimation in R

[mmosalman]

I want to find ML estimates of a model using mle2 in bbmle package. When I
insert new parameters (for new covariates) in model the log-likelihood value
does not change and the estimated value is exactly the initial value that I
determined. What's the problem? This is the code and the result:



As you see the estimated values for b2 , b3 and b4 are the initial values of
them. The log-likelihood value did not change!



--
View this message in context: 
http://r.789695.n4.nabble.com/maximum-likelihood-estimation-in-R-tp4649226.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] maximum likelihood estimation in R

[David Winsemius]

On Nov 10, 2012, at 9:22 PM, mmosalman wrote:

 I want to find ML estimates of a model using mle2 in bbmle package. When I
 insert new parameters (for new covariates) in model the log-likelihood value
 does not change and the estimated value is exactly the initial value that I
 determined. What's the problem? This is the code and the result:

Nope. The code is not here. It may be visible in the Nabble universe but in the 
Rhelp universe it is completely missing.

 As you see the estimated values for b2 , b3 and b4 are the initial values of
 them. The log-likelihood value did not change!
 
 --
 View this message in context: 
 http://r.789695.n4.nabble.com/maximum-likelihood-estimation-in-R-tp4649226.html
 Sent from the R help mailing list archive at Nabble.com.

Wrong. Nabble lies. It is not the Rhelp Archive and this message was sent from 
the Mailing list server, not from Nabble.

Please to pay attention to this set of resources.
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

Well you apparently did try, so that's something.

-- 

David Winsemius, MD
Alameda, CA, USA

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] Maximum Likelihood estimation of KB distribution

[chamilka]

Hi, The following distribution is known as Kumaraswamy binomial Distribution.
http://r.789695.n4.nabble.com/file/n4636782/kb.png 

For a given data I need to estimate the paramters (alpha and beta) of this
distribution(Known as Kumaraswamy binomial Distribution, A Binomial Like
Distribution). For that, in order to use *optim()*, I first declared the
Negative Log-likelihood of this distribution as follows;

*Loglik.newdis2-function(x,a,b,n,imax) { 
  term-0 
  for (i in 0:imax) { 
   term=term+(((-1)**i)*(choose(b-1,i))*(beta(x+a+a*i,n-x+1))) 
   } 
  dens=a*b*choose(n,x)*term 
  KBLL2-sum(log(dens)) 
  return(-KBLL2) 
  }  *

since there is an infinite convergent series in this PMF, I decided to
specify a maximum value as imax instead of infinity without loss of any
information, and n is the binomial trials. 

Please tell me whether the declared negative loglikelihood (NLL) is correct
for this distribution?  

***I couldn't get the result of this
paper(http://aps.ecnu.edu.cn/EN/abstract/abstract8722.shtml) with my NLL

Thanks in advance. 

--
View this message in context: 
http://r.789695.n4.nabble.com/Maximum-Likelihood-estimation-of-KB-distribution-tp4636782.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] Maximum Likelihood Estimation Poisson distribution mle {stats4}

[chamilka]

Hi everyone!
I am using the mle {stats4} to estimate the parameters of distributions by
MLE method. I have a problem with the examples they provided with the
mle{stats4} html files. Please check the example and my question below!
*Here is the mle html help file * 
http://stat.ethz.ch/R-manual/R-devel/library/stats4/html/mle.html
http://stat.ethz.ch/R-manual/R-devel/library/stats4/html/mle.html 

*In the example provided with the help *
 od - options(digits = 5)
 x - 0:10*#generating Poisson counts*
 y - c(26, 17, 13, 12, 20, 5, 9, 8, 5, 4, 8)  *#generating the
 frequesncies*
 
 ## Easy one-dimensional MLE:
 nLL - function(lambda) -sum(stats::dpois(y, lambda, log=TRUE))  *#they
 define the Poisson negative loglikelihood*
 fit0 - mle(nLL, start = list(lambda = 5), nobs = NROW(y)) * #they
 estimate the Poisson parameter using mle*
 fit0  *#they call the output*

Call:
mle(minuslogl = nLL, start = list(lambda = 5), nobs = NROW(y))

Coefficients:
lambda 
11.545   * #this is their estimated Lambda Vallue.*

*Now my question is in a Poisson distribution the Maximum Likelihood
estimator of the mean parameter lambda is  the sample mean, so if we
calculate the sample mean of that generated Poisson distribution manually
using R we get  the below!*

 sample.mean- sum(x*y)/sum(y)
 sample.mean
[1] 3.5433

*This is the contradiction!! *
Here I am getting the estimate as 3.5433(which is reasonable as most of the
values are clustered around 3), but mle code gives the estimate 11.545(which
may not be correct as this is out side the range 0:10)

Why this contradiction? 

--
View this message in context: 
http://r.789695.n4.nabble.com/Maximum-Likelihood-Estimation-Poisson-distribution-mle-stats4-tp4635464.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Maximum Likelihood Estimation Poisson distribution mle {stats4}

[S Ellison]

 -Original Message-
  sample.mean- sum(x*y)/sum(y)
  sample.mean
 [1] 3.5433
 
 *This is the contradiction!! *
 Here I am getting the estimate as 3.5433(which is reasonable 
 as most of the values are clustered around 3), but mle code 
 gives the estimate 11.545(which may not be correct as this is 
 out side the range 0:10)
 
 Why this contradiction? 
 
 

Ermm.. the sample mean is mean(y) which is 11.545

I deduce that you and the help page have different views on what the sample y 
represents.

i also note that x does not figure at all in the help page's log-likelihood, 
suggesting that y is a simple set of counts, whereas you have interpreted x as 
the number of instances of each y. That appears not to be the case.

S


 

***
This email and any attachments are confidential. Any use...{{dropped:8}}

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Maximum Likelihood Estimation Poisson distribution mle {stats4}

[peter dalgaard]

On Jul 5, 2012, at 10:48 , chamilka wrote:

 Hi everyone!
 I am using the mle {stats4} to estimate the parameters of distributions by
 MLE method. I have a problem with the examples they provided with the
 mle{stats4} html files. Please check the example and my question below!
 *Here is the mle html help file * 
 http://stat.ethz.ch/R-manual/R-devel/library/stats4/html/mle.html
 http://stat.ethz.ch/R-manual/R-devel/library/stats4/html/mle.html 
 
 *In the example provided with the help *
 od - options(digits = 5)
 x - 0:10*#generating Poisson counts*
 y - c(26, 17, 13, 12, 20, 5, 9, 8, 5, 4, 8)  *#generating the
 frequesncies*
 
 ## Easy one-dimensional MLE:
 nLL - function(lambda) -sum(stats::dpois(y, lambda, log=TRUE))  *#they
 define the Poisson negative loglikelihood*
 fit0 - mle(nLL, start = list(lambda = 5), nobs = NROW(y)) * #they
 estimate the Poisson parameter using mle*
 fit0  *#they call the output*
 
 Call:
 mle(minuslogl = nLL, start = list(lambda = 5), nobs = NROW(y))
 
 Coefficients:
 lambda 
 11.545   * #this is their estimated Lambda Vallue.*
 
 *Now my question is in a Poisson distribution the Maximum Likelihood
 estimator of the mean parameter lambda is  the sample mean, so if we
 calculate the sample mean of that generated Poisson distribution manually
 using R we get  the below!*
 
 sample.mean- sum(x*y)/sum(y)
 sample.mean
 [1] 3.5433
 
 *This is the contradiction!! *
 Here I am getting the estimate as 3.5433(which is reasonable as most of the
 values are clustered around 3), but mle code gives the estimate 11.545(which
 may not be correct as this is out side the range 0:10)
 
 Why this contradiction? 

You misunderstand the setup. y are 11 Poisson counts, x is a dose variable used 
further down in the example. So

 mean(y)
[1] 11.54545

Your sample.mean would be relevant if there were 127 count measurements, 26 of 
which were 0, 17 were 1, etc. But then the likelihood would be different:

 x - 0:10
 nLL - function(lambda) -sum(y*stats::dpois(x, lambda, log=TRUE))
 mle(nLL, start = list(lambda = 5), nobs = NROW(y))

Call:
mle(minuslogl = nLL, start = list(lambda = 5), nobs = NROW(y))

Coefficients:
 lambda 
3.54328 




 
 --
 View this message in context: 
 http://r.789695.n4.nabble.com/Maximum-Likelihood-Estimation-Poisson-distribution-mle-stats4-tp4635464.html
 Sent from the R help mailing list archive at Nabble.com.
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

-- 
Peter Dalgaard, Professor
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Email: pd@cbs.dk  Priv: pda...@gmail.com

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Maximum Likelihood Estimation Poisson distribution mle {stats4}

[chamilka]

Thank you S Ellison-2 for your reply. I will understand it with Prof.Peter
Dalgaard's answer.. 

--
View this message in context: 
http://r.789695.n4.nabble.com/Maximum-Likelihood-Estimation-Poisson-distribution-mle-stats4-tp4635464p4635484.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Maximum Likelihood Estimation Poisson distribution mle {stats4}

[chamilka]

Thank you very much Professor .Peter Dalgaard for your kind explanations..
This made my work easy.. I am struggling with this for more than 2 days and
now I got the correct reply. 
Thank again. 

--
View this message in context: 
http://r.789695.n4.nabble.com/Maximum-Likelihood-Estimation-Poisson-distribution-mle-stats4-tp4635464p4635485.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] Maximum likelihood estimation in R

[Maomao Luo]

Dear R-helper,
 
I am trying to do maximum likelihood estimation in R.  I use the optim 
function.  Since I have no prior information on the true values of the 
parameters, I just randomly select different sets of starting values to feed 
into the program.  Each time, I get the following error message: Error in 
optim(theta0, lf, method = BFGS, hessian = T, Y = Y, X = X,  : non-finite 
finite-difference value [6]
In addition: There were 50 or more warnings (use warnings() to see the first 50)

I guess this has something to do with the starting values that I chose.  But I 
have no idea as to how to pick better starting values in a case like mine or 
how to get rid of this error message.  Could you give me some suggestion?  
Thanks a lot for your kind help.
  
 
My R code is as follows.  If more information is needed, please let me know.

 




lf-function(theta,Y,X,G)   {
N-26*26
l-matrix(1,26,1)

alpha-theta[1:4]
betad1-t(t(theta[5:6]))
betad2-t(t(theta[7:8]))
betad3-t(t(theta[9:10]))
betad4-t(t(theta[11:12]))
betao1-t(t(theta[13:14]))
betao2-t(t(theta[15:16]))
betao3-t(t(theta[17:18]))
betao4-t(t(theta[19:20]))
gamma-theta[21:24]
pd-theta[25]
po-theta[26]
pw-theta[27]

E1-  Y-alpha[1]*l%*%t(l)-X%*%betad1%*%t(l)-l%*%t(betao1)%*%t(X)-gamma[1]*G
E2- WY-alpha[2]*l%*%t(l)-X%*%betad2%*%t(l)-l%*%t(betao2)%*%t(X)-gamma[2]*G
E3- YW-alpha[3]*l%*%t(l)-X%*%betad3%*%t(l)-l%*%t(betao3)%*%t(X)-gamma[3]*G
E4-WYW-alpha[4]*l%*%t(l)-X%*%betad4%*%t(l)-l%*%t(betao4)%*%t(X)-gamma[4]*G
 
Q-matrix(0,4,4)
Q[1,1]-sum(E1*t(E1))
Q[1,2]-sum(E1*t(E2))
Q[1,3]-sum(E1*t(E3))
Q[1,4]-sum(E1*t(E4))
Q[2,1]-sum(E2*t(E1))
Q[2,2]-sum(E2*t(E2))
Q[2,3]-sum(E2*t(E3))
Q[2,4]-sum(E2*t(E4))
Q[3,1]-sum(E3*t(E1))
Q[3,2]-sum(E3*t(E2))
Q[3,3]-sum(E3*t(E3))
Q[3,4]-sum(E3*t(E4))
Q[4,1]-sum(E4*t(E1))
Q[4,2]-sum(E4*t(E2))
Q[4,3]-sum(E4*t(E3))
Q[4,4]-sum(E4*t(E4))
 
T-matrix(c(1, -pd, -po, -pw),nrow=4)
 
logl- log(det(diag(1,N)-pd*wd - po*wo -pw*ww))-N/2*log(t(T)%*%Q%*%T)
return(-logl)

   }
 
theta0-c(0,0,0,0, 0,0,0,0,0,0,0,0,  0,0,0,0,0,0,0,0,   0,0,0,0, 0.2, 0.1, 0)

p-optim(theta0,lf,method=BFGS, hessian=T, Y=Y, X=X, G=G)
 
OI-solve(p$hessian)
 




Hotmail is redefining busy with tools for the New Busy. Get more from your 
inbox. See how.
_
Hotmail is redefining busy with tools for the New Busy. Get more from your 
inbox.

N:WL:en-US:WM_HMP:042010_2
[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Maximum Likelihood Estimation in R

[Henkep]

that worked out well, thank you again! 

I also tried to use lm, and as expected in this case, I almost got the same
estimates of the parameters as in the MLE-case. 

Best Regards

Henrik
-- 
View this message in context: 
http://r.789695.n4.nabble.com/Maximum-Likelihood-Estimation-in-R-tp2018822p2022807.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Maximum Likelihood Estimation in R

[Henkep]

Abhishek:
Thank you!

Thomas:
that worked out well, thank you again! 

I also tried to use lm, and as expected in this case, I almost got the same
estimates of the parameters as in the MLE-case. 





Best Regards 

Henrik 
-- 
View this message in context: 
http://r.789695.n4.nabble.com/Maximum-Likelihood-Estimation-in-R-tp2018822p2022840.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Maximum Likelihood Estimation in R

[Henkep]

Thank you! 

Best Regards 

Henrik
-- 
View this message in context: 
http://r.789695.n4.nabble.com/Maximum-Likelihood-Estimation-in-R-tp2018822p2022832.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] Maximum Likelihood Estimation in R

[Henkep]

Dear R-Help,
 
I also send the following post by e-mail to you, however I try to post it
here aswell. My name is Henrik and I am currently trying to solve a Maximum
Likelihood optimization problem in R. Below you can find the output from R,
when I use the BFGS method:
 
The problem is that the parameters that I get are very unreasonable, I would
expect the absolute value of each parameter to be bounded by say 5.
(furthermore the variable stdev should be greater than zero).
 
One of the problems seems to be that I need to bound the stdev-variable from
below by zero to avoid the  NaN:s produced. I unfortunately do not know how
to do that.
 
Below y is the dataset, to which, I want to fit the parameters. I.e. y is
the vector (or R equivalent) of observations. 
I have also multiplied the log-likelihood function by -1, since I know that
R by default minimizes the objective function. 
 
I would be very happy if you can come up with some Ideas on what is going
wrong in the code below.
 
Thank you for your time!
 
Henrik
 
 
 data-read.table(file=c:/users/oukal.csv,header=TRUE)
 y-data$V1
 loglik-function(estimate,y)
+ {
+ lap-estimate[1]
+ stdev-estimate[2]
+ rev-estimate[3]
+ n-length(y)
+ 
+ 0.5*n*log(2*pi)+ 0.5*n*log(stdev) +
+ (1/(2*stdev*stdev))*sum(y-(rev/12)-lag(y)*exp(-lap/12))
+ }
 
 optim(c(4,4,4),loglik,y=y,method=BFGS)
$par
[1] -4884.34155   445.5235088.53777
$value
[1] -1.910321e+174
$counts
function gradient 
 2907 
$convergence
[1] 0
$message
NULL
There were 50 or more warnings (use warnings() to see the first 50)

-- 
View this message in context: 
http://n4.nabble.com/Maximum-Likelihood-Estimation-in-R-tp2018822p2018822.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Maximum Likelihood Estimation in R

[Thomas Stewart]

Two possible problems:

(a) If you're working with a normal likelihood---and it seems that you
are---the exponent should be squared.  As in:
... + (1/(2*stdev*stdev))*sum(  ( y-(rev/12)-lag(y)*exp(-lap/12) )^2  )

(b) lag may not be working like you think it should.  Consider this silly
example:
y-c(1,2,3,5,34,2,3,5,2)
lag(y)-y
lag(ts(y))-ts(y)

Good luck,
-tgs



On Wed, Apr 21, 2010 at 8:14 AM, Henkep flyerhe...@hotmail.com wrote:


 Dear R-Help,

 I also send the following post by e-mail to you, however I try to post it
 here aswell. My name is Henrik and I am currently trying to solve a Maximum
 Likelihood optimization problem in R. Below you can find the output from R,
 when I use the BFGS method:

 The problem is that the parameters that I get are very unreasonable, I
 would
 expect the absolute value of each parameter to be bounded by say 5.
 (furthermore the variable stdev should be greater than zero).

 One of the problems seems to be that I need to bound the stdev-variable
 from
 below by zero to avoid the  NaN:s produced. I unfortunately do not know how
 to do that.

 Below y is the dataset, to which, I want to fit the parameters. I.e. y is
 the vector (or R equivalent) of observations.
 I have also multiplied the log-likelihood function by -1, since I know that
 R by default minimizes the objective function.

 I would be very happy if you can come up with some Ideas on what is going
 wrong in the code below.

 Thank you for your time!

 Henrik


  data-read.table(file=c:/users/oukal.csv,header=TRUE)
  y-data$V1
  loglik-function(estimate,y)
 + {
 + lap-estimate[1]
 + stdev-estimate[2]
 + rev-estimate[3]
 + n-length(y)
 +
 + 0.5*n*log(2*pi)+ 0.5*n*log(stdev) +
 + (1/(2*stdev*stdev))*sum(y-(rev/12)-lag(y)*exp(-lap/12))
 + }
 
  optim(c(4,4,4),loglik,y=y,method=BFGS)
 $par
 [1] -4884.34155   445.5235088.53777
 $value
 [1] -1.910321e+174
 $counts
 function gradient
 2907
 $convergence
 [1] 0
 $message
 NULL
 There were 50 or more warnings (use warnings() to see the first 50)

 --
 View this message in context:
 http://n4.nabble.com/Maximum-Likelihood-Estimation-in-R-tp2018822p2018822.html
 Sent from the R help mailing list archive at Nabble.com.

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.


[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Maximum Likelihood Estimation in R

[Abhishek Pratap]

Hey Henrik

I dont do MLE myself but this recent blog might be helpful.

http://www.johnmyleswhite.com/notebook/2010/04/21/doing-maximum-likelihood-estimation-by-hand-in-r/

-A

On Wed, Apr 21, 2010 at 10:02 AM, Thomas Stewart tgstew...@gmail.com wrote:
 Two possible problems:

 (a) If you're working with a normal likelihood---and it seems that you
 are---the exponent should be squared.  As in:
 ... + (1/(2*stdev*stdev))*sum(  ( y-(rev/12)-lag(y)*exp(-lap/12) )^2  )

 (b) lag may not be working like you think it should.  Consider this silly
 example:
 y-c(1,2,3,5,34,2,3,5,2)
 lag(y)-y
 lag(ts(y))-ts(y)

 Good luck,
 -tgs



 On Wed, Apr 21, 2010 at 8:14 AM, Henkep flyerhe...@hotmail.com wrote:


 Dear R-Help,

 I also send the following post by e-mail to you, however I try to post it
 here aswell. My name is Henrik and I am currently trying to solve a Maximum
 Likelihood optimization problem in R. Below you can find the output from R,
 when I use the BFGS method:

 The problem is that the parameters that I get are very unreasonable, I
 would
 expect the absolute value of each parameter to be bounded by say 5.
 (furthermore the variable stdev should be greater than zero).

 One of the problems seems to be that I need to bound the stdev-variable
 from
 below by zero to avoid the  NaN:s produced. I unfortunately do not know how
 to do that.

 Below y is the dataset, to which, I want to fit the parameters. I.e. y is
 the vector (or R equivalent) of observations.
 I have also multiplied the log-likelihood function by -1, since I know that
 R by default minimizes the objective function.

 I would be very happy if you can come up with some Ideas on what is going
 wrong in the code below.

 Thank you for your time!

 Henrik


  data-read.table(file=c:/users/oukal.csv,header=TRUE)
  y-data$V1
  loglik-function(estimate,y)
 + {
 + lap-estimate[1]
 + stdev-estimate[2]
 + rev-estimate[3]
 + n-length(y)
 +
 + 0.5*n*log(2*pi)+ 0.5*n*log(stdev) +
 + (1/(2*stdev*stdev))*sum(y-(rev/12)-lag(y)*exp(-lap/12))
 + }
 
  optim(c(4,4,4),loglik,y=y,method=BFGS)
 $par
 [1] -4884.34155   445.52350    88.53777
 $value
 [1] -1.910321e+174
 $counts
 function gradient
     290        7
 $convergence
 [1] 0
 $message
 NULL
 There were 50 or more warnings (use warnings() to see the first 50)

 --
 View this message in context:
 http://n4.nabble.com/Maximum-Likelihood-Estimation-in-R-tp2018822p2018822.html
 Sent from the R help mailing list archive at Nabble.com.

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.


        [[alternative HTML version deleted]]

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.


__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Maximum Likelihood Estimation in R

[Henkep]

Thank you Thomas. 

(a) an embarrassing mistake by me. Of course it should be squared. Thank you
for pointing that out. 
(b) Do you possibly have any suggestions on how to solve this issue?  I
presume that there is no reason in trying to create a lagged vector
manually?

Best Regards

Henrik
-- 
View this message in context: 
http://n4.nabble.com/Maximum-Likelihood-Estimation-in-R-tp2018822p2019701.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Maximum Likelihood Estimation in R

[Thomas Stewart]

Henrik-

A coding solutions may be

... + (1/(2*stdev*stdev))*sum(  ( y-(rev/12)- c(0,y[-n]) *exp(-lap/12) )^2
 )

where n is the number of observations in y.

Personally, I would use lm.  Your model can be written as a linear function.


Let x=c(0,y[-n]).  Then run  lm(y~x).  The parameter estimates in y~x
correspond to the parameters in your model (with some very minor
arithmetic).

The lm() approach is quick and easy.  If, however, your goal is to gain
experience with MLEs by hand, I fully support the manual option.  You could
do both and compare the answers.

-tgs

On Wed, Apr 21, 2010 at 5:12 PM, Henkep flyerhe...@hotmail.com wrote:


 Thank you Thomas.

 (a) an embarrassing mistake by me. Of course it should be squared. Thank
 you
 for pointing that out.
 (b) Do you possibly have any suggestions on how to solve this issue?  I
 presume that there is no reason in trying to create a lagged vector
 manually?

 Best Regards

 Henrik
 --
 View this message in context:
 http://n4.nabble.com/Maximum-Likelihood-Estimation-in-R-tp2018822p2019701.html
 Sent from the R help mailing list archive at Nabble.com.

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.


[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] Maximum Likelihood Estimation

[Prof. John C Nash]

What do you mean by the estimates were very bad? In nearly 40 years of
working with optimization, I've seen badly set-up  functions cause
troubles, I've seen multiple minima situations, I've seen comparisons of
results from one data set to the estimates for another, and I've seen
optimization produce the proper results when published ones were just
plain wrong.

So there are many possibilities.

And for Cobb-Douglas functions, I've seen least squares estimates based
on the log of the production function give estimates with different
signs from those using a nonlinear least squares objective on the
original exponential form.

You are not the first to step in the cow pie of Cobb Douglas.

We need some numbers please. If you send a file or files (I don't think
you can attach files to Rhelp, but can do so off-list) with a runnable
script and the data, I'll be willing to try it out.

JN





 Message: 84
 Date: Tue, 3 Nov 2009 19:49:17 +
 From: Andre Barbosa Oliveira andreab...@hotmail.com
 Subject: [R] Maximum Likelihood Estimation
 To: r-help@r-project.org
 Message-ID: blu109-w2d5a97dbfd61d3976e589c2...@phx.gbl
 Content-Type: text/plain
 
 
 Hi,
 
 I would like estimate a model for function of production's Coob-Douglas using 
 maximum likelihood. The model is log(Y)= 
 beta[1]+beta[2]*log(L)+beta[3]*log(K). I tried estimate this model using the 
 tools nlm ( ) and optim ( ) using the log-likelihood function below:

  mloglik - function (beta, Y, L, K) {
 + n - length(Y)
 + sum ( (log(Y)- beta[1]-beta[2]*log(L)-beta[3]*log(K))2)/2*beta[4]2 + 
 n/2*log(2*pi)+ n*log(beta[4])
 + }
 Then I did estimates the parameters using nlm ( ) and optim ( ), but the 
 estimates were very bad. I used these codes:

  mlem - nlm (mloglik, c(1,1,1,1), Y=Y, L=L, K=K)
 
  mlem2 - optim(c(1,1,1,1), mloglik, Y=Y, L=L, K=K, method=BFGS)
 
 How I improve the estimates What's the best and more simple form for 
 estimate a modelo using the  maximum likelihood's method???
 
 Best regards,
 
 Andre' Barbosa Oliveira
 
 Student of Master in Economics at University Federal of Rio Grande do Sul - 
 Brazil 
 
 _
 Novo windowslive.com.br. Descubra como juntar a galera com os produtos 
 Windows Live.
 
 
   [[alternative HTML version deleted]]
 


__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] Maximum Likelihood Estimation

[Andre Barbosa Oliveira]

Hi,

 

 

I would like estimate a model for function of production's Coob-Douglas using 
maximum likelihood. The model is log(Y)= beta[1]+beta[2]*log(L)+beta[3]*log(K). 
I tried estimate this model using the tools nlm ( ) and optim ( ) using the 
log-likelihood function below:

 

 mloglik - function (beta, Y, L, K) {
+ n - length(Y)
+ sum ( (log(Y)- beta[1]-beta[2]*log(L)-beta[3]*log(K))^2)/2*beta[4]^2 + 
n/2*log(2*pi)+ n*log(beta[4])
+ }

 

Then I did estimates the parameters using nlm ( ) and optim ( ), but the 
estimates were very bad. I used these codes:

 

 mlem - nlm (mloglik, c(1,1,1,1), Y=Y, L=L, K=K)

 

 mlem2 - optim(c(1,1,1,1), mloglik, Y=Y, L=L, K=K, method=BFGS)

 

How I improve the estimates What's the best and more simple form for 
estimate a modelo using the  maximum likelihood's method???

 

Best regards,

 

André Barbosa Oliveira

Student of Master in Economics at University Federal of Rio Grande do Sul - 
Brazil 
  
_
Novo windowslive.com.br. Descubra como juntar a galera com os produtos Windows 
Live.


[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] Maximum likelihood estimation of parameters make no biological sense

[Luis Ridao Cruz]

R-help,

I'm trying to estimate some parameters using the Maximum Likehood method.
The model describes fish growth using a sigmoidal-type of curve:

fn_w - function(params) {
Winf   -  params[1]
k   -  params[2]
t0  -  params[3]
b   -  params[4]
   sigma -  params[5]
what  -  Winf * (1-exp(- k *(tt - t0)))^b
logL   -  -sum(dnorm(log(wobs),log(what),sqrt(sigma),TRUE))
return(logL)
}

tt  - 4:14
wobs - c(1.545, 1.920, 2.321 ,2.591, 3.676, 4.425 ,5.028, 5.877, 6.990, 6.800 
,6.900)

An then the optimization method:

OPT -optim(c(8, .1, 0, 3, 1), fn_w, method=L-BFGS-B
,lower=c(0.0, 0.001, 0.001,0.001, 0.01), upper = rep(Inf, 5), hessian=TRUE, 
control=list(trace=1))

which gives:

$par Winf   k t0 b  
sigma
[1] 24.27206813  0.04679844  0.0010  1.61760492  0.0100

$value
[1] -11.69524

$counts
function gradient 
 143  143 

$convergence
[1] 0

$message
[1] CONVERGENCE: REL_REDUCTION_OF_F = FACTR*EPSMCH

$hessian
  [,1]  [,2] [,3]  [,4]  [,5]
[1,]  1.867150e+00  1.262763e+03-7.857719 -5.153276e+01 -1.492850e-05
[2,]  1.262763e+03  8.608461e+05 -5512.469266 -3.562137e+04  9.693180e-05
[3,] -7.857719e+00 -5.512469e+0341.670222  2.473167e+02 -5.356813e+01
[4,] -5.153276e+01 -3.562137e+04   247.316675  1.535086e+03 -1.464370e-03
[5,] -1.492850e-05  9.693180e-05   -53.568127 -1.464370e-03  1.730462e+04

after iteration number 80.

From the biological point of view Winf =24(hipothesized asimptotical maximum 
weight)
makes no sense while the b parameter is no nearly close to b=3 leading to a 
non-sigmoidal
curve.

However using the least-squares method provide with more sensible parameter 
estimates

$par Winf  k  t0b
[1] 10.3827256  0.0344187  3.1751376  2.9657368

$value
[1] 2.164403

$counts
function gradient 
  48   48 

$convergence
[1] 0

$message
[1] CONVERGENCE: REL_REDUCTION_OF_F = FACTR*EPSMCH

Is there anything wrong with my MLE function and parameters?
I have tried with distinct initial parameters also.

Can anyone give me a clue on this?

Thanks in advance.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Maximum likelihood estimation of parameters make no biological sense

[dave fournier]

Hi,

Your results are do to using an unstable parameterization
of the Von Bertalanffy growth curve, combined with the unreliable
optimization methods supplied with R.  I coded up your model in
AD Model Builder which supplies exact derivatives through
AD.

I used your starting values and ran the model with no optimization steps
just to se that we had the same value for the -log-likelihood
Results are
# Number of parameters = 5  Objective function value = -11.6954  Maximum
gradien
t component = 0.0
# winf:
24.2720681300
# k:
0.046798440
# t0:
0.001000
# vhat:
0.01000
# b:
1.61760492000

However the R routine is stuck. When I let the ADMB code run it produced

# Number of parameters = 5  Objective function value = -13.8515  Maximum
gradient component = 9.41643e-05
# winf:
15.7188821203
# k:
0.118198731245
# t0:
-32.9089295327
# vhat:
0.00471832483493
# b:
184.999879271

Note that b-- infinity. I have it bounded at 185.
t0-- -infinity  so that the model is only using a small part of the
growth curve which happens to fit the data better.

The estimated correlation matrix for the parameter estimates tells the story

index name  valuestd dev   1   2   3   4   5
 1   winf 1.5719e+01 5.1252e+00   1.
 2   k1.1820e-01 2.7849e-02  -0.9832  1.
 3   t0  -3.2909e+01 7.6867e+00  -0.9748  0.9990  1.
 4   vhat 4.7183e-03 2.0119e-03   0.  0.  0.  1.
 5   b1.8500e+02 1.6374e+00  -0.0002  0.0003 -0.0094  0.  1.

You can see that several of the parameters are highly confounded.

Also the eigenvalues of the Hessian are

  0.01691149331  0.02045399106 963.2994413  2255.900979  4225373.963

So you have a condition number of about 10^8. Very difficult to work
with such a function with only approximate derivatives.

I think the moral of the story is that you should use a more stable
parameterization or an industrial strength estimation system or maybe
both.

  Cheers,

   Dave

   Cheers,

Dave







-- 
David A. Fournier
P.O. Box 2040,
Sidney, B.C. V8l 3S3
Canada
Phone/FAX 250-655-3364
http://otter-rsch.com

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] Maximum likelihood estimation of a truncated regression model

[Srinivasan Parthasarathy]

Hi,

I have a quick question regarding estimation of a truncation
regression model (truncated above at 1) using MLE in R. I will be most
grateful to you if you can help me out.

The model is linear and the relationship is dhat = bhat0+Z*bhat+e,
where dhat is the dependent variable 0 and upper truncated at 1;
bhat0 is the intercept; Z is the independent variable and is a uniform
random variate between 0 and 1; and e is the error term. I realised
that R doesn't have a built-in function to estimate truncated
regression models as does STATA, LIMDEP etc. I tried the survival and
FEAR packages and couldn't fit it for my case. So I wrote the log
likelihood function of the truncated regression model and
reparametrised it using Olsen (1978) so that the function is globally
concave and has an unique maximiser. I used a quasi-Newton method
(BFGS) to maximise my function. I also used Newton-Raphson method
(maxNR) to maximise my function. The (naive) code can be seen below.

sw1-function(theta,dhat,z)
{
gamma-theta[1:2]
eta-theta[3]
d1-dhat*eta-z%*%gamma
d2- eta-z%*%gamma
logl- (-n*0.5*log(2*pi))+(n*log(eta))-(0.5*sum(d1*d1))-(sum(log(pnorm(d2
return(-logl)
}
try(j-optim(c(a,b,c),sw1,method=BFGS,dhat=dhat,z=z),silent=TRUE)
q[i,]=abs(j$par[2]/j$par[3])

I am trying to iterate the above procedure 2000 times and compute the
RMSE of the estimated bhat. In the above code, a, b, and c are the
intercept, bhat and 1/sigma estimates from OLS of the linear model [I
get the standard error estimate from OLS and obtain the sigma by
multiplying it by sqrt(n) where n=200]. I realised that in many of the
iterations, the above code doesn't optimize as the bhat estimate I get
from the above code is worser than the one I get from STATA although
optim (or maxNR) says the function has converged successfully. I
realise that the starting parameter values could be a reason for these
sub-optimal solutions.

1. Is there an optimization routine that doesn't need starting parameter values?
2. Can you please tell me what could be the possible error in the above code?
3. Is there a routine built in R that could do truncated regression
for my situation?

Thank you very much for your help.

Best Regards,

Srini

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Maximum likelihood estimation

[Ravi Varadhan]

Hi,


You should re-write your function `fr' as follows:


fr - function(par){
a - par[1]
b - par[2]
p - par[3]
lambda - par[4]

l - 0.5*(lambda + b*p + (1-p)*(lambda-b))
l^2  lambda*b*p
delta - sqrt(abs(l^2 - b*p*lambda))
mt - a/p*(1-exp(-l*t)*cosh(delta*t)-(l-b*p)*exp(-t)*sinh(delta*t)/delta)
logl - sum(diff(y)*log(diff(mt))-diff(mt)-lfactorial(diff(y)))
return(-logl)
}

However, I don't understand what the following fragment is doing in `fr':

l^2  lambda*b*p

Can you clarfy that?

Ravi.


---

Ravi Varadhan, Ph.D.

Assistant Professor, The Center on Aging and Health

Division of Geriatric Medicine and Gerontology 

Johns Hopkins University

Ph: (410) 502-2619

Fax: (410) 614-9625

Email: [EMAIL PROTECTED]

Webpage:  http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html

 





-Original Message-
From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On
Behalf Of toh
Sent: Thursday, September 04, 2008 9:15 PM
To: r-help@r-project.org
Subject: Re: [R] Maximum likelihood estimation


Yes I'm trying to optimize the parameters a, b, p and lambda where a  0, b
 0 and 0  p  1. I attached the error message that I got when I run mle. 


 t - c(1:90)
 y -
 c(5,10,15,20,26,34,36,43,47,49,80,84,108,157,171,183,191,200,204,211,2
 17,226,230,
+
234,236,240,243,252,254,259,263,264,268,271,277,290,309,324,331,346,367,375,
381,
+
401,411,414,417,425,430,431,433,435,437,444,446,446,448,451,453,460,463,463,
464,
+
464,465,465,465,466,467,467,467,468,469,469,469,469,470,472,472,473,473,473,
473,
+ 473,473,473,475,475,475,475)
 
 library(stats4)
 fr - function(a, b, p, lambda){
+ l - 0.5*(lambda + b*p + (1-p)*(lambda-b))
+ l^2  lambda*b*p
+ delta - sqrt(abs(l^2 - b*p*lambda))
+ mt - 
+ a/p*(1-exp(-l*t)*cosh(delta*t)-(l-b*p)*exp(-t)*sinh(delta*t)/delta)
+ logl - sum(diff(y)*log(diff(mt))-diff(mt)-lfactorial(diff(y)))
+ return(-logl)
+ }
 
 mle(start=list(a=12,b=0.01,p=0.5,lambda=0.01),fr, method=L-BFGS-B,
+ lower = c(0.002, 0.002, 0.002, 0.002), upper = c(Inf, Inf, 0.999,
Inf),control=list(fnscale=-1))
Error in optim(start, f, method = method, hessian = TRUE, ...) : 
  non-finite finite-difference value [3]












Prof Brian Ripley wrote:
 
From ?optim
 
fn: A function to be minimized (or maximized), with first
argument the vector of parameters over which minimization is
to take place.  It should return a scalar result.
 
 I think you intended to optimize over c(a,b,p, lambda), so you need to 
 specify them as a single vector.
 
 You may be making life unnecessarily hard for yourself: see function 
 mle() in package stats4.
 
 Showing your code without a verbal description of what you are doing 
 nor the error message you got is less helpful than we need.
 
 On Wed, 3 Sep 2008, toh wrote:
 

 Hi R-experts,
 I'm new to R in mle. I tried to do the following but just couldn't 
 get it right. Hope someone can point out the mistakes. thanks a lot.

 t - c(1:90)
 y -
 c(5,10,15,20,26,34,36,43,47,49,80,84,108,157,171,183,191,200,204,211,
 217,226,230,

 234,236,240,243,252,254,259,263,264,268,271,277,290,309,324,331,346,3
 67,375,381,

 401,411,414,417,425,430,431,433,435,437,444,446,446,448,451,453,460,4
 63,463,464,


464,465,465,465,466,467,467,467,468,469,469,469,469,470,472,472,473,473,473,
473,
  473,473,473,475,475,475,475)
 fr - function(a, b, p, lambda){
 l - 0.5*(lambda + b*p + (1-p)*(lambda-b))
 l^2  lambda*b*p
 delta - sqrt(abs(l^2 - b*p*lambda))
 mt - 
 a/p*(1-exp(-l*t)*cosh(delta*t)-(l-b*p)*exp(-t)*sinh(delta*t)/delta)
 logl - sum(diff(y)*log(diff(mt))-diff(mt)-lfactorial(diff(y)))
 return(-logl)
 }
 optim(c(15,0.01,0.5,0.01),fr, method=L-BFGS-B, lower = c(0.002, 
 0.002, 0.002, 0.002), upper = c(Inf, Inf, 0.999,
 Inf),control=list(fnscale=-1))

 --
 View this message in context:
 http://www.nabble.com/Maximum-likelihood-estimation-tp19304249p193042
 49.html Sent from the R help mailing list archive at Nabble.com.

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

 
 -- 
 Brian D. Ripley,  [EMAIL PROTECTED]
 Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
 University of Oxford, Tel:  +44 1865 272861 (self)
 1 South Parks Road, +44 1865 272866 (PA)
 Oxford OX1 3TG, UKFax:  +44 1865 272595
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.
 
 

--
View this message

[R] Maximum likelihood estimation

[toh]

Hi R-experts,
I'm new to R in mle. I tried to do the following but just couldn't get it
right. Hope someone can point out the mistakes. thanks a lot.

t - c(1:90)
y -
c(5,10,15,20,26,34,36,43,47,49,80,84,108,157,171,183,191,200,204,211,217,226,230,

234,236,240,243,252,254,259,263,264,268,271,277,290,309,324,331,346,367,375,381,

401,411,414,417,425,430,431,433,435,437,444,446,446,448,451,453,460,463,463,464,

464,465,465,465,466,467,467,467,468,469,469,469,469,470,472,472,473,473,473,473,
473,473,473,475,475,475,475)
fr - function(a, b, p, lambda){
l - 0.5*(lambda + b*p + (1-p)*(lambda-b))
l^2  lambda*b*p
delta - sqrt(abs(l^2 - b*p*lambda))
mt - a/p*(1-exp(-l*t)*cosh(delta*t)-(l-b*p)*exp(-t)*sinh(delta*t)/delta)
logl - sum(diff(y)*log(diff(mt))-diff(mt)-lfactorial(diff(y)))
return(-logl)
}
optim(c(15,0.01,0.5,0.01),fr, method=L-BFGS-B, 
lower = c(0.002, 0.002, 0.002, 0.002), upper = c(Inf, Inf, 0.999,
Inf),control=list(fnscale=-1))

-- 
View this message in context: 
http://www.nabble.com/Maximum-likelihood-estimation-tp19304249p19304249.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Maximum likelihood estimation

[Prof Brian Ripley]

From ?optim


  fn: A function to be minimized (or maximized), with first
  argument the vector of parameters over which minimization is
  to take place.  It should return a scalar result.

I think you intended to optimize over c(a,b,p, lambda), so you need to 
specify them as a single vector.


You may be making life unnecessarily hard for yourself: see function mle() 
in package stats4.


Showing your code without a verbal description of what you are doing nor 
the error message you got is less helpful than we need.


On Wed, 3 Sep 2008, toh wrote:



Hi R-experts,
I'm new to R in mle. I tried to do the following but just couldn't get it
right. Hope someone can point out the mistakes. thanks a lot.

t - c(1:90)
y -
c(5,10,15,20,26,34,36,43,47,49,80,84,108,157,171,183,191,200,204,211,217,226,230,

234,236,240,243,252,254,259,263,264,268,271,277,290,309,324,331,346,367,375,381,

401,411,414,417,425,430,431,433,435,437,444,446,446,448,451,453,460,463,463,464,

464,465,465,465,466,467,467,467,468,469,469,469,469,470,472,472,473,473,473,473,
473,473,473,475,475,475,475)
fr - function(a, b, p, lambda){
l - 0.5*(lambda + b*p + (1-p)*(lambda-b))
l^2  lambda*b*p
delta - sqrt(abs(l^2 - b*p*lambda))
mt - a/p*(1-exp(-l*t)*cosh(delta*t)-(l-b*p)*exp(-t)*sinh(delta*t)/delta)
logl - sum(diff(y)*log(diff(mt))-diff(mt)-lfactorial(diff(y)))
return(-logl)
}
optim(c(15,0.01,0.5,0.01),fr, method=L-BFGS-B,
lower = c(0.002, 0.002, 0.002, 0.002), upper = c(Inf, Inf, 0.999,
Inf),control=list(fnscale=-1))

--
View this message in context: 
http://www.nabble.com/Maximum-likelihood-estimation-tp19304249p19304249.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.



--
Brian D. Ripley,  [EMAIL PROTECTED]
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel:  +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UKFax:  +44 1865 272595

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Maximum likelihood estimation

[toh]

Yes I'm trying to optimize the parameters a, b, p and lambda where a  0, b 
0 and 0  p  1. I attached the error message that I got when I run mle. 


 t - c(1:90)
 y -
 c(5,10,15,20,26,34,36,43,47,49,80,84,108,157,171,183,191,200,204,211,217,226,230,
+
234,236,240,243,252,254,259,263,264,268,271,277,290,309,324,331,346,367,375,381,
+
401,411,414,417,425,430,431,433,435,437,444,446,446,448,451,453,460,463,463,464,
+
464,465,465,465,466,467,467,467,468,469,469,469,469,470,472,472,473,473,473,473,
+ 473,473,473,475,475,475,475)
 
 library(stats4)
 fr - function(a, b, p, lambda){
+ l - 0.5*(lambda + b*p + (1-p)*(lambda-b))
+ l^2  lambda*b*p
+ delta - sqrt(abs(l^2 - b*p*lambda))
+ mt - a/p*(1-exp(-l*t)*cosh(delta*t)-(l-b*p)*exp(-t)*sinh(delta*t)/delta)
+ logl - sum(diff(y)*log(diff(mt))-diff(mt)-lfactorial(diff(y)))
+ return(-logl)
+ }
 
 mle(start=list(a=12,b=0.01,p=0.5,lambda=0.01),fr, method=L-BFGS-B, 
+ lower = c(0.002, 0.002, 0.002, 0.002), upper = c(Inf, Inf, 0.999,
Inf),control=list(fnscale=-1))
Error in optim(start, f, method = method, hessian = TRUE, ...) : 
  non-finite finite-difference value [3]












Prof Brian Ripley wrote:
 
From ?optim
 
fn: A function to be minimized (or maximized), with first
argument the vector of parameters over which minimization is
to take place.  It should return a scalar result.
 
 I think you intended to optimize over c(a,b,p, lambda), so you need to 
 specify them as a single vector.
 
 You may be making life unnecessarily hard for yourself: see function mle() 
 in package stats4.
 
 Showing your code without a verbal description of what you are doing nor 
 the error message you got is less helpful than we need.
 
 On Wed, 3 Sep 2008, toh wrote:
 

 Hi R-experts,
 I'm new to R in mle. I tried to do the following but just couldn't get it
 right. Hope someone can point out the mistakes. thanks a lot.

 t - c(1:90)
 y -
 c(5,10,15,20,26,34,36,43,47,49,80,84,108,157,171,183,191,200,204,211,217,226,230,

 234,236,240,243,252,254,259,263,264,268,271,277,290,309,324,331,346,367,375,381,

 401,411,414,417,425,430,431,433,435,437,444,446,446,448,451,453,460,463,463,464,

 464,465,465,465,466,467,467,467,468,469,469,469,469,470,472,472,473,473,473,473,
  473,473,473,475,475,475,475)
 fr - function(a, b, p, lambda){
 l - 0.5*(lambda + b*p + (1-p)*(lambda-b))
 l^2  lambda*b*p
 delta - sqrt(abs(l^2 - b*p*lambda))
 mt - a/p*(1-exp(-l*t)*cosh(delta*t)-(l-b*p)*exp(-t)*sinh(delta*t)/delta)
 logl - sum(diff(y)*log(diff(mt))-diff(mt)-lfactorial(diff(y)))
 return(-logl)
 }
 optim(c(15,0.01,0.5,0.01),fr, method=L-BFGS-B,
 lower = c(0.002, 0.002, 0.002, 0.002), upper = c(Inf, Inf, 0.999,
 Inf),control=list(fnscale=-1))

 -- 
 View this message in context:
 http://www.nabble.com/Maximum-likelihood-estimation-tp19304249p19304249.html
 Sent from the R help mailing list archive at Nabble.com.

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

 
 -- 
 Brian D. Ripley,  [EMAIL PROTECTED]
 Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
 University of Oxford, Tel:  +44 1865 272861 (self)
 1 South Parks Road, +44 1865 272866 (PA)
 Oxford OX1 3TG, UKFax:  +44 1865 272595
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.
 
 

-- 
View this message in context: 
http://www.nabble.com/Maximum-likelihood-estimation-tp19304249p19323140.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Maximum likelihood estimation

[Peter Dalgaard]

Jurica Brajković wrote:
 Hello, 

 I am struggling for some time now to estimate AR(1) process for commodity 
 price time series. I did it in STATA but cannot get a result in R.

 The equation I want to estimate is: p(t)=a+b*p(t-1)+error 
 Using STATA I get 0.92 for a, and 0.73 for b. 

 Code that I use in R is: 
 p-matrix(data$p) # price at time t
 lp-cbind(1,data$lp) # price at time t-1

  
 mle - function(theta) {
   sigma2-theta[1]
   b- theta[-1]
   n-length(p)
   e-p-lp%*%b
   logl- -(n/2)*log(sigma2)-((t(e)%*%e)/(2*sigma2))
   return(-logl)
   }


 out - optim(c(0,0,0),mle, method = L-BFGS-B,
 lower = c(0, -Inf, -Inf),
 upper = c(Inf, Inf, Inf))

 The result I get is:  Error in optim(c(0, 0, 0), mle, method = L-BFGS-B, 
 lower = c(0, -Inf,:L-BFGS-B needs finite values of 'fn'

 Can somebody spot the mistake?

 Many thanks, 

 Jurica Brajkovic 

   
As far as I can see, the first element of the vector of starting values
(sigma2) is 0 and you are taking log of it

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.
   


-- 
   O__   Peter Dalgaard Øster Farimagsgade 5, Entr.B
  c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K
 (*) \(*) -- University of Copenhagen   Denmark  Ph:  (+45) 35327918
~~ - ([EMAIL PROTECTED])  FAX: (+45) 35327907

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] Maximum likelihood estimation

[Jurica Brajković]

Hello, 

I am struggling for some time now to estimate AR(1) process for commodity price 
time series. I did it in STATA but cannot get a result in R.

The equation I want to estimate is: p(t)=a+b*p(t-1)+error 
Using STATA I get 0.92 for a, and 0.73 for b. 

Code that I use in R is: 
p-matrix(data$p) # price at time t
lp-cbind(1,data$lp) # price at time t-1

 
mle - function(theta) {
  sigma2-theta[1]
  b- theta[-1]
  n-length(p)
  e-p-lp%*%b
  logl- -(n/2)*log(sigma2)-((t(e)%*%e)/(2*sigma2))
  return(-logl)
  }


out - optim(c(0,0,0),mle, method = L-BFGS-B,
lower = c(0, -Inf, -Inf),
upper = c(Inf, Inf, Inf))

The result I get is:  Error in optim(c(0, 0, 0), mle, method = L-BFGS-B, 
lower = c(0, -Inf,:L-BFGS-B needs finite values of 'fn'

Can somebody spot the mistake?

Many thanks, 

Jurica Brajkovic 

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] Maximum Likelihood Estimation

[Todd Brauer]

Using R, I would like to calculate algorithms to estimate coefficients á and â 
within the gamma function: f(costij)=((costij)^á)*exp(â*costij).  I have its 
logarithmic diminishing line data (Logarithmic Diminishing Line Data Table) and 
have installed R¢s Maximum Likelihood Estimation package; however, I am unsure 
which method to apply in order to calculate the algorisms (i.e., Newton-Raphson 
Maximization, Nelder-Mead Maximization, etc.)  Any guidance you all could 
provide would be appreciated.
 
Logarithmic Diminishing Line Data Table
1  0.983385666
2  0.578408021
3  0.421101402
4  0.334555838
5  0.278826347
6  0.239521701
7  0.210100828
8  0.187133124
9  0.168633944
10   0.15336976
11   0.140530624
12   0.129561113
13   0.120066673
14   0.111758725
15   0.104420943
16   0.097887721
17   0.092029955
18   0.086745393
19   0.081951911
20   0.077582726
21   0.073582928
22   0.069906916
23   0.066516481
24   0.06337934
25   0.060468017
26   0.057758964
27   0.055231875
28   0.052869138
29   0.050655392
30   0.048577177
31   0.046622642
32   0.044781306
33   0.043043865
34   0.041402028
35   0.03984838
36   0.038376265
37   0.036979694
38   0.03565326
39   0.034392065
40   0.033191667
41   0.03204802
42   0.030957435
43   0.029916538
44   0.02892224
45   0.027971702
46   0.027062313
47   0.026191669
48   0.025357547
49   0.024557893
50   0.023790804
51   0.023054514
52   0.022347383
53   0.021667883
54   0.021014591
55   0.02038618
56   0.01978141
57   0.019199121
58   0.018638226
59   0.018097707
60   0.017576607
61   0.017074029
62   0.016589127
63   0.016121106
64   0.015669216
65   0.015232751
66   0.014811043
67   0.014403462
68   0.014009411
69   0.013628326
70   0.013259673
71   0.012902944
72   0.01255766
73   0.012223364
74   0.011899623
75   0.011586025
76   0.011282179
77   0.010987712
78   0.01070227
79   0.010425513
80   0.01015712
81   0.009896783
82   0.009644208
83   0.009399116
84   0.009161239
85   0.008930321
86   0.008706116
87   0.008488392
88   0.008276923
89   0.008071496
90   0.007871904
91   0.00767795
92   0.007489446
93   0.007306209
94   0.007128067
95   0.006954853
96   0.006786404
97   0.006622569
98   0.006463198
99   0.00630815
100 0.006157287
 
Thanks,
Todd


  
[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Maximum Likelihood Estimation

[Ben Bolker]

Todd Brauer toddbrauer at yahoo.com writes:

 
 Using R, I would like to calculate algorithms to estimate coefficients á and â
within the gamma function:
 f(costij)=((costij)^á)*exp(â*costij).  I have its logarithmic diminishing line
data
 (Logarithmic Diminishing Line Data Table) and have installed R¢s Maximum
Likelihood Estimation
 package; however, I am unsure which method to apply in order to calculate the
algorisms (i.e.,
 Newton-Raphson Maximization, Nelder-Mead Maximization, etc.)  Any guidance you
all could provide
 would be appreciated.

  For a simple, well-behaved problem (like this one) it doesn't
really matter that much what algorithm you pick. optim()'s default
is Nelder-Mead (typically slightly more robust and less efficient).
Here's some sample code to do a least-squares fit:

x - read.table(gammafit.txt)
names(x) - c(x,y)
plot(x)
sum(x)
gfun - function(p) {
a - p[1]
b - p[2]
exp.y - x$x^b*exp(-a*x$x)
sum((exp.y-x$y)^2)
}
opt1 - optim(gfun,par=c(a=0.1,b=1))
v - opt1$par
curve(x^v[b]*exp(-v[a]*x),add=TRUE,col=2)

  However, if you really want to do MLE you may need to specify
your problem a little more carefully -- sum of squares will
find a reasonable answer but likelihood-based inference won't
be right.  Is this a survival problem (in which case you
might want to use the survival package)?  If you construct
your own likelihood function, you might need
to use pgamma(...,lower.tail=FALSE).

  Ben Bolker

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Maximum likelihood estimation in R with censored Data

[Thomas Lumley]

Try survreg(), in the survival package.

 -thomas

On Fri, 13 Jun 2008, Bluder Olivia wrote:


Hello,



I'm trying to calculate the Maximum likelihood estimators for a dataset
which contains censored data.



I started by using the function nlm, but isn't there a separate method
for doing this for e.g. the weibull and the log-normal distribution?



Thanks,

Olivia




[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.



Thomas Lumley   Assoc. Professor, Biostatistics
[EMAIL PROTECTED]   University of Washington, Seattle

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] Maximum likelihood estimation in R with censored Data

[Bluder Olivia]

Hello,

 

I'm trying to calculate the Maximum likelihood estimators for a dataset
which contains censored data.

 

I started by using the function nlm, but isn't there a separate method
for doing this for e.g. the weibull and the log-normal distribution?

 

Thanks,

Olivia

 


[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Maximum likelihood estimation in R with censored Data

[Ben Bolker]

Bluder Olivia olivia.bluder at k-ai.at writes:

 
 Hello,
 
 I'm trying to calculate the Maximum likelihood estimators for a dataset
 which contains censored data.
 
 I started by using the function nlm, but isn't there a separate method
 for doing this for e.g. the weibull and the log-normal distribution?
 
 Thanks,
 
 Olivia

  This is not *quite* enough detail about what you
want to do.  Can you (as the posting guide suggests!)
give us a small example of what you want to do?  You may be able
to do this via the survreg() command in the survival
package, or you may want to do it yourself by constructing
a log-likelihood function with dweibull() for uncensored
data and pweibull() for censored data [or dlnorm/plnorm].

  Ben Bolker

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Maximum likelihood estimation in R with censored Data

[Vincent Goulet]

Le ven. 13 juin à 13:55, Ben Bolker a écrit :


Bluder Olivia olivia.bluder at k-ai.at writes:



Hello,

I'm trying to calculate the Maximum likelihood estimators for a  
dataset

which contains censored data.

I started by using the function nlm, but isn't there a separate  
method
for doing this for e.g. the weibull and the log-normal  
distribution?


Thanks,

Olivia


 This is not *quite* enough detail about what you
want to do.  Can you (as the posting guide suggests!)
give us a small example of what you want to do?  You may be able
to do this via the survreg() command in the survival
package, or you may want to do it yourself by constructing
a log-likelihood function with dweibull() for uncensored
data and pweibull() for censored data [or dlnorm/plnorm].


If you want to go the second route, function coverage() in package  
actuar will build the censored density function for you. You can then  
feed this function to fitdistr() just like for usual ML estimation.


HTH  Vincent




 Ben Bolker

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.