最大似然估计实例 | Fitting a Model by Maximum Likelihood (MLE)
参考：Fitting a Model by Maximum Likelihood
最大似然估计是用于估计模型参数的，首先我们必须选定一个模型，然后比对有给定的数据集，然后构建一个联合概率函数，因为给定了数据集，所以该函数就是以模型参数为自变量的函数，通过求导我们就能得到使得该函数值（似然值）最大的模型参数了。
Maximum-Likelihood Estimation (MLE) is a statistical technique for estimating model parameters. It basically sets out to answer the question: what model parameters are most likely to characterise a given set of data? First you need to select a model for the data. And the model must have one or more (unknown) parameters. As the name implies, MLE proceeds to maximise a likelihood function, which in turn maximises the agreement between the model and the data.
Most illustrative examples of MLE aim to derive the parameters for a probability density function (PDF) of a particular distribution. In this case the likelihood function is obtained by considering the PDF not as a function of the sample variable, but as a function of distribution’s parameters. For each data point one then has a function of the distribution’s parameters. The joint likelihood of the full data set is the product of these functions. This product is generally very small indeed, so the likelihood function is normally replaced by a log-likelihood function. Maximising either the likelihood or log-likelihood function yields the same results, but the latter is just a little more tractable!
Fitting a Normal Distribution
Let’s illustrate with a simple example: fitting a normal distribution. First we generate some data.
> set.seed(1001)>> N <- 100>> x <- rnorm(N, mean = 3, sd = 2)>> mean(x)[1] 2.998305> sd(x)[1] 2.288979 |
Then we formulate the log-likelihood function.
> LL <- function(mu, sigma) {+ R = dnorm(x, mu, sigma)+ #+ -sum(log(R))+ } |
And apply MLE to estimate the two parameters (mean and standard deviation) for which the normal distribution best describes the data.
> library(stats4)>> mle(LL, start = list(mu = 1, sigma=1))Call:mle(minuslogl = LL, start = list(mu = 1, sigma = 1))Coefficients: mu sigma2.998305 2.277506Warning messages:1: In dnorm(x, mu, sigma) : NaNs produced2: In dnorm(x, mu, sigma) : NaNs produced3: In dnorm(x, mu, sigma) : NaNs produced |
Those warnings are a little disconcerting! They are produced when negative values are attempted for the standard deviation.
> dnorm(x, 1, -1) [1] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN [30] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN [59] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN [88] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN |
There are two ways to sort this out. The first is to apply constraints on the parameters. The mean does not require a constraint but we insist that the standard deviation is positive.
> mle(LL, start = list(mu = 1, sigma=1), method = "L-BFGS-B", lower = c(-Inf, 0), upper = c(Inf, Inf))Call:mle(minuslogl = LL, start = list(mu = 1, sigma = 1), method = "L-BFGS-B", lower = c(-Inf, 0), upper = c(Inf, Inf))Coefficients: mu sigma2.998304 2.277506 |
This works because mle() calls optim(), which has a number of optimisation methods. The default method is BFGS. An alternative, the L-BFGS-B method, allows box constraints.
The other solution is to simply ignore the warnings. It’s neater and produces the same results.
> LL <- function(mu, sigma) {+ R = suppressWarnings(dnorm(x, mu, sigma))+ #+ -sum(log(R))+ }>> mle(LL, start = list(mu = 1, sigma=1))Call:mle(minuslogl = LL, start = list(mu = 1, sigma = 1))Coefficients: mu sigma2.998305 2.277506 |
The maximum-likelihood values for the mean and standard deviation are damn close to the corresponding sample statistics for the data. Of course, they do not agree perfectly with the values used when we generated the data: the results can only be as good as the data. If there were more samples then the results would be closer to these ideal values.
A note of caution: if your initial guess for the parameters is too far off then things can go seriously wrong!
> mle(LL, start = list(mu = 0, sigma=1))Call:mle(minuslogl = LL, start = list(mu = 0, sigma = 1))Coefficients: mu sigma 51.4840 226.8299 |
Fitting a Linear Model
Now let’s try something a little more sophisticated: fitting a linear model. As before, we generate some data.
> x <- runif(N)> y <- 5 * x + 3 + rnorm(N) |
We can immediately fit this model using least squares regression.
> fit <- lm(y ~ x)>> summary(fit)Call:lm(formula = y ~ x)Residuals: Min 1Q Median 3Q Max-1.96206 -0.59016 -0.00166 0.51813 2.43778Coefficients: Estimate Std. Error t value Pr(>|t|)(Intercept) 3.1080 0.1695 18.34 <2e-16 ***x 4.9516 0.2962 16.72 <2e-16 ***---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 0.8871 on 98 degrees of freedomMultiple R-squared: 0.7404, Adjusted R-squared: 0.7378F-statistic: 279.5 on 1 and 98 DF, p-value: < 2.2e-16 |
The values for the slope and intercept are very satisfactory. No arguments there. We can superimpose the fitted line onto a scatter plot.
> plot(x, y)> abline(fit, col = "red") |
Pushing on to the MLE for the linear model parameters. First we need a likelihood function. The model is not a PDF, so we can’t proceed in precisely the same way that we did with the normal distribution. However, if you fit a linear model then you want the residuals to be normally distributed. So the likelihood function fits a normal distribution to the residuals.
LL <- function(beta0, beta1, mu, sigma) { # Find residuals # R = y - x * beta1 - beta0 # # Calculate the likelihood for the residuals (with mu and sigma as parameters) # R = suppressWarnings(dnorm(R, mu, sigma)) # # Sum the log likelihoods for all of the data points # -sum(log(R))} |
One small refinement that one might make is to move the logarithm into the call to dnorm().
LL <- function(beta0, beta1, mu, sigma) { R = y - x * beta1 - beta0 # R = suppressWarnings(dnorm(R, mu, sigma, log = TRUE)) # -sum(R)} |
It turns out that the initial guess is again rather important and a poor choice can result in errors. We will return to this issue a little later.
> fit <- mle(LL, start = list(beta0 = 3, beta1 = 1, mu = 0, sigma=1))Error in solve.default(oout$hessian) : system is computationally singular: reciprocal condition number = 3.01825e-22> fit <- mle(LL, start = list(beta0 = 5, beta1 = 3, mu = 0, sigma=1))Error in solve.default(oout$hessian) : Lapack routine dgesv: system is exactly singular: U[4,4] = 0 |
But if we choose values that are reasonably close then we get a decent outcome.
> fit <- mle(LL, start = list(beta0 = 4, beta1 = 2, mu = 0, sigma=1))> fitCall:mle(minuslogl = LL, start = list(beta0 = 4, beta1 = 2, mu = 0, sigma = 1))Coefficients: beta0 beta1 mu sigma 3.5540217 4.9516133 -0.4459783 0.8782272 |
The maximum-likelihood estimates for the slope (beta1) and intercept (beta0) are not too bad. But there is a troubling warning about NANs being produced in the summary output below.
> summary(fit)Maximum likelihood estimationCall:mle(minuslogl = LL, start = list(beta0 = 4, beta1 = 2, mu = 0, sigma = 1))Coefficients: Estimate Std. Errorbeta0 3.5540217 NaNbeta1 4.9516133 0.2931924mu -0.4459783 NaNsigma 0.8782272 0.0620997-2 log L: 257.8177Warning message:In sqrt(diag(object@vcov)) : NaNs produced |
It stands to reason that we actually want to have the zero mean for the residuals. We can apply this constraint by specifying mu as a fixed parameter. Another option would be to simply replace mu with 0 in the call to dnorm(), but the alternative is just a little more flexible.
> fit <- mle(LL, start = list(beta0 = 2, beta1 = 1.5, sigma=1), fixed = list(mu = 0), nobs = length(y))> summary(fit)Maximum likelihood estimationCall:mle(minuslogl = LL, start = list(beta0 = 2, beta1 = 1.5, sigma = 1), fixed = list(mu = 0), nobs = length(y))Coefficients: Estimate Std. Errorbeta0 3.1080423 0.16779428beta1 4.9516164 0.29319233sigma 0.8782272 0.06209969-2 log L: 257.8177 |
The resulting estimates for the slope and intercept are rather good. And we have standard errors for these parameters as well.
How about assessing the overall quality of the model? We can look at the Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC). These can be used to compare the performance of different models for a given set of data.
> AIC(fit)[1] 263.8177> BIC(fit)[1] 271.6332> logLik(fit)'log Lik.' -128.9088 (df=3) |
Returning now to the errors mentioned above. Both of the cases where the call to mle() failed resulted from problems with inverting the Hessian Matrix. With the implementation of mle() in the stats4 package there is really no way to get around this problem apart from having a good initial guess. In some situations though, this is just not feasible. There are, however, alternative implementations of MLE which circumvent this problem. The bbmle package has mle2() which offers essentially the same functionality but includes the option of not inverting the Hessian Matrix.
> library(bbmle)> > fit <- mle2(LL, start = list(beta0 = 3, beta1 = 1, mu = 0, sigma = 1))> > summary(fit)Maximum likelihood estimationCall:mle2(minuslogl = LL, start = list(beta0 = 3, beta1 = 1, mu = 0, sigma = 1))Coefficients: Estimate Std. Error z value Pr(z) beta0 3.054021 0.083897 36.4019 <2e-16 ***beta1 4.951617 0.293193 16.8886 <2e-16 ***mu 0.054021 0.083897 0.6439 0.5196 sigma 0.878228 0.062100 14.1421 <2e-16 ***---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1-2 log L: 257.8177 |
Here mle2() is called with the same initial guess that broke mle(), but it works fine. The summary information for the optimal set of parameters is also more extensive.
Fitting a linear model is just a toy example. However, Maximum-Likelihood Estimation can be applied to models of arbitrary complexity. If the model residuals are expected to be normally distributed then a log-likelihood function based on the one above can be used. If the residuals conform to a different distribution then the appropriate density function should be used instead of dnorm().
