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.

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.