[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 at 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