STA 506 - Midterm Examination

Midterm Exam Objectives

• Understand the definition and relationship between PDFs and CDFs, including their non-parametric estimators: the empirical distribution function and kernel density estimation (KDE).

• Estimate sampling distributions using simulation-based methods, specifically the bootstrap.

• Derive point estimates of parameters using the method of moments and maximum likelihood estimation (MLE).

• Describe the asymptotic (normal) and bootstrap sampling distributions of maximum likelihood estimators.

• Apply all the above inferential procedures in a programming environment to perform numerical data analysis.

Policies of Using AI Tools

• Policy on AI Tool Use: Students must adhere to the AI tool policy specified in the course syllabus. The direct copying of AI-generated content is strictly prohibited. All submitted work must reflect your own understanding; where external tools are consulted, content must be thoroughly rephrased and synthesized in your own words.

• Code Inclusion Requirement: Any code included in your essay must be properly commented to explain the purpose and/or expected output of key code lines. Submitting AI-generated code without meaningful, student-added comments will not be accepted.

Policies of the Midterm Exam

• This exam must be completed independently. Collaboration with others is not permitted.

• You may consult class notes and online resources. However, please ensure that your use of AI tools complies with the AI policy established for this course.

• All work must be submitted by the deadline. Submissions (including partial work) received after the deadline will not be accepted.

Submission Preparation

Since this exam involves derivations, you can use the following approaches to prepare your submission and save time:

1. Handwrite the derivation on paper and take a photo. From there, you can either use AI tools to convert it into LaTeX code and paste it into RMarkdown, or insert the image directly into your RMarkdown document, following the example shown in class this week.

2. Prepare the formulas using other non-programmatic editors and take a screenshot to embed the image in your RMarkdown document using knitr::include_graphics.

3. Compile all handwritten derivations into a single file and submit it as an attachment.

Please follow the same submission procedure as used for the weekly assignments.

Exam Questions

Caution: Please follow the suggested expressions and guided steps to complete the exam. Other approaches, such as transformations that trivialize the problems, will not meet the exam objectives.

Log-normal Distribution

The log-normal distributions have been widely used in many different fields such as finance and economics, engineering and reliability, environmental science, medical and biology etc.

If \(X\) follows a lognormal distribution:

\[ f(x; \mu, \sigma^2) = \frac{1}{x\sqrt{2\pi \sigma^2}} \exp\left[-\frac{(\ln x - \mu)^2}{2\sigma^2}\right], \quad x > 0 \]

where \(\mu \in \mathbb{R}\) is the mean of \(\ln X\) and \(\sigma^2 > 0\) is the variance of \(\ln X\).

Caution: \(\mu\) and \(\sigma^2\) are not population mean and variance of the lognormal distribution!. Instead, we will use \(\mu_{LN}\) and \(\sigma_{LN}^2\) to denote the lognormal mean and variance respectively.

For an i.i.d. random sample \(x_1, x_2, \dots, x_n\):

\[\begin{align*} \ell(\mu, \sigma^2) &= \sum_{i=1}^n \ln f(x_i; \mu, \sigma^2) \\ &= -\frac{n}{2} \ln(2\pi) - \frac{n}{2} \ln \sigma^2 - \sum_{i=1}^n \ln x_i - \frac{1}{2\sigma^2} \sum_{i=1}^n (\ln x_i - \mu)^2 \end{align*}\]

The \(k\)-th moment of a lognormal distribution given by

\[ \boxed{ \mathbb{E}[X^k] = \exp\left( k\mu + \frac{1}{2} k^2 \sigma^2 \right). } \]

Particularly, the mean of the above lognormal distribution, denoted by \(\mu_{LN}\), is given by

\[ \mu_{LN} = \mathbb{E}[X] = \exp\left( \mu + \frac{1}{2} \sigma^2 \right). \]

Once the parameters were estimated, say \((\overline{\mu}, \overline{\sigma^2})\), we can use the plug-in Principle to estimate the lognormal mean by

\[ \boxed{\overline{\mu_{LN}} = \exp\left( \overline{\mu} + \frac{1}{2} \overline{\sigma^2} \right).} \]

