Consider the model \(f_6(t; a_1, a_2) = 100 + a_1t + a_2(1-e^{-0.0003t})\) where \(t \geq 0\). The function \(f_6\) models the brightness of a lightbulb, measured as a percent of the original intensity of the lightbulb, given the number of hours the lightbulb as been on, \(t\).
In order to fit \(f_6\) to a list of 44 measurements, \((t_i, y_i)\) we need to find the first and second partial derivatives of the loglikelihood function for the residuals. Assuming the residuals (or errors) are independent and normally distributed (with mean 0 and standard deviation 1), the loglikelihood function for these errors is \(\ell_6(a_1,a_2; \mathbf{t},\mathbf{y}) = 44\ln\left(\frac{1}{\sqrt{2\pi}}\right) + \sum_{i=1}^{44} \left(- \frac{1}{2}\right)(y_i - 100 - a_1t_i - a_2(1-e^{-0.0003t_i}))^2\).
We compute \(\frac{\partial \ell_6}{\partial a_1}\) by treating \(a_2\) as a constant and differentiate with respect to \(a_1\) using the derivative rules.
\[\begin{align*} \frac{\partial \ell_6}{\partial a_1} &= \frac{\partial}{\partial a_1}\left(44\ln\left(\frac{1}{\sqrt{2\pi}}\right)\right) + \sum_{i=1}^{44} \frac{\partial}{\partial a_1}\left(\left(- \frac{1}{2}\right)(y_i - 100 - a_1t_i - a_2(1-e^{-0.0003t_i}))^2\right) &\text{(sum/difference rule)}\\ &= 0 + \sum_{i=1}^{44} \left(- \frac{1}{2}\right)\frac{\partial}{\partial a_1}\left((y_i - 100 - a_1t_i - a_2(1-e^{-0.0003t_i}))^2\right) &\text{(constant and constant multiple rules)}\\ &= \sum_{i=1}^{44} \left(- \frac{1}{2}\right)\left(2(y_i - 100 - a_1t_i - a_2(1-e^{-0.0003t_i}))(-t_i)\right) &\text{(chain rule)}\\ &= \sum_{i=1}^{44} t_i(y_i - 100 - a_1t_i - a_2(1-e^{-0.0003t_i})) &\text{(simplify; multiply $-\frac{1}{2}$, 2, and $-t_i$)}\\ \end{align*}\]
We can simplify \(\frac{\partial \ell_6}{\partial a_1}\) using the properties of sums.
\[\begin{align*} \frac{\partial \ell_6}{\partial a_1} &= \sum_{i=1}^{44} t_i(y_i - 100 - a_1t_i - a_2(1-e^{-0.0003t_i})) & \\ &= \sum_{i=1}^{44} t_i(y_i - 100) - a_1t_i^2 - a_2t_i(1-e^{-0.0003t_i}) &\text{(distribution inside the sum)}\\ &= \sum_{i=1}^{44} t_i(y_i - 100) + \sum_{i=1}^{44} -a_1t_i^2 + \sum_{i=1}^{44} -a_2t_i(1-e^{-0.0003t_i}) &\text{(break the sum into three sums)}\\ &= \sum_{i=1}^{44} t_i(y_i - 100) - a_1\sum_{i=1}^{44}t_i^2 - a_2\sum_{i=1}^{44}t_i(1-e^{-0.0003t_i}) &\text{(pull constant multiples out of the sums)}\\ \end{align*}\]
We see \(\frac{\partial \ell_6}{\partial a_1}\) is of the form a constant, subtract a constant multiple of \(a_1\), subtract a constant multiple of \(a_2\). \[\frac{\partial \ell_6}{\partial a_1} = \left(\sum_{i=1}^{44} t_i(y_i - 100)\right) - a_1\left(\sum_{i=1}^{44}t_i^2\right) - a_2\left(\sum_{i=1}^{44}t_i(1-e^{-0.0003t_i})\right)\]
These constants depend on our data and can be computed in R. Let \(A = \sum_{i=1}^{44} t_i(y_i - 100) \approx 83738.9\), \(B = \sum_{i=1}^{44}t_i^2 \approx 328767530\), and \(C = \sum_{i=1}^{44}t_i(1-e^{-0.0003t_i}) \approx 59520.5\).
library(data4led)
bulb <- led_bulb(1,seed=2021)
t <- bulb$hours
y <- bulb$percent_intensity
A <- sum((y-100)*t)
A## [1] 83738.93B <- sum(t^2)
B## [1] 328767530C <- sum((1-exp(-0.0003*t))*t)
C## [1] 59520.54C4.6 <- -sum((y-100)*(1-exp(-0.0003*t)))
C5.6 <- sum((1-exp(-0.0003*t))^2)We compute \(\frac{\partial \ell_6}{\partial a_2}\) by treating \(a_1\) as a constant and differentiate with respect to \(a_2\) using the derivative rules.
\[\begin{align*} \frac{\partial \ell_6}{\partial a_2} &= \frac{\partial}{\partial a_2}\left(44\ln\left(\frac{1}{\sqrt{2\pi}}\right)\right) + \sum_{i=1}^{44} \frac{\partial}{\partial a_2}\left(\left(- \frac{1}{2}\right)(y_i - 100 - a_1t_i - a_2(1-e^{-0.0003t_i}))^2\right) &\text{(sum/difference rule)}\\ &= 0 + \sum_{i=1}^{44} \left(- \frac{1}{2}\right)\frac{\partial}{\partial a_2}\left((y_i - 100 - a_1t_i - a_2(1-e^{-0.0003t_i}))^2\right) &\text{(constant and constant multiple rules)}\\ &= \sum_{i=1}^{44} \left(- \frac{1}{2}\right)\left(2(y_i - 100 - a_1t_i - a_2(1-e^{-0.0003t_i}))(-(1-e^{-0.0003t_i}))\right) &\text{(chain rule)}\\ &= \sum_{i=1}^{44} (1-e^{-0.0003t_i})(y_i - 100 - a_1t_i - a_2(1-e^{-0.0003t_i})) &\text{(simplify; multiply $-\frac{1}{2}$, 2, and $-(1-e^{-0.0003t_i})$)}\\ \end{align*}\]
We can simplify \(\frac{\partial \ell_6}{\partial a_2}\) using the properties of sums.
\[\begin{align*} \frac{\partial \ell_6}{\partial a_2} &= \sum_{i=1}^{44} (1-e^{-0.0003t_i})(y_i - 100 - a_1t_i - a_2(1-e^{-0.0003t_i})) & \\ &= \sum_{i=1}^{44} (1-e^{-0.0003t_i})(y_i - 100) - a_1t_i(1-e^{-0.0003t_i}) - a_2(1-e^{-0.0003t_i})^2 &\text{(distribution inside the sum)}\\ &= \sum_{i=1}^{44} (1-e^{-0.0003t_i})(y_i - 100) + \sum_{i=1}^{44} -a_1t_i(1-e^{-0.0003t_i}) + \sum_{i=1}^{44} -a_2(1-e^{-0.0003t_i})^2 &\text{(break the sum into three sums)}\\ &= \sum_{i=1}^{44} (1-e^{-0.0003t_i})(y_i - 100) - a_1\sum_{i=1}^{44}t_i(1-e^{-0.0003t_i}) - a_2\sum_{i=1}^{44}(1-e^{-0.0003t_i})^2 &\text{(pull constant multiples out of the sums)}\\ \end{align*}\]
We see \(\frac{\partial \ell_6}{\partial a_2}\) is of the form a constant, subtract a constant multiple of \(a_1\), subtract a constant multiple of \(a_2\). \[\frac{\partial \ell_6}{\partial a_2} = \left(\sum_{i=1}^{44} (1-e^{-0.0003t_i})(y_i - 100)\right) - a_1\left(\sum_{i=1}^{44}t_i(1-e^{-0.0003t_i})\right) - a_2\left(\sum_{i=1}^{44}(1-e^{-0.0003t_i})^2\right)\]
These constants depend on our data and can be computed in R. Let \(D = \sum_{i=1}^{44} (1-e^{-0.0003t_i})(y_i - 100) \approx -15.691\), \(C = \sum_{i=1}^{44}t_i(1-e^{-0.0003t_i}) \approx 59520.5\), and \(E = \sum_{i=1}^{44}(1-e^{-0.0003t_i})^2 \approx 10.998\). Notice that \(C\) is a constant in \(\frac{\partial \ell_6}{\partial a_2}\) and \(\frac{\partial \ell_6}{\partial a_1}\).
D <- -sum((y-100)*(1-exp(-0.0003*t)))
D## [1] -15.6911C <- sum((1-exp(-0.0003*t))*t)
C## [1] 59520.54E <- sum((1-exp(-0.0003*t))^2)
E## [1] 10.99773We will use the first partials of \(\ell_6\) to calculate the second partials of \(\ell_6\). We only need to compute \(\frac{\partial^2 \ell_6}{\partial a_1^2}\), \(\frac{\partial^2 \ell_6}{\partial a_2^2}\) and \(\frac{\partial^2 \ell_6}{\partial a_1a_2}\) since the mix partials are equal for \(\ell_6\).
We compute \(\frac{\partial^2 \ell_6}{\partial a_1^2}\) by taking the derivative of \(\frac{\partial \ell_6}{\partial a_1}\) with respect to \(a_1\). \[\frac{\partial^2 \ell_6}{\partial a_1^2} = \frac{\partial}{\partial a_1} \left(\left(\sum_{i=1}^{44} t_i(y_i - 100)\right) - a_1\left(\sum_{i=1}^{44}t_i^2\right) - a_2\left(\sum_{i=1}^{44}t_i(1-e^{-0.0003t_i})\right)\right) = 0 - \left(\sum_{i=1}^{44}t_i^2\right) - 0\] We see \(\frac{\partial^2 \ell_6}{\partial a_1^2} = -\sum_{i=1}^{44}t_i^2 \approx -99542\).
-B## [1] -328767530We compute \(\frac{\partial^2 \ell_6}{\partial a_2^2}\) by taking the derivative of \(\frac{\partial \ell_6}{\partial a_2}\) with respect to \(a_2\). \[\frac{\partial^2 \ell_6}{\partial a_2^2} = \frac{\partial}{\partial a_2} \left(\left(\sum_{i=1}^{44} (1-e^{-0.0003t_i})(y_i - 100)\right) - a_1\left(\sum_{i=1}^{44}t_i(1-e^{-0.0003t_i})\right) - a_2\left(\sum_{i=1}^{44}(1-e^{-0.0003t_i})^2\right)\right) = 0 - 0 - \sum_{i=1}^{44}(1-e^{-0.0003t_i})^2\] We see \(\frac{\partial^2 \ell_6}{\partial a_2^2} = - \sum_{i=1}^{44}(1-e^{-0.0003t_i})^2 \approx -10.998\).
-E## [1] -10.99773We compute \(\frac{\partial^2 \ell_6}{\partial a_1a_2}\) by taking the derivative of \(\frac{\partial \ell_6}{\partial a_2}\) with respect to \(a_1\). \[\frac{\partial^2 \ell_6}{\partial a_1a_2} = \frac{\partial}{\partial a_1} \left(\left(\sum_{i=1}^{44} (1-e^{-0.0003t_i})(y_i - 100)\right) - a_1\left(\sum_{i=1}^{44}t_i(1-e^{-0.0003t_i})\right) - a_2\left(\sum_{i=1}^{44}(1-e^{-0.0003t_i})^2\right)\right) = 0 - \sum_{i=1}^{44}t_i(1-e^{-0.0003t_i}) - 0\] We see \(\frac{\partial^2 \ell_6}{\partial a_2^2} = - \sum_{i=1}^{44}t_i(1-e^{-0.0003t_i}) \approx -59520.5\).
-C## [1] -59520.54