Let \(\mathcal{F}\) be a parametric statistical model for the data \(y\), where \(f(y;\theta)\) is the density or probability function and \(\theta \in \Theta\) is a p-dimensional parameter. By considering \(f(y;\theta)\) as function of \(\theta\) with \(y\) fixed to the observed value, then the likelihood function of \(\theta\) based on \(y\), \(L: \Theta \rightarrow \mathbb{R}^{+}\), is \[L(\theta) = L(\theta; y) \propto f_Y(y; \theta) \,\,, \] where the proportionality is related to a constant \(c(y)>0\) that does not depend on \(\theta\).

On the basis of the data \(y\), \(\theta \in \Theta\) is more likely than \(\theta' \in \Theta\) as an index of the data generating model if \(L(\theta) > L(\theta')\)

• Let \(y=(y_{1},\ldots, y_{n})\) be a sample of i.i.d. values from a Bernoulli distribution, \(Y \sim \mbox{Be}(p)\), with \(n = 100\) and \(p = 0.6\). The log-likelihood function is \[\ell(p) = \sum_{i=1}^{n} y_i \log (p) + \Big(n - \sum_{i=1}^{n} y_i\Big) \log (1-p)\] and the MLE is available in closed form, \(\hat p = \bar y\): this allows us to check the output of the numerical optimisers.
• Logit reparametrisation: \(\psi = {\rm logit}(p) = \log\{p/(1-p)\}\) maps \((0, 1)\) into \(\mathbb{R}\), with inverse \(p(\psi) = \exp(\psi)/\{1+\exp(\psi)\}\). By the invariance of the MLE, \(\hat \psi = {\rm logit}(\hat p)\).

In the Binomial example, the MLE is available in closed form and we obtain numerical values just by plugging the data into the formula. However, remember that the MLE for a parameter may not have an analytic solution. Many times we do not have a closed form for MLE estimates and we may need numerical optimisation. R provides various functions for performing numerical methods:

• nlm(): minimizes a function using a Newton-type algorithm. It needs a starting value and does not allow constraints on the parameters. It is usually fast and reliable. It returns the ML estimate \(\hat{\theta}\) (estimate), the minimum \(-\ell(\hat{\theta})\) of the negative log-likelihood (minimum) and the hessian (hessian), if hessian = TRUE.

• optim(): minimizes a function using Nelder-Mead, quasi-Newton and conjugate gradients algorithms. It includes an option for box-constrained optimization, and it requires a starting value. It returns the ML estimate \(\hat{\theta}\) (par), the minimum of the negative log-likelihood \(-\ell(\hat{\theta})\) (value) and the hessian (hessian), if hessian = TRUE. You also can maximize the function by using fnscale = -1 in control argument.

• nlminb(): often is more stable, robust and reliable than optim (in particular with “nasty” functions). It performs only minimization and does not yield numerical derivatives as output. It returns the ML estimate \(\hat{\theta}\) (par) and the minimum of the negative log-likelihood \(-\ell(\hat{\theta})\) (objective).

• optimize(): by using a combination of golden section search and successive parabolic interpolation, searches in an interval for a minimum or a maximum (if maximum=TRUE) of a function. It returns \(\hat{\theta}\) (maximum if maximum=TRUE, otherwise minimum) and the value of the function there (objective). Drawback: suited only for one-dimensional parameter.

Here, we are simply using the Binomial example to show how (some of) these functions work.

Note: In the following, we consider and implement the negative log-likelihood function, since some numerical optimisers that we will see are only able to perform minimisation of functions.

In addition, we implement the score and the observed information (used for comparison with the optimiser output)

To check if the score and the observed information are correctly implemented, we can compute the numerical gradient and hessian of the log-likelihood function, as implemented in the {numDeriv} R package.

Let’s start with the first optimiser nlm (considering as starting value 0.5, which is quite close to the MLE; however, see the warnings printed in the R console)

The previous optimiser does not allow constraints on the parameters; however, there are some replacements in the algorithm that allow you to reach the result. If you want to optimise a function with constrained parameters, it is better to use an optimiser allowing for box-constrained optimisation, as for instance optim.

• The optim() function provides a lot of numerical methods (Nelder-Mead, quasi-Newton, conjugate-gradient methods and simulated annealing). As a drawback, the user has to set up carefully the initial guesses and the method, since the final solution may be quite sensitive to these choices…To compute the observed information, the function optimHess() computes numerical derivatives of generic functions.

• As an alternative, we could consider again the reparametrisation introduced above. Generally, for numerical purposes, it is convenient to reparameterize the model in such a way that the new parameter space is unbounded. Then, we explore this by using nlm()

• Waiting times (minutes) between consecutive orders, \(Y_i \stackrel{iid}\sim\) Exponential with mean \(\theta\): \(f(y; \theta) = \frac{1}{\theta} e^{-y/\theta}, \ y > 0\) (R uses the rate: rexp(n, rate = 1/theta))
• Log-likelihood, up to an additive constant: \(\ell(\theta) = -n \log \theta - S/\theta\), with \(S = \sum_{i=1}^n y_i\)
• Data: \(n = 8\) waiting times (see below); PART 1, 2 and 3 of the R script are completed in teams

• Compute ll_exp for the two samples at \(\theta = 3, 5, 8\). What do you notice?

• True \(\theta = 5\), samples of size \(n = 8\), \(R = 10000\) replications (each row of Y is one sample)
• 2a. Evaluate the score and the observed information at the true value, in each sample

• By plugging in: \(\hat \lambda = 1/\bar y\). Is it unbiased? (true \(\lambda = 0.2\))

• LRT: \(W(\theta_0)\); Wald: \((\hat \theta - \theta_0)^2 J(\hat \theta)\); Score: \(U(\theta_0)^2/\mathcal{I}(\theta_0)\), with \(\mathcal{I}(\theta) = n/\theta^2\)
• Under \(H_0\), the three statistics are approximately \(\chi^2_1\)