Working Dataset: A reservoir engineer at an oil exploration company has recently drilled 50 new wells in a shale formation. According to geological theory, the size of oil reserves in this type of formation is often modeled using a log-normal distribution. The engineer is responsible for estimating the typical reserve size and quantifying the uncertainty of these estimates for a presentation to senior management. The estimated ultimate recovery (EUR) for the 50 wells, measured in thousands of barrels (Mbbl), is as follows:

17.3, 24.8, 8.2, 31.5, 14.1, 42.7, 11.9, 55.3, 21.4, 9.7, 36.2, 18.6, 63.1, 13.2, 28.9, 
47.5, 19.8, 7.5, 33.4, 52.1, 16.0, 25.9, 38.7, 10.3, 44.2, 22.5, 58.6, 12.8, 30.3, 48.9, 
20.1, 35.4, 15.7, 60.2, 26.7, 41.3, 18.1, 53.8, 23.9, 46.2, 29.4, 37.6, 14.9, 50.5, 32.8, 
19.3, 56.7, 11.2, 39.5, 27.1

This assignment focuses on finding the asymptotic sampling distributions of the MLE \(\hat{\mu}\) and \(\hat{\sigma^2}\) corresponding to \(\mu\) and \(\sigma^2\) based on an environmental study data.

Question 1: Method of Moments Estimation (MME) of \(\mu\) and \(\sigma^2\)

Estimate the parameters \(\mu\) and \(\sigma^2\) (caution: not \(\sigma\)!) from the raw oil reserve data using the given \(k\)-th moment function. Organize your work into three parts:

1. Derive the system of two equations for the method of moments estimators (MMEs) of \(\mu\) and \(\sigma^2\).

The \(k\)-th moment of lognormal random variable is given by \[ \boxed{ \mathbb{E}[X^k] = \exp\left( k\mu + \frac{1}{2} k^2 \sigma^2 \right). } \]

Since there are two unknown parameters, \(\mu\) and \(\sigma^2\), we use the first two moments.

Use the first moment to create the first equation

For \(k=1\), \[ \mathbb{E}[X] = \exp\left(\mu + \frac{\sigma^2}{2}\right) \]

Set the population moment equal to the first sample moment.

\[ \frac{1}{n}\sum_{i=1}^{n} X_i = \exp\left(\mu + \frac{\sigma^2}{2}\right) \]

Use the second moment to create the second equation

For \(k=2\), \[ \mathbb{E}[X^2] = \exp\left(2\mu + 2\sigma^2\right) \]

Set the population second moment equal to the second sample moment.

\[ \frac{1}{n}\sum_{i=1}^{n} X_i^2 = \exp\left(2\mu + 2\sigma^2\right) \]

These two equations form the method of moments system

\[ \frac{1}{n}\sum_{i=1}^{n} X_i = \exp\left(\mu + \frac{\sigma^2}{2}\right) \]

\[ \frac{1}{n}\sum_{i=1}^{n} X_i^2 = \exp\left(2\mu + 2\sigma^2\right) \]

2. Derive the closed-form expressions for the MMEs of \(\mu\), \(\sigma^2\), and \(\mu_{LN}\), denoted by \(\widetilde{\mu}\), \(\widetilde{\sigma^2}\), and \(\widetilde{\mu_{LN}}\), respectively.

We now solve the system from part (a) to obtain closed-form MMEs.

Let

\[ m_1 = \frac{1}{n}\sum_{i=1}^{n}X_i \]

\[ m_2 = \frac{1}{n}\sum_{i=1}^{n}X_i^2 \]

These represent the first and second sample moments.

Take the natural log of both equations to remove the exponential terms.

\[ \ln(m_1) = \mu + \frac{\sigma^2}{2} \]

\[ \ln(m_2) = 2\mu + 2\sigma^2 \]

Multiply the first equation by 2.

\[ 2\ln(m_1) = 2\mu + \sigma^2 \]

Subtract the equations to isolate \(\sigma^2\).

\[ \widetilde{\sigma^2} = \ln(m_2) - 2\ln(m_1) \]

Substitute this result back into the first equation to solve for \(\mu\).

\[ \widetilde{\mu} = \ln(m_1) - \frac{1}{2}\widetilde{\sigma^2} \]

Using the plug-in principle, the estimator for the lognormal mean is

\[ \widetilde{\mu_{LN}} = \exp\left(\widetilde{\mu} + \frac{1}{2}\widetilde{\sigma^2}\right) \]

