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 log normal 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\).

Answer to Question 1 Part A:

We can use the the formula for the \(k\)-th moment of the log normal distribution to determine the 1st and 2nd population moments:

$$ \[\begin{aligned} \mu_1 &= \mathbb{E}[X] = \exp \left(\mu+\frac{1}{2}\sigma^2 \right) \\\\ \mu_2 &= \mathbb{E}[X^2]= \exp \left(2\mu+2\sigma^2\right) \end{aligned}\]

$$

We can then determine the corresponding sample moments:

$$ \[\begin{aligned} m_1 &= \frac{1}{n}\sum_{i=1}^{n}{X_i} \\\\ m_2 &= \frac{1}{n}\sum_{i=1}^{n}{X_i^2} \end{aligned}\]

$$

Therefore the system of two equations for the method of moments estimator of \(\mu\) and \(\sigma^2\) are:

$$ \[\begin{aligned} \exp \left(\mu+\frac{1}{2}\sigma^2 \right) &= m_1\\\\ \exp \left(2\mu+2\sigma^2\right) &= m_2 \end{aligned}\]

$$

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.

Answer to Question 1 Part B:

We can solve the above system of equations to derive the closed-form expressions for the MMEs of \(\mu\), \(\sigma^2\) and \(\mu_{LN}\), denoted by \(\tilde{\mu}\), \(\tilde{\sigma^2}\), and \(\tilde{\mu_{LN}}\):

First we will solve both equations for \(2\mu\) so that we can solve for \(\sigma^2\):

$$ \[\begin{aligned} \exp \left(\mu+\frac{1}{2}\sigma^2 \right) &= m_1 \Rightarrow \mu+\frac{1}{2}\sigma^2=\ln{m_1} \Rightarrow 2\mu=2\ln{m_1}-\sigma^2\\\\ \exp \left(2\mu+2\sigma^2\right) &= m_2 \Rightarrow2\mu+2\sigma^2=\ln{m_2} \Rightarrow 2\mu=\ln{m_2}-2\sigma^2 \end{aligned}\]

$$

We can then set the equations equal to each other and solve for \(\sigma^2\):

\[ 2\ln{m_1}-\sigma^2=\ln{m_2}-2\sigma^2 \Rightarrow\sigma^2=\ln{m_2}-2\ln{m_1} \]

Having determined \(\sigma^2\) we can plug this into one of our previous equations to determine \(\mu\):

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

Finally, we can plug our values for \(\mu\) and \(\sigma^2\) into the equation:

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

To get:

\[ \mu_{LN}=\exp\left(2\ln{m_1}-\frac{1}{2}\ln{m_2}+\frac{1}{2}(\ln{m_2}-2\ln{m_1})\right)=\exp(\ln{m_1})=m_1 \]

Therefore, the closed-form expressions for the MMEs of \(\mu\), \(\sigma^2\) and \(\mu_{LN}\), denoted by \(\tilde{\mu}\), \(\tilde{\sigma^2}\), and \(\tilde{\mu_{LN}}\) are:

$$ \[\begin{aligned} \tilde{\mu} &= 2\ln{m_1}-\frac{1}{2}\ln{m_2} \\\\ \tilde{\sigma^2} &= \ln{m_2}-2\ln{m_1} \\\\ \tilde{\mu_{LN}} &= m_1 \end{aligned}\]

$$

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

Answer to Question 1 Part C:

oil = 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) #Entering dataset

momfunction <- function(data){ #Defining function use later
m1 <- mean(data)
m2 <- mean(data**2) #Used to get first and second moment

mom_mu <- 2*log(m1)-(1/2)*log(m2) #Used to get mu-tilde

mom_sigma_sq <- log(m2)-2*log(m1) #Used to get sigma_squared-tilde

mom_mu_ln <- m1 #Used to get mu_LN-tilde

return(list(mu_tilde = mom_mu, sigma_squared_tilde = mom_sigma_sq, mu_ln_tilde = mom_mu_ln)) #Returns values
}

momfunction(oil) #Using function on dataset

$mu_tilde
[1] 3.301269

$sigma_squared_tilde
[1] 0.2339652

