1.2.1. Computing the direction d_k

A desirable property for d_k is that d_k ensures a descent f(x_k + 1) < f(x_k). Newton methods are such that d_k minimizes a local quadratic approximation of f based on a Taylor expansion, that is $q_f(d) = f(x_k) + g(x_k)^Td +\frac{1}{2} d^T H(x_k) d$ where g denotes the gradient and H denotes the Hessian.

The consists in using the exact solution of local minimization problem d_k = −H(x_k)⁻¹g(x_k).
In practice, other methods are preferred (at least to ensure positive definiteness). The method approximates the Hessian by a matrix H_k as a function of H_k − 1, x_k, f(x_k) and then d_k solves the system H_kd = −g(x_k). Some implementation may also directly approximate the inverse of the Hessian W_k in order to compute d_k = −W_kg(x_k). Using the Sherman-Morrison-Woodbury formula, we can switch between W_k and H_k.

To determine W_k, first it must verify the secant equation H_ky_k = s_k or y_k = W_ks_k where y_k = g_k + 1 − g_k and s_k = x_k + 1 − x_k. To define the n(n − 1) terms, we generally impose a symmetry and a minimum distance conditions. We say we have a rank 2 update if H_k = H_k − 1 + auu^T + bvv^T and a rank 1 update if $H_k = H_{k-1} + a u u^T $. Rank n update is justified by the spectral decomposition theorem.

There are two rank-2 updates which are symmetric and preserve positive definiteness

• DFP minimizes min ||H − H_k||_F such that H = H^T: $$ H_{k+1} = \left (I-\frac {y_k s_k^T} {y_k^T s_k} \right ) H_k \left (I-\frac {s_k y_k^T} {y_k^T s_k} \right )+\frac{y_k y_k^T} {y_k^T s_k} \Leftrightarrow W_{k+1} = W_k + \frac{s_k s_k^T}{y_k^{T} s_k} - \frac {W_k y_k y_k^T W_k^T} {y_k^T W_k y_k} . $$
  
• BFGS minimizes min ||W − W_k||_F such that W = W^T: $$ H_{k+1} = H_k - \frac{ H_k y_k y_k^T H_k }{ y_k^T H_k y_k } + \frac{ s_k s_k^T }{ y_k^T s_k } \Leftrightarrow W_{k+1} = \left (I-\frac {y_k s_k^T} {y_k^T s_k} \right )^T W_k \left (I-\frac { y_k s_k^T} {y_k^T s_k} \right )+\frac{s_k s_k^T} {y_k^T s_k} . $$

In R, the so-called BFGS scheme is implemented in optim.

Another possible method (which is initially arised from quadratic problems) is the nonlinear conjugate gradients. This consists in computing directions (d₀, …, d_k) that are conjugate with respect to a matrix close to the true Hessian H(x_k). Directions are computed iteratively by d_k = −g(x_k) + β_kd_k − 1 for k > 1, once initiated by d₁ = −g(x₁). β_k are updated according a scheme:

• $\beta_k = \frac{ g_k^T g_k}{g_{k-1}^T g_{k-1} }$: Fletcher-Reeves update,
• $\beta_k = \frac{ g_k^T (g_k-g_{k-1} )}{g_{k-1}^T g_{k-1}}$: Polak-Ribiere update.

There exists also three-term formula for computing direction d_k = −g(x_k) + β_kd_k − 1 + γ_kd_t for t < k. A possible scheme is the Beale-Sorenson update defined as $\beta_k = \frac{ g_k^T (g_k-g_{k-1} )}{d^T_{k-1}(g_{k}- g_{k-1})}$ and $\gamma_k = \frac{ g_k^T (g_{t+1}-g_{t} )}{d^T_{t}(g_{t+1}- g_{t})}$ if k > t + 1 otherwise γ_k = 0 if k = t. See Yuan (2006) for other well-known schemes such as Hestenses-Stiefel, Dixon or Conjugate-Descent. The three updates (Fletcher-Reeves, Polak-Ribiere, Beale-Sorenson) of the (non-linear) conjugate gradient are available in optim.

1.2.2. Computing the stepsize t_k