3. Write an R function that outputs the MMEs \(\widetilde{\mu}\), \(\widetilde{\sigma^2}\), and \(\widetilde{\mu_{LN}}\) based on the given data set.

# Compute the first two sample moments
mme_lognormal <- function(x) {
  m1 <- mean(x)
  m2 <- mean(x^2)
  
  # Compute MME of sigma^2
  sigma2_tilde <- log(m2) - 2 * log(m1)
  
  # Compute MME of mu
  mu_tilde <- log(m1) - 0.5 * sigma2_tilde
  
  # Compute plug-in estimate of lognormal mean
  mu_LN_tilde <- exp(mu_tilde + 0.5 * sigma2_tilde)
  
  # Return all three estimates
  return(list(
    mu_tilde = mu_tilde,
    sigma2_tilde = sigma2_tilde,
    mu_LN_tilde = mu_LN_tilde
  ))
}

# Store the data
oil_data <- c(
  17.3, 24.8, 8.2, 31.5, 14.1, 42.7, 11.9, 55.3, 21.4, 9.7,
  36.2, 18.6, 63.1, 13.2, 28.9, 47.5, 19.8, 7.5, 33.4, 52.1,
  16.0, 25.9, 38.7, 10.3, 44.2, 22.5, 58.6, 12.8, 30.3, 48.9,
  20.1, 35.4, 15.7, 60.2, 26.7, 41.3, 18.1, 53.8, 23.9, 46.2,
  29.4, 37.6, 14.9, 50.5, 32.8, 19.3, 56.7, 11.2, 39.5, 27.1
)

# Compute the estimates
mme_lognormal(oil_data)

$mu_tilde
[1] 3.301269

$sigma2_tilde
[1] 0.2339652

$mu_LN_tilde
[1] 30.516

Question 2: Finding the MLE of \(\mu\) and \(\sigma^2\)

This question involves two parts:

(1). Derive the score functions (i.e., the gradient of the log-likelihood) for the parameters of the log-normal distribution, and set them to zero to obtain the system of score equations. Begin with the full log-likelihood and compute its partial derivatives with respect to \(\mu\) and \(\sigma^2\). [Hint: letting \(\beta = \sigma^2\) may simplify differentiation.]

Assume that \(\{x_1, x_2, \cdots, x_n \} \to \text{ LN }(\mu, \sigma^2)\) with density function given by

\[ \ell(\mu, \sigma^2) = -\frac{n}{2} \ln(2\pi) - n \ln \sigma - \sum_{i=1}^n \ln x_i - \frac{1}{2\sigma^2} \sum_{i=1}^n (\ln x_i - \mu)^2 \]

Let \(\beta = \sigma^2\)

\[ \ell(\mu,\beta) = -\frac{n}{2}\ln(2\pi) -\frac{n}{2}\ln\beta -\sum_{i=1}^n \ln x_i -\frac{1}{2\beta}\sum_{i=1}^n(\ln x_i-\mu)^2 \]

Score function for \(\mu\)

\[ \frac{\partial \ell}{\partial \mu} = -\frac{1}{2\beta} \frac{\partial}{\partial \mu} \left[ \sum_{i=1}^n(\ln x_i-\mu)^2 \right] \]

\[ = -\frac{1}{2\beta} \sum_{i=1}^n 2(\ln x_i-\mu)(-1) \]

\[ = \frac{1}{\beta}\sum_{i=1}^n(\ln x_i-\mu) \]

Set score equation equal to zero

\[ \frac{1}{\beta}\sum_{i=1}^n(\ln x_i-\mu)=0 \]

\[ \sum_{i=1}^n(\ln x_i-\mu)=0 \]

Score function for \(\beta\)

\[ \frac{\partial \ell}{\partial \beta} = -\frac{n}{2\beta} - \frac{\partial}{\partial \beta} \left[ \frac{1}{2\beta}\sum_{i=1}^n(\ln x_i-\mu)^2 \right] \]

\[ = -\frac{n}{2\beta} + \frac{1}{2\beta^2}\sum_{i=1}^n(\ln x_i-\mu)^2 \]

Set score equation equal to zero

\[ -\frac{n}{2\beta} + \frac{1}{2\beta^2}\sum_{i=1}^n(\ln x_i-\mu)^2 =0 \]

Score equations

\[ \sum_{i=1}^n \ln x_i - n\mu = 0 \]