$mu_ln_tilde
[1] 30.516

Based on our output we can determine that for our given data set \(\tilde{\mu}=3.301269\), \(\tilde{\sigma^2}=0.2339652\), and \(\tilde{\mu}_{LN}=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 \]

Answer to Question 2 Part 1:

First we will take the partial derivative of the log-likelihood function with respect to \(\mu\):

$$ \[\begin{aligned} \frac{\partial\ell}{\partial\mu} &= \frac{\partial}{\partial\mu}\left[-\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 \right]\\\\ &= \frac{1}{\sigma^2}\sum_{i=1}^{n}{(\ln x_i - \mu)} \end{aligned}\]

$$

Then we will take the partial derivative with respect to \(\sigma^2\):

$$ \[\begin{aligned} \frac{\partial\ell}{\partial\sigma^2} &= \frac{\partial}{\partial\sigma^2}\left[-\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 \right]\\\\ &= \frac{\partial}{\partial\sigma^2}\left[-\frac{n}{2} \ln(2\pi) - n \ln (\sigma^2)^{1/2} - \sum_{i=1}^n \ln x_i - \frac{1}{2\sigma^2} \sum_{i=1}^n (\ln x_i - \mu)^2 \right] \\\\ &= -\frac{n}{2\sigma^2}+\frac{1}{2\sigma^4}\sum_{i=1}^n (\ln x_i - \mu)^2 \end{aligned}\]

$$

Finally we will set the derivatives equal to zero to obtain the system of score equations:

$$ \[\begin{aligned} \frac{\partial\ell}{\partial\mu} &= \frac{1}{\sigma^2}\sum_{i=1}^{n}{(\ln x_i - \mu)} = 0\\\\ \frac{\partial\ell}{\partial\sigma^2} &= -\frac{n}{2\sigma^2}+\frac{1}{2\sigma^4}\sum_{i=1}^n (\ln x_i - \mu)^2 = 0 \end{aligned}\]

$$

(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}\).

#Answer to Question 2 Part 2:

First we will solve the system of score equations for the \(\mu\) to obtain \(\hat{\mu}\):

$$ \[\begin{aligned} \frac{\partial\ell}{\partial\mu}=0 &\Rightarrow \frac{1}{\sigma^2}\sum_{i=1}^{n}{(\ln x_i - \mu)} = 0\\\\ &\Rightarrow \sum_{i=1}^{n}{(\ln x_i - \mu)} = 0 \\\\ &\Rightarrow \sum_{i=1}^{n}{(\ln x_i)}-n\mu=0 \\\\ &\Rightarrow \frac{1}{n}\sum_{i=1}^{n}{(\ln x_i)}=\mu \end{aligned}\]

$$

Then we will solve for \(\sigma^2\) to obtain \(\hat{\sigma^2}\):

$$ \[\begin{aligned} \frac{\partial\ell}{\partial\sigma^2}=0 &\Rightarrow -\frac{n}{2\sigma^2}+\frac{1}{2\sigma^4}\sum_{i=1}^n (\ln x_i - \mu)^2 = 0\\\\ &\Rightarrow -n\sigma^2+\sum_{i=1}^n (\ln x_i - \mu)^2=0\\\\ &\Rightarrow \frac{1}{n}\sum_{i=1}^n (\ln x_i - \mu)^2=\sigma^2 \end{aligned}\]

$$

Finally we can plug in our previously determined formula for \(\mu\):

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

Therefore, the closed form of the MLE of \(\mu\) and \(\sigma^2\), denoted by \(\widehat{\mu}\) and \(\widehat{\sigma^2}\) is:

$$ \[\begin{aligned} \hat{\mu} &= \frac{1}{n}\sum_{i=1}^{n}{(\ln x_i)}\\\\ \hat{\sigma^2} &= \frac{1}{n}\sum_{i=1}^n \left(\ln x_i - \frac{1}{n}\sum_{i=1}^{n}{\ln x_i}\right)^2 \end{aligned}\]

$$