Let ϕ_k(t) = f(x_k + td_k) for a given direction/iterate (d_k, x_k). We need to find conditions to find a satisfactory stepsize t_k. In literature, we consider the descent condition: ϕ_k′(0) < 0 and the Armijo condition: ϕ_k(t) ≤ ϕ_k(0) + tc₁ϕ_k′(0) ensures a decrease of f. Nocedal & Wright (2006) presents a backtracking (or geometric) approach satisfying the Armijo condition and minimal condition, i.e. Goldstein and Price condition.

• set t_k, 0 e.g. 1, 0 < α < 1,
• Repeat until Armijo satisfied,
  • t_{k, i + 1} = α × t_k, i.
• end Repeat

This backtracking linesearch is available in optim.

2.1.1. Theoretical value

The density of the beta distribution is given by $$ f(x; \delta_1,\delta_2) = \frac{x^{\delta_1-1}(1-x)^{\delta_2-1}}{\beta(\delta_1,\delta_2)}, $$ where β denotes the beta function, see the NIST Handbook of mathematical functions https://dlmf.nist.gov/. We recall that β(a, b) = Γ(a)Γ(b)/Γ(a + b). There the log-likelihood for a set of observations (x₁, …, x_n) is $$ \log L(\delta_1,\delta_2) = (\delta_1-1)\sum_{i=1}^n\log(x_i)+ (\delta_2-1)\sum_{i=1}^n\log(1-x_i)+ n \log(\beta(\delta_1,\delta_2)) $$ The gradient with respect to a and b is $$ \nabla \log L(\delta_1,\delta_2) = \left(\begin{matrix} \sum\limits_{i=1}^n\ln(x_i) - n\psi(\delta_1)+n\psi( \delta_1+\delta_2) \\ \sum\limits_{i=1}^n\ln(1-x_i)- n\psi(\delta_2)+n\psi( \delta_1+\delta_2) \end{matrix}\right), $$ where ψ(x) = Γ′(x)/Γ(x) is the digamma function, see the NIST Handbook of mathematical functions https://dlmf.nist.gov/.

2.1.2. R implementation

As in the fitdistrplus package, we minimize the opposite of the log-likelihood: we implement the opposite of the gradient in grlnL. Both the log-likelihood and its gradient are not exported.

lnL <- function(par, fix.arg, obs, ddistnam) 
  fitdistrplus:::loglikelihood(par, fix.arg, obs, ddistnam) 
grlnlbeta <- fitdistrplus:::grlnlbeta

3.1.1. Theoretical value

The p.m.f. of the Negative binomial distribution is given by $$ f(x; m,p) = \frac{\Gamma(x+m)}{\Gamma(m)x!} p^m (1-p)^x, $$ where Γ denotes the beta function, see the NIST Handbook of mathematical functions https://dlmf.nist.gov/. There exists an alternative representation where μ = m(1 − p)/p or equivalently p = m/(m + μ). Thus, the log-likelihood for a set of observations (x₁, …, x_n) is $$ \log L(m,p) = \sum_{i=1}^{n} \log\Gamma(x_i+m) -n\log\Gamma(m) -\sum_{i=1}^{n} \log(x_i!) + mn\log(p) +\sum_{i=1}^{n} {x_i}\log(1-p) $$ The gradient with respect to m and p is $$ \nabla \log L(m,p) = \left(\begin{matrix} \sum_{i=1}^{n} \psi(x_i+m) -n \psi(m) + n\log(p) \\ mn/p -\sum_{i=1}^{n} {x_i}/(1-p) \end{matrix}\right), $$ where ψ(x) = Γ′(x)/Γ(x) is the digamma function, see the NIST Handbook of mathematical functions https://dlmf.nist.gov/.

3.1.2. R implementation

As in the fitdistrplus package, we minimize the opposite of the log-likelihood: we implement the opposite of the gradient in grlnL.

grlnlNB <- function(x, obs, ...)
{
  m <- x[1]
  p <- x[2]
  n <- length(obs)
  c(sum(psigamma(obs+m)) - n*psigamma(m) + n*log(p),
    m*n/p - sum(obs)/(1-p))
}