\[ \sum_{i=1}^n(\ln x_i-\mu)^2 - n\sigma^2 = 0 \] (2). Perform some algebra to find the closed form of the MLE of \(\mu\) and \(\sigma^2\), denoted by \(\widehat{\mu}\) and \(\widehat{\sigma^2}\).

Solve the score equations

\[ \sum_{i=1}^{n} \ln x_i - n\mu = 0 \]

Solve for \(\mu\):

\[ \sum_{i=1}^{n} \ln x_i = n\mu \]

\[ \widehat{\mu} = \frac{1}{n}\sum_{i=1}^{n}\ln x_i \]

Now use the second score equation:

\[ \sum_{i=1}^{n}(\ln x_i-\mu)^2 - n\sigma^2 = 0 \]

Substitute \(\widehat{\mu}=\frac{1}{n}\sum_{i=1}^{n}\ln x_i\) into the equation:

\[ \sum_{i=1}^{n}(\ln x_i-\widehat{\mu})^2 - n\sigma^2 = 0 \]

Solve for \(\sigma^2\):

\[ \sum_{i=1}^{n}(\ln x_i-\widehat{\mu})^2 = n\sigma^2 \]

\[ \widehat{\sigma^2} = \frac{1}{n}\sum_{i=1}^{n}(\ln x_i-\widehat{\mu})^2 \]

(3). Recall that lognormal population mean (i.e., the first moment) is

\[ \mu_{LN} = \mathbb{E}[X] = \exp\left( \mu + \frac{1}{2} \sigma^2 \right). \]

Using the relationship between variance and the first and second moments, \(\text{Var}(X) = \mathbb{E}[X^2] - \left(\mathbb{E}[X]\right)^2\), we derive the population variance of the lognormal distribution, denoted as \(\sigma_{\text{LN}}^2\). This result will be used in applying the Central Limit Theorem in subsequent questions

Lognormal population mean

\[ \mu_{LN} = E[X] = \exp\left(\mu + \frac{1}{2}\sigma^2\right) \]

Second moment of a lognormal variable

\[ E[X^2] = \exp(2\mu + 2\sigma^2) \]

Use the variance identity

\[ \operatorname{Var}(X) = E[X^2] - (E[X])^2 \]

\[ \sigma_{LN}^2 = \exp(2\mu + 2\sigma^2) - \left[\exp\left(\mu + \frac{1}{2}\sigma^2\right)\right]^2 \]

Simplify the second term

\[ \left[\exp\left(\mu + \frac{1}{2}\sigma^2\right)\right]^2 = \exp(2\mu + \sigma^2) \]

Factor common term

\[ \sigma_{LN}^2 = \exp(2\mu + 2\sigma^2) - \exp(2\mu + \sigma^2) \]

\[ \sigma_{LN}^2 = \exp(2\mu + \sigma^2)\left(\exp(\sigma^2)-1\right) \]

Question 3: Sampling Distribution of Lognormal Sample Mean \(\widehat{\mu_{LN}}\)

We have derived two different estimators of lognormal population mean \(\mu_{LN}\). This Question focuses on the sampling distribution of lognormal sample sample means based on MLE. The following are detailed instructions.

(1). Write an R function to estimate the MLE of \(\mu\) and \(\sigma^2\), denoted by \(\widehat{\mu}\) and \(\widehat{\sigma^2}\) and return the MLE of \(\mu_{LN}\), denoted by \(\widehat{\mu_{LN}}\).

Compute MLEs from one sample

# Compute MLEs from one sample
mle_lognormal <- function(x) {

  log_x <- log(x)
  n <- length(x)

  # MLE of mu (derived in Question 2)
  mu_hat <- mean(log_x)

  # MLE of sigma^2 (derived in Question 2)
  sigma2_hat <- sum((log_x - mu_hat)^2) / n

  # Plug-in estimate of lognormal mean
  mu_LN_hat <- exp(mu_hat + 0.5 * sigma2_hat)

  return(list(
    mu_hat = mu_hat,
    sigma2_hat = sigma2_hat,
    mu_LN_hat = mu_LN_hat
  ))
}

Apply the function to the oil reserve data

mle_lognormal(oil_data)

$mu_hat
[1] 3.268536

$sigma2_hat
[1] 0.3271989

$mu_LN_hat
[1] 30.94264