(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

Answer to Question 2 Part 3:

Since we know \(\mathbb{E}[X] = \exp\left( \mu + \frac{1}{2} \sigma^2 \right)\) and we previously determined that \(\mathbb{E}[X^2]= \exp \left(2\mu+2\sigma^2\right)\) we can plug \(\mathbb{E}[X]\) and \(\mathbb{E}[X^2]\) into \(\text{Var}(X) = \mathbb{E}[X^2] - \left(\mathbb{E}[X]\right)^2\):

$$ \[\begin{aligned} \sigma_{\text{LN}}^2=\text{Var}(X) &= \mathbb{E}[X^2] - \left(\mathbb{E}[X]\right)^2 \\\\ &= \exp \left(2\mu+2\sigma^2\right)-\left[\exp\left( \mu + \frac{1}{2} \sigma^2 \right)\right]^2 \\\\ &=\exp \left(2\mu+2\sigma^2\right)-\exp\left(2\mu+\sigma^2 \right) \end{aligned}\]

$$

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}}\).

Answer to Question 3 Part 1:

mlefunction <- function(data){ #Defining function used later

n <- length(data) #Used to get n
  
mle_mu <- mean(log(data)) #Used to get mu-hat

mle_sigma_sq <- (1/n)*sum((log(data)-mle_mu)^2) #Used to get sigma^2-hat

mle_mu_ln <- exp(mle_mu+(1/2)*mle_sigma_sq) #Used to get mu_LN-hat

return(list(mu_hat = mle_mu, sigma_squared_hat = mle_sigma_sq, mu_ln_hat = mle_mu_ln)) #Used to return values
}

mlefunction(oil) #Using function on datset

$mu_hat
[1] 3.268536

$sigma_squared_hat
[1] 0.3271989

$mu_ln_hat
[1] 30.94264

Based on our output we can determine that for our given data set \(\hat{\mu}=3.268536\), \(\hat{\sigma^2}=0.3271989\), and \(\hat{\mu}_{LN}=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)}\).

Answer to Question 3 Part 2:

n <- length(oil) #Sets sample size for bootstrap equal to length of original sample
B <- 1000 #How many times we should bootstrap sample

bootstrap.mu_ln_hat <- numeric(B) #Resets 'bootstrap.mu_ln_hat' 

for (i in 1:B){
  boot.sample <- sample(oil, size=n, replace=TRUE)
  result <- mlefunction(boot.sample) #Uses mlefunction created above and stores results in 'result'
  bootstrap.mu_ln_hat[i] <- result$mu_ln_hat
  #Used to randomly sample (with replacement) from 'oil' 'B' times to get sample mu_ln_hats
}

(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)}\).

hist(bootstrap.mu_ln_hat, prob=TRUE, #Creates histogram of bootstrap mu_ln_hats
     main = "Bootstrap Distribution of Mu_LN-Hats", #Title
     xlab = "Mu_LN-Hats") #x-label
lines(density(bootstrap.mu_ln_hat, bw=1), col="steelblue", lwd=2) #Creates density curve using kernel density estimation 

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

\[ \bar{X} \to N \left(\mu_{LN}, \frac{\sigma_{LN}}{\sqrt{n}}\right) \]

where \(\mu_{LN}\) and \(\sigma_{LN}\) can be estimated by using the plug-in principle of MLE. Based on the asymptotic sampling distribution of lognormal sample mean \(\bar{X}\), add the asymptotic density curve to the above Bootstrap kernel density curve.

Answer to Question 3 Part 4:

set.seed(123)

mleoutput <- mlefunction(oil) #Used to define function used later

mu <- mleoutput$mu_hat #Retrieves previously found mu-hat
  
sigma_sqr <- mleoutput$sigma_squared_hat #Retrieves previously found sigma^2-hat

n <- length(oil) #Used to get n

asympt_mu_ln <- exp(mu+(1/2)*sigma_sqr) #Calculates asympt_mu_ln
  
asympt_sigma_ln <- sqrt(exp(2*mu+2*sigma_sqr)-exp(2*mu+sigma_sqr))/(sqrt(n)) #Calculates asympt_sigma_ln

x <- seq(min(bootstrap.mu_ln_hat), max(bootstrap.mu_ln_hat), length=1000) #Used for normal curve later

