Maximum likelihood

From an old post to the R mailing list by Peter Dalgaard:

“the possibility of finding maximum likelihood estimates by writing up the likelihood and maximizing it is often overlooked…”

R really is magical sometimes. Suppose you want to fit a distribution, M. All you need is to maximise \prod_{i=1}^N P(x_i | M), or equivalently, \sum_{i=1}^N \log P(x_i | M). Here’s an example of fitting a Gaussian, starting by breaking a fairly good first guess…

> x = rnorm(1000, 100, 15)
> f = function(p) -2*sum(dnorm(x, p[1], p[2], log=T))
> optim(c(mean(x)-50, sd(x)+15), f)
$par
[1] 100 15

$value
[1] 8193

$counts
function gradient
69 NA

$convergence
[1] 0

$message
NULL

(Well actually -2 times the log-likelihood.) Now to have a look at your estimate:

hist(x, probability=T)
curve(dnorm(x,100,15), min(x), max(x), add=T)

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s