(2). Take 1000 Bootstrap samples from the raw oil data set. For each bootstrap sample, call the above R function to find the bootstrap MLE of \(\mu_{LN}\), denoted by \(\widehat{\mu_{LN}^*}^{(1)}, \widehat{\mu_{LN}^*}^{(2)}, \cdots, \widehat{\mu_{LN}^*}^{(1000)}\).

# Number of bootstrap samples
B <- 1000

# Sample size
n <- length(oil_data)

# Store bootstrap MLEs of mu_LN
boot_mu_LN <- numeric(B)

for(i in 1:B) {

  # Draw bootstrap sample with replacement
  boot_sample <- sample(oil_data, size = n, replace = TRUE)

  # Compute MLEs from the bootstrap sample
  boot_fit <- mle_lognormal(boot_sample)

  # Store bootstrap MLE of mu_LN
  boot_mu_LN[i] <- boot_fit$mu_LN_hat
}

#Verify that 1000 bootstrap estimates were generated
length(boot_mu_LN)

[1] 1000

(3). Using kernel density estimation approach to estimate the bootstrap density curves based on the bootstrap MLEs: \(\widehat{\mu_{LN}^*}^{(1)}, \widehat{\mu_{LN}^*}^{(2)}, \cdots, \widehat{\mu_{LN}^*}^{(1000)}\).

Using the bootstrap estimates \(\widehat{\mu}_{LN}^{*(1)}, \widehat{\mu}_{LN}^{*(2)}, \ldots, \widehat{\mu}_{LN}^{*(1000)}\) obtained in part (2), we estimate the bootstrap sampling density using kernel density estimation.

# Estimate the bootstrap density using kernel density estimation
boot_density <- density(boot_mu_LN)

# Plot the bootstrap kernel density estimate
plot(boot_density,
     main = "Bootstrap Density of Lognormal Mean",
     xlab = expression(hat(mu[LN])),
     ylab = "Density",
     lwd = 2)

(4). In order to find the asymptotic sampling distribution using the central limit theorem (CLT). That is, when sample size is large,

By the central limit theorem, when the sample size is large,

\[ \bar{X} \approx N\left(\mu_{LN}, \frac{\sigma_{LN}^2}{n}\right). \]

From Question 2, the plug-in estimate of the lognormal mean is

\[ \widehat{\mu}_{LN} = \exp\left(\widehat{\mu} + \frac{1}{2}\widehat{\sigma}^2\right). \]

From the question page, the lognormal population variance is

\[ \sigma_{LN}^2 = \mathbb{E}[X^2] - (\mathbb{E}[X])^2. \]

Since

\[ \mathbb{E}[X] = \exp\left(\mu + \frac{1}{2}\sigma^2\right) \quad\text{and}\quad \mathbb{E}[X^2] = \exp(2\mu + 2\sigma^2), \]

the lognormal variance is

\[ \sigma_{LN}^2 = \exp(2\mu + 2\sigma^2) - \exp(2\mu + \sigma^2). \]

Applying the MLE plug-in principle,

\[ \widehat{\sigma}_{LN}^2 = \exp(2\widehat{\mu} + 2\widehat{\sigma}^2) - \exp(2\widehat{\mu} + \widehat{\sigma}^2). \]

The asymptotic sampling distribution is therefore estimated by

\[ \bar{X} \approx N\left(\widehat{\mu}_{LN}, \frac{\widehat{\sigma}_{LN}^2}{n}\right). \]

(5). Write a short summary to compare the two sampling distributions of the sample mean of the lognormal distribution.

The bootstrap sampling distribution of the lognormal sample mean was estimated using 1000 bootstrap MLEs of \(\mu_{LN}\). The resulting kernel density curve is centered near the estimated lognormal mean \(\widehat{\mu}_{LN} \approx 30.94\). The asymptotic sampling distribution obtained from the central limit theorem is a normal distribution with mean \(\widehat{\mu}_{LN}\) and variance \(\widehat{\sigma}_{LN}^2/n\). In the plot, the asymptotic normal density closely follows the bootstrap density in both location and spread. This indicates that the normal approximation provides a reasonable description of the sampling distribution of the lognormal sample mean for this sample size.

Question 4: Extra Bonus Credit