hist(bootstrap.mu_ln_hat, prob=TRUE, #Creates histogram of bootstrap mu_ln_hats
     main = "Asymptotic Density Curve vs. Bootstrap Kernel Density Curve",
     xlab = "Mu_LN-Hats")
lines(density(bootstrap.mu_ln_hat, bw=1), col="steelblue", lwd=2) #Creates density curve using kernel density estimation approach
lines(x, dnorm(x, mean=asympt_mu_ln, sd=asympt_sigma_ln), #Creates normal curve using aympt_mu_ln and asympt_sigma_ln
      col="red",
      lwd=2)

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

Answer to Question 3 Part 5:

Both the bootstrap distribution (in blue) and the asymptotic sampling distribution (in red) are fairly symmetric as shown by the figure above. Both distributions appear to roughly center around 30, which makes sense given that our estimate for \(\hat{\mu}_{LN}\) was 30.94264. The asymptotic sampling distribution seems to have a slightly wider spread (given that it has a lower peak and faintly heavier tails). However, overall the asymptotic sampling distribution closely follows the bootstrap distributions indicating that the CLT provides a fairly accurate approximation of the sampling distribution.

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().

We previously found that the partial derivatives of the log likelihood function are:

$$ \[\begin{aligned} \frac{\partial\ell}{\partial\mu} &= \frac{1}{\sigma^2}\sum_{i=1}^{n}{(\ln x_i - \mu)} \\\\ \frac{\partial\ell}{\partial\sigma^2} &= -\frac{n}{2\sigma^2}+\frac{1}{2\sigma^4}\sum_{i=1}^n (\ln x_i - \mu)^2 \end{aligned}\]

$$

We can use these as the starting point to calculate the 2nd derivatives used for the Hessian matrix:

\[ \mathcal{H}(\hat{\mu}, \hat{\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} \]

To fill in the Hessian matrix we will calculate and fill in the partial derivatives:

$$ \[\begin{aligned} \frac{\partial^2\ell}{\partial\mu^2} &= \frac{\partial^2}{\partial\mu^2}\left[ \frac{1}{\hat{\sigma^2}}\sum_{i=1}^{n}{(\ln x_i - \hat{\mu})}\right]= -\frac{n}{\hat{\sigma^2}}\\\\ \frac{\partial^2\ell}{\partial\mu\partial\sigma^2} &= \frac{\partial^2}{\partial\mu\partial\sigma^2}\left[\frac{1}{\hat{\sigma^2}}\sum_{i=1}^{n}{(\ln x_i - \hat{\mu})} \right]=-\frac{1}{\hat{\sigma^4}}\sum_{i=1}^{n}{(\ln x_i - \hat{\mu})}\\\\ \frac{\partial^2\ell}{\partial(\sigma^2)^2} &= \frac{\partial^2}{\partial(\sigma^2)^2}\left[-\frac{n}{2\hat{\sigma^2}}+\frac{1}{2\hat{\sigma^4}}\sum_{i=1}^n (\ln x_i - \hat{\mu})^2 \right]=\frac{n}{2\hat{\sigma^4}}-\frac{1}{\hat{\sigma^6}}\sum_{i=1}^n (\ln x_i - \hat{\mu})^2 \end{aligned}\]

$$

Therefore the Hessian Matrix is:

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

We can use our derived Hessian matrix formula to create a function that can be used to find the Hessian matrix for our data set:

hessian <- function(h.mu, h.sigma_sqr, data){ #Defines function
  n <- length(data) #Used to get n
  
  h1 <- -n/h.sigma_sqr #Gets (1,1)
  h2 <- -(1/(h.sigma_sqr**2))*sum(log(data)-h.mu) #Gets (1,2) and (2,1)
  h3 <- (n/(2*h.sigma_sqr**2))-(1/h.sigma_sqr**3)*sum((log(data)-h.mu)**2) #Gets (2,2)
  
  H <- matrix(c(h1, h2, 
                h2, h3),
  nrow = 2, byrow = TRUE) #Makes matrix
  
  return(H) #Returns function
}

hessian(mleoutput$mu_hat, mleoutput$sigma_squared_hat, oil) #Uses function on oil data set

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