Derive the Hessian matrix by computing the second-order partial derivatives of the log-likelihood function (with respect to \(\mu\) and \(\sigma^2\)). Then, translate the derived Hessian into an R function that takes the MLEs of the parameters and the data as inputs and returns the Hessian matrix. Finally, compare your analytically derived Hessian with the numerical Hessian returned by optim().

(1) Analytical Hessian Matrix

The Hessian matrix is the matrix of second derivatives of the log-likelihood function.

\[ H(\mu,\sigma^2) = \begin{bmatrix} \frac{\partial^2 \ell}{\partial \mu^2} & \frac{\partial^2 \ell}{\partial \mu \partial \sigma^2} \\ \frac{\partial^2 \ell}{\partial \sigma^2 \partial \mu} & \frac{\partial^2 \ell}{\partial (\sigma^2)^2} \end{bmatrix} \]

From Question 2,

\[ \frac{\partial \ell}{\partial \mu} = \frac{1}{\sigma^2}\sum_{i=1}^{n}(\ln x_i-\mu) \]

Taking another derivative,

\[ \frac{\partial^2 \ell}{\partial \mu^2} = -\frac{n}{\sigma^2} \]

The mixed derivative is

\[ \frac{\partial^2 \ell}{\partial \mu \partial \sigma^2} = -\frac{1}{(\sigma^2)^2}\sum_{i=1}^{n}(\ln x_i-\mu) \]

From Question 2,

\[ \frac{\partial \ell}{\partial \sigma^2} = -\frac{n}{2\sigma^2} + \frac{1}{2(\sigma^2)^2} \sum_{i=1}^{n}(\ln x_i-\mu)^2 \]

Taking another derivative,

\[ \frac{\partial^2 \ell}{\partial (\sigma^2)^2} = \frac{n}{2(\sigma^2)^2} - \frac{1}{(\sigma^2)^3} \sum_{i=1}^{n}(\ln x_i-\mu)^2 \]

Therefore the Hessian matrix is

\[ H(\mu,\sigma^2) = \begin{bmatrix} -\frac{n}{\sigma^2} & -\frac{1}{(\sigma^2)^2}\sum_{i=1}^{n}(\ln x_i-\mu) \\ -\frac{1}{(\sigma^2)^2}\sum_{i=1}^{n}(\ln x_i-\mu) & \frac{n}{2(\sigma^2)^2} - \frac{1}{(\sigma^2)^3}\sum_{i=1}^{n}(\ln x_i-\mu)^2 \end{bmatrix} \]

(2) Analytical Hessian Matrix in R

# Analytical Hessian matrix based on formulas from part (1)
hessian_lognormal <- function(mu, sigma2, x){

  log_x <- log(x)
  n <- length(x)

  H11 <- -n / sigma2
  H12 <- -(1 / (sigma2^2)) * sum(log_x - mu)
  H21 <- H12
  H22 <- (n / (2 * sigma2^2)) - (1 / (sigma2^3)) * sum((log_x - mu)^2)

  matrix(c(H11, H12,
           H21, H22),
         nrow = 2,
         byrow = TRUE)
}

# Compute Hessian using the MLE estimates
mle_fit <- mle_lognormal(oil_data)

H_analytical <- hessian_lognormal(
  mle_fit$mu_hat,
  mle_fit$sigma2_hat,
  oil_data
)

H_analytical

              [,1]          [,2]
[1,] -1.528123e+02  1.037019e-13
[2,]  1.037019e-13 -2.335158e+02

(3) Numerical Hessian from optim()

# Negative log-likelihood function
neg_loglik <- function(par, x){

  mu <- par[1]
  sigma2 <- par[2]

  log_x <- log(x)
  n <- length(x)

  -( -n/2*log(2*pi)
     -n/2*log(sigma2)
     -sum(log_x)
     -(1/(2*sigma2))*sum((log_x - mu)^2) )
}

# Compute numerical Hessian
optim_fit <- optim(
  par = c(mle_fit$mu_hat, mle_fit$sigma2_hat),
  fn = neg_loglik,
  x = oil_data,
  hessian = TRUE
)

H_numerical <- optim_fit$hessian

H_numerical

         [,1]     [,2]
[1,] 152.8123   0.0000
[2,]   0.0000 233.5289

The analytical Hessian matrix and the numerical Hessian obtained from optim() have the same magnitudes but opposite signs because optim() computes the Hessian of the negative log-likelihood. This confirms that the analytical Hessian derivation is correct.-