---
title: "Masked Data Likelihood Model: Components with Exponentially Distributed Lifetimes Arranged In Series Configuration"
output:
    rmarkdown::html_vignette:
        toc: true
vignette: >
  %\VignetteIndexEntry{Masked Data Likelihood Model: Components with Exponentially Distributed Lifetimes Arranged In Series Configuration}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
header-includes:
- \renewcommand{\v}[1]{\boldsymbol{#1}}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  echo = TRUE,
  comment = "#>")

# Set to TRUE to regenerate long-running simulation results
run_long <- FALSE

# maskedcauses re-exports generics from likelihood.model
# (loglik, score, hess_loglik, fit, assumptions, fim)
library(maskedcauses)
# md_encode_matrix and md_boolean_matrix_to_charsets are provided by this package
old_opts <- options(digits = 4)
```


The R package `maskedcauses` is a framework for
estimating the parameters of latent component lifetimes from *masked data*
in a series system.

Exponentially Distributed Component Lifetimes {#expo}
=====================================================
Consider a series system in which the components have exponentially distributed
lifetimes.
The $j$\textsuperscript{th} component of the $i$\textsuperscript{th} has a
lifetime distribution given by
$$
    T_{i j} \sim \operatorname{EXP}(\lambda_j)
$$
for $j=1,\ldots,m$.
Thus, $\lambda = \bigl(\lambda_1,\ldots,\lambda_m\bigr)$.
The random variable $T_{i j}$ has a reliability, pdf, and hazard function
given respectively by
\begin{align}
    \label{eq:expo_reliability}
    R_j(t|\lambda_j)   &= \exp(-\lambda_j t),\\
    \label{eq:expo_pdf}
    f_j(t|\lambda_j)   &= \lambda_j \exp(-\lambda_j t),\\
    \label{eq:expo_haz}
    h_j(\cdot|\lambda_j) &= \lambda_j
\end{align}
where $t > 0$ is the lifetime and $\lambda_j > 0$ is the failure rate of the $j$-th component.

The lifetime of the series system composed of $m$ components with exponentially
distributed lifetimes has a reliability function given by
\begin{equation}
\label{eq:sys_expo_reliability_function}
R_{T_i}(t_i|\v\lambda) = \exp \biggl(-\sum_{j=1}^{m}{\lambda_j} t_i\biggr)
\end{equation}
where $t > 0$.
\begin{proof}
By the series system reliability theorem,
$$
    R_{T_i}(t_i;\v\lambda) = \prod_{j=1}^{m} R_j(t_i;\lambda_j).
$$
Plugging in the component reliability functions given by
Equation \eqref{eq:expo_reliability} obtains the result
\begin{align*}
R_{T_i}(t_i;\v\lambda) &= \prod_{j=1}^{m} \exp(-\lambda_j t_i)\\
    &= \exp \biggl(-\sum_{j=1}^{m}{\lambda_j} t_i\biggr).
\end{align*}
\end{proof}

A series system with exponentially distributed lifetimes is also 
exponentially distributed.
\begin{theorem}
\label{thm:expo_series_family}
The random lifetime $T_i$ of a series system composed of $m$ components with 
exponentially distributed lifetimes is exponentially distributed with a failure 
rate that is the sum of the component failure rates,
$$
    T_i \sim \operatorname{EXP} \biggl(\sum_{j=1}^{m} \lambda_j\biggr).
$$
\end{theorem}
\begin{proof}
By Equation \eqref{eq:sys_expo_reliability_function}, the series system has a
reliability function in the family of the exponential distribution with a
failure rate that is the sum of the component failure rates.
\end{proof}

The series system's failure rate function is given by
\begin{equation}
\label{eq:expo_sys_failure_rate}
h(\cdot|\v\lambda) = \sum_{j=1}^{m} \lambda_j
\end{equation}
whose proof follows from the series system failure rate theorem.

We see that the system failure rate $\lambda_{\text{sys}} = \sum_{j=1}^m \lambda_j$ is *constant*,
consistent with the exponential distribution being the only continuous
distribution that has a constant failure rate.

The pdf of the series system is given by
\begin{equation}
\label{eq:expo_sys_pdf}
f_{T_i}(t_i|\v\lambda) = \biggl( \sum_{j=1}^{m} {\lambda_j} \biggr) \exp
    \biggl(-\sum_{j=1}^{m}{\lambda_j} t_i\biggr)
\end{equation}
where $t_i > 0$ is the lifetime of the system.
\begin{proof}
By definition,
$$
f_{T_i}(t_i|\v\lambda) = h_{T_i}(t_i|\v\lambda) R_{T_i}(t_i|\v\lambda).
$$
Plugging in the failure rate and reliability functions given respectively by
Equations \eqref{eq:expo_sys_failure_rate} and \eqref{eq:expo_reliability} completes
the proof.
\end{proof}

The conditional probability that component $k$ is the cause of a system failure
at time $t$ is given by
\begin{equation}
\label{eq:expo_prob_K_given_S}
f_{K_i|T_i}(k|t,\v\lambda) = f_{K_i}(k|\v\lambda) =
    \frac{\lambda_k}{\sum_{p=1}^{m} \lambda_p}
\end{equation}
where $k \in \{1,\ldots,m\}$ and $t > 0$.
\begin{proof}
By the conditional cause-of-failure theorem,
$$
f_{K_i|T_i}(k|t,\v\lambda) = \frac{h_k(t|\lambda_k)}{h_{T_i}(t;\v\lambda)}.
$$
Plugging in the failure rate of the component indexed by $k$ and the failure
rate of the system given respectively by Equations \eqref{eq:expo_sys_failure_rate}
and \eqref{eq:expo_haz} completes the proof.
\end{proof}

Due to the constant failure rates of the components, $K_i$ and $T_i$ are
mutually independent.
The joint pdf of $K_i$ and $T_i$ is given by
\begin{equation}
\label{eq:expo_joint_k_s}
f_{K_i,T_i}(k,t|\v\lambda) = \lambda_k \exp \biggl(-\sum_{j=1}^{m}{\lambda_j} t\biggr)
\end{equation}
where $k \in \{1,\ldots,m\}$ and $t > 0$.
\begin{proof}
By definition,
$$
f_{K_i,T_i}(k,t|\v\lambda) =
    f_{K_i|T_i}(k|t,\v\lambda) f_{T_i}(t|\v\lambda).
$$
Plugging in the conditional probability and the marginal probability given
respectively by Equations \eqref{eq:expo_prob_K_given_S} and
\eqref{eq:expo_sys_pdf} completes the proof.
\end{proof}

Likelihood Model {#likelihood-model}
====================================

In this study, the system is a series system with $m$
components. The true DGP for the system lifetime is in the
exponential series system family, i.e., the component lifetimes are
exponentially and independently distributed and we denote the
true parameter value by $\theta$.

The principal object of study is $\v\lambda$, which in the case of
the exponential series system family consists of $m$ rate (scale)
parameters for each component lifetime, $\v\lambda = (\lambda_1, \ldots, \lambda_m)'$.

We are interested in estimating the $\theta$ from masked data.
The masking comes in two independent forms:

- Censored system failure times, e.g., right-censoring. The system failure
time is the minimum of the component lifetimes, and it is right-censored
if the system failure time does not occur during the observation period,
$$
    T_i = \min\{\tau_i, T_{i 1}, \ldots, T_{i m}\},
$$
where $\tau_i$ is the right-censoring time for the $i$\textsuperscript{th}
observation and $T_{i 1},\ldots,T_{i m}$ are the component lifetimes
for the $i$th system.

- The cause of failure, the failed component, is masked. This masking
comes in the form of a candidate set $\mathcal{C}_i$ that, on average,
conveys information about the component cause of failure.

The candidate set $\mathcal{C}_i$ is a random variable that is a subset of
$\{1,\ldots,m\}$. The true DGP for the candidate set model has a general form
that may be denoted by
$$
    \Pr\{\mathcal{C}_i=c_i|T_1=j,\ldots,T_m,\theta,\text{other factors}\}.
$$

This is a pretty complicated looking model, and we are not even interested
in the DGP for candidate sets, except to the extent that it affects the
sampling distribution of the MLE for $\theta$.

In theory, given some candidate set model, we could construct a joint likelihood
function for the full model and jointly estimate the parameters of both the
candidate set model and $\theta$. In practice, however, this could be a very
challenging task unless we make some simplifying assumptions about the DGP for
candidate sets.

## Candidate set models
In every model we consider, we assume that the candidate set $\mathcal{C}_i$
is only a function of the component lifetimes $T_{i 1},\ldots,T_{i m}$, $\theta$,
and the right-censoring time $\tau_i$. That is, the candidate set $\mathcal{C}_i$
is independent of any other factors (or held constant for the duration of the
experiment), like ambient temperature, and these factors also have a negligible
effect on the series system lifetime and thus we can ignore them.

### Reduced likelihood model
In the Bernoulli candidate set model, we make the following assumptions about
how candidate sets are generated:

  - $C_1$: The index of the failed component is in the candidate set,
  i.e., $\Pr\{K_i \in \mathcal{C}_i\} = 1$, where
  $K_i = \arg\min_j \{ T_{i j} : j = 1,\ldots,m\}$.

  - $C_2$: The probability of $C_i$ given $K_i$ and $T_i$ is equally
  probable when the failed component varies over the components in the candidate
  set, i.e., $\Pr\{\mathcal{C}_i=c_i|K_i=j,T_i=t_i,\theta\} = \Pr\{C_i=c_i|K_i=j',T_i=t_i\}$ for
  any $j,j' \in c_i$.
    
  - $C_3$: The masking probabilities are conditionally independent of $\theta$
  given $K_i$ and $T_i$, i.e., $\Pr\{\mathcal{C}_i=c_i|K_i=j,T_i=t_i\}$ is not a
  function of $\theta$.

Using these simplifying assumptions, we can arrive at a reduced likelihood
function that only depends on $\theta$ and the observed data and as long as our
candidate set satisfies conditions $C_1$, $C_2$, and $C_3$, our reduced
likelihood function obtains the same MLEs as the full likelihood function.

We see that
$$
    \Pr\{\mathcal{C}_i=c_i,|K_i=j,T_i=t_i\} = g(c_i,t_i),
$$
since the probability cannot depend on $j$ by condition $C_2$ and cannot depend
on $\theta$ by condition $C_3$. Thus, we can write the likelihood function as
$$
    L(\theta) = \prod_{i=1}^n f_{T_i}(t_i|\theta) g(c_i,t_i).
$$

We show that $g(c_i,t_i)$ is proportional to
$$
    g(c_i,t_i) \propto \sum_{j \in c_i} f_j(t_i|\theta_j) \prod_{l=j,l \neq j}^m R_l(t_i|\theta_l),
$$
and thus the reduced likelihood is proportional to the full likelihood, yielding the same MLEs.

Note, however, that different ways in which the conditions are met will yield
MLEs with different sampling distributions, e.g., more or less efficient estimators.

### Bernoulli candidate set model #1
This is a special case of the reduced likelihood model.
In this model, we satisfy conditions $C_1$, $C_2$, and $C_3$, but we include 
each of the non-failed components with a fixed probability $p$, $0 < p < 1$.

In the simplest case, $p = 0.5$, and candidate set $c_i$ has a probability
given by
$$
\Pr\{\mathcal{C}_i=c_i|K_i=j,T_i=t_i\} =
\begin{cases}
        (1/2)^{m-1} & \text{if $j \in c_i$ and $c_i \subseteq \{1,\ldots,m\}$} \\
        0 & \text{if $j \notin c_i$}.
\end{cases}
$$
\begin{proof}
Since there are $m-1$ non-failed components (the failed component $j$ is 
in $c_i$ with probability $1$), there are $2^{m-1}$ possible
candidate sets (the size of the power set of the non-failed component indexes).
Each of these candidate sets has equal probability of occurring, and thus
the probability of any particular candidate set is $1/2^{m-1}$.
\end{proof}

### Bernoulli candidate set model #2

Now, we remove condition $C_2$. We still assume conditions $C_1$ and $C_3$,
but we allow $C_i$ to depend on the failed component $K_i$, i.e.,
$$
    \Pr\{\mathcal{C}_i=c_i|K_i=j,T_i=t_i,\theta\} \neq \Pr\{C_i=c_i|K_i=j',T_i=t_i\}
$$
for $j,j' \in c_i$.

In this case, we can write the likelihood function as
$$
    L(\theta) = \prod_{i=1}^n f_{T_i}(t_i|\theta) \prod_{j=1}^m \Pr\{K_i=j|T_i=t_i\} \prod_{c_i \in \mathcal{C}_i} g(c_i,t_i,j).
$$


Simulation {#simulation}
========================

The most straightforward series system to estimate is the series system with
exponentially distributed component lifetimes.

Suppose an exponential series system with $m$ components is parameterized by
the following R code:

```{r}
theta <- c(1,     # component 1 failure rate
           1.1,   # 3
           0.95,  # 5
           1.15,  # 6
           1.1)   # 7

m <- length(theta)
```

So, in our study, $\theta = (`r theta`)'$.
The component assigned to index $j$ has an exponentially distributed
lifetime with a failure rate $\theta_j$, e.g., $\theta_2 = `r theta[2]`$ is the
failure rate of the component indexed by $2$.

Let's simulate generating the lifetimes of the $m = `r m`$ components for this series
system:
```{r}
set.seed(7231) # set seed for reproducibility
n <- 7500
comp_times <- matrix(nrow=n,ncol=m)
for (j in 1:m)
    comp_times[,j] <- rexp(n,theta[j])
comp_times <- md_encode_matrix(comp_times,"t")
head(comp_times, 4)
```

Next, we use the function `md_series_lifetime_right_censoring` to decorate the
masked data with the right-censor-censoring time chosen by the probability
$\Pr\{T_i > \tau\} = 0.75$:
```{r right-censoring}
q <- 0.25
tau <- rep(-(1/sum(theta))*log(q),n)
data <- comp_times |> md_series_lifetime_right_censoring(tau)
latent <- attr(data, "latent")
head(data[, !colnames(data) %in% latent], 4)
```

## Masked component cause of failure
We simulate candidate sets using the Bernoulli candidate model with an
appropriate set of parameters to satisfy conditions $C_1$, $C_2$, and $C_3$:
```{r bernoulli-cand, warning=F, message=F}
p <- .3
data <- data |> md_bernoulli_cand_c1_c2_c3(p)
head(data[, paste0("q", 1:m)], 4)
```

Now, to generate candidate sets, we sample from these probabilities:
```{r cand-sampler}
data <- data |> md_cand_sampler()
data$omega <- ifelse(data$delta, "exact", "right")
display <- md_boolean_matrix_to_charsets(data, drop_set = TRUE)
latent <- attr(display, "latent")
head(display[, !colnames(display) %in% latent], 6)
```

We see that after dropping latent (unobserved) columns, we only have the right
censoring time, right censoring indicator, and the candidate sets. (Note that
this time we showed the candidate sets in a more friendly way using
`md_boolean_matrix_to_charsets`.)

Constructing the Likelihood Model {#constructing-model}
=================================
The likelihood model is a statistical model that describes the distribution of
the observed data as a function of the parameters of interest.

We construct a likelihood model for the masked data model with exponentially
distributed component lifetimes with the following code:

```{r likelihood-model}
model <- exp_series_md_c1_c2_c3()
```


Maximum likelihood estimation {#mle}
====================================
The log-likelihood for our masked data model under masking conditions
C1 (failed component in candidate set), C2 (uniform candidate set
probability), and C3 (masking independent of parameters) is given by
\begin{equation}
\label{eq:loglik_masked}
\ell(\lambda) =
    \sum_{i=1}^{n} (1-\delta_i) \log \biggl(\sum_{j \in c_i} \lambda_j \biggr) -
    \biggl( \sum_{i=1}^{n} s_i \biggr)
    \biggl( \sum_{j=1}^{m} \lambda_j \biggr).
\end{equation}
\begin{proof}
By the general masked-data log-likelihood (summing over the survival and
candidate-set contributions for each observation),
$$
\ell(\lambda) = \sum_{i=1}^n \log R(s_i;\lambda) + \sum_{i=1}^n (1-\delta_i)
    \log \biggl\{ \sum_{k\in c_i} h_k(s_i;{\lambda_k}) \biggr\}.
$$
Plugging in the component failure rate and system reliability functions given
respectively by Equations \eqref{eq:expo_haz} and
\eqref{eq:sys_expo_reliability_function} and simplifying completes the proof.
\end{proof}

The set of solutions to the MLE equations must be stationary points, i.e., a
point at which the score function of type $\mathbb{R}^m \mapsto \mathbb{R}^m$
is zero.
The $j$-th component of the output of the score function is given by
\begin{equation}
\label{eq:score_expo_j}
\frac{\partial \ell}{\partial \lambda_p} = 
  \sum_{i=1}^{n} \biggl( \sum_{j \in c_i} \lambda_j \biggr)^{-1}
  1_{\{p \in c_i \text{ and } \delta_i = 0\}} - \sum_{i=1}^{n} s_i.
\end{equation}

We may find an MLE by solving the score equation,
i.e., finding stationary points satisfying
$$
\frac{\partial \ell}{\partial \lambda_j}\Biggr|_{\hat\lambda_j} = 0
$$
for $j=1,\ldots,m$.
We approximate a solution to this problem by using the iterative
Newton-Raphson method.

The Newton-Raphson method needs the observed information matrix, which is a
function of $\lambda$ of type $\mathbb{R}^m \mapsto \mathbb{R}^{m \times m}$.
The $(j,k)$-th element of $J(\lambda)$ is given by
\begin{equation}
\label{eq:info_expo}
\frac{\partial^2 \ell}{\partial \lambda_j \partial \lambda_k} = 
  \sum_{i=1}^{n} \biggl( \sum_{j \in c_i} \lambda_j \biggr)^{-2}
  1_{\{j,k \in c_i \text{ and } \delta_i = 0\}}.
\end{equation}

## Log-likelihood of $\theta$ given masked data

The reduced log-likelihood function (the log of the kernel of the likelihood
function) is given by
$$
\ell(\theta|\text{data}) =
    -\left(\sum_{i=1}^{n} t_i\right)
    \left(\sum_{j=1}^{m} \theta_j\right) +
    \sum_{i=1}^{n} (1-\delta_i)\log\left(\sum_{j \in c_i} \theta_j\right).
$$

We compute the log-likelihood function using generic dispatch:
```{r loglike-function}
ll <- loglik(model)
```

The returned function `ll(df, par)` evaluates the log-likelihood at parameter
`par` given data `df`. For example, at the true parameter value:
```{r loglike-eval}
ll(data, theta)
```

Note that the implementation uses minimally sufficient statistics, which
improves computational efficiency.

The log-likelihood function contains the maximum amount of information
about parameter $\theta$ given the sample of masked data `data` satisfying
conditions $C_1$, $C_2$, and $C_3$.

Suppose we do not know that $\theta = (`r theta`)'$.
With the log-likelihood, we may estimate $\theta$ with $\hat\theta$ by solving
$$
\hat{\theta} = \operatorname{argmax}_{\theta \in \Omega} \ell(\theta),
$$
i.e., finding the point that *maximizes* the log-likelihood on
the observed sample `data`.
This is known as *maximum likelihood estimation* (MLE).
We typically solve for the MLE by solving
$$
\nabla \ell|_{\theta=\hat{\theta}} = 0.
$$

A popular choice is gradient ascent, which is an iterative method
based on the update rule
$$
\theta^{(n+1)} = \theta^n + \eta \ell(\theta^n),
$$
where $\eta$ is the learning rate.

We can also obtain the score (gradient) function via generic dispatch:
```{r score-function, message=FALSE, warning=FALSE}
grad <- score(model)
```

The score at the true parameter should be close to zero (at the MLE, it is exactly zero):
```{r score-eval}
grad(data, theta)
```

The `likelihood.model` framework provides analytical score and Hessian
implementations when available, falling back to numerical differentiation
otherwise.

In what follows, we use [algebraic.mle](https://github.com/queelius/algebraic.mle)
to help solve the MLE equations and display various properties of the solution.

To solve the MLE equation, we use the generic `fit()` function, which dispatches
to `fit.likelihood_model` for any object with `"likelihood_model"` in its class.
The `fit()` function returns a solver that uses `optim` internally with the
BFGS method by default.

```{r fit-model, warning=FALSE}
# Get the solver from the model using generic dispatch
solver <- fit(model)

# Solve for MLE with initial guess
theta0 <- rep(1, m)
estimate <- solver(data, par = theta0, method = "Nelder-Mead")
```

The result is an `mle` object from the `algebraic.mle` package with rich
accessor methods:

```{r mle-results}
# Print summary with confidence intervals
print(estimate)
```

We can access specific components of the MLE:
```{r mle-accessors}
# Point estimate
theta.hat <- estimate$par
cat("MLE:", round(theta.hat, 4), "\n")

# Standard errors (sqrt of diagonal of variance-covariance matrix)
cat("SE:", round(sqrt(diag(estimate$vcov)), 4), "\n")

# Log-likelihood at MLE
cat("Log-likelihood:", round(estimate$loglik, 4), "\n")
```

Recall that the true parameter is $\theta = (`r theta`)'$.

Due to sampling variability, different runs of the experiment
will result in different outcomes, i.e., $\hat{\theta}$ has a
sampling distribution.
We see that $\hat{\theta} \neq \theta$, but it is reasonably
close.
We may measure this sampling variability using the variance-covariance
matrix, bias, mean squared error (MSE), and confidence intervals.


Observation Types and Censoring {#observation-types}
====================================================

The likelihood model supports four observation types. Each arises from a
different monitoring scheme:

- **Exact** ($\omega = \text{exact}$): The system failure time $t$ is observed
  directly, e.g., continuous monitoring.
- **Right-censored** ($\omega = \text{right}$): The system is known to have
  survived past $t$, e.g., study ends before failure.
- **Left-censored** ($\omega = \text{left}$): The system is known to have
  failed before $t$, e.g., a single inspection finds it failed.
- **Interval-censored** ($\omega = \text{interval}$): The failure occurred
  in $(t, t_{\text{upper}})$, e.g., periodic inspections bracket the failure.

For the exponential series system, the individual-observation log-likelihood
contributions are:
\begin{align}
\text{Exact:} \quad & \log \lambda_c - \lambda_{\text{sys}} t \\
\text{Right:} \quad & -\lambda_{\text{sys}} t \\
\text{Left:} \quad & \log \lambda_c + \log(1 - e^{-\lambda_{\text{sys}} \tau})
  - \log \lambda_{\text{sys}} \\
\text{Interval:} \quad & \log \lambda_c - \lambda_{\text{sys}} a
  + \log(1 - e^{-\lambda_{\text{sys}}(b - a)}) - \log \lambda_{\text{sys}}
\end{align}
where $\lambda_c = \sum_{j \in c_i} \lambda_j$ and
$\lambda_{\text{sys}} = \sum_j \lambda_j$.

The left-censored contribution can be interpreted as
$\log w_c + \log F(\tau)$, where $w_c = \lambda_c / \lambda_{\text{sys}}$
is the candidate cause weight and $F(\tau)$ is the system CDF.

**The exponential model is the only one where all four types have fully
analytical loglik, score, AND Hessian.** The Weibull models require
numerical integration or `numDeriv` for left/interval types.

## Observe functors

The package provides composable observation functors for generating data under
different monitoring schemes:

```{r observe-functors}
# Periodic inspection every 0.5 time units, study ends at tau = 5
obs_periodic <- observe_periodic(delta = 0.5, tau = 5)

# Single inspection at tau = 3
obs_left <- observe_left_censor(tau = 3)

# Mix of continuous and periodic monitoring
obs_mixed <- observe_mixture(
  observe_right_censor(tau = 5),
  observe_left_censor(tau = 3),
  weights = c(0.7, 0.3)
)
```

## Generating mixed-censoring data

```{r mixed-censoring-data}
gen <- rdata(model)

# Periodic inspections
set.seed(7231)
df_periodic <- gen(theta, n = 500, p = 0.3,
                   observe = observe_periodic(delta = 0.5, tau = 5))
cat("Periodic inspection observation types:\n")
print(table(df_periodic$omega))

# Mixed monitoring
set.seed(7231)
df_mixed <- gen(theta, n = 500, p = 0.3,
                observe = observe_mixture(
                  observe_right_censor(tau = 5),
                  observe_left_censor(tau = 3),
                  weights = c(0.7, 0.3)
                ))
cat("\nMixed monitoring observation types:\n")
print(table(df_mixed$omega))
```

## Likelihood evaluation on mixed-censoring data

All four observation types contribute analytically to the loglik, score, and
Hessian for the exponential model:

```{r mixed-censoring-eval}
ll_fn <- loglik(model)
scr_fn <- score(model)
hess_fn <- hess_loglik(model)

# Evaluate at true parameters
ll_val <- ll_fn(df_periodic, theta)
scr_val <- scr_fn(df_periodic, theta)
hess_val <- hess_fn(df_periodic, theta)

cat("Log-likelihood (periodic):", round(ll_val, 4), "\n")
cat("Score (periodic):", round(scr_val, 4), "\n")
cat("Hessian eigenvalues:", round(eigen(hess_val)$values, 4), "\n")
```

We verify that the analytical score is consistent with numerical
differentiation:

```{r score-verify}
scr_numerical <- numDeriv::grad(
  func = function(th) ll_fn(df_periodic, th),
  x = theta
)
cat("Max |analytical - numerical| score:",
    formatC(max(abs(scr_val - scr_numerical)), format = "e", digits = 2), "\n")
```


## Monte Carlo simulation study

To understand the sampling properties of the MLE, we conduct a Monte Carlo
simulation study. We repeatedly:

1. Generate masked data from the true model
2. Fit the likelihood model to obtain $\hat{\theta}$
3. Compute asymptotic confidence intervals

From the collection of estimates, we can compute:

- **Bias**: $\text{E}[\hat{\theta}] - \theta$
- **Variance**: $\text{Var}[\hat{\theta}]$
- **MSE**: $\text{E}[(\hat{\theta} - \theta)^2] = \text{Bias}^2 + \text{Var}$
- **Coverage probability**: Proportion of CIs containing the true value

```{r load-precomputed, include=FALSE, eval=!run_long}
# Load pre-computed simulation results when not re-running
list2env(readRDS("precomputed_results.rds"), envir = environment())
```

```{r monte-carlo-setup, cache=TRUE, eval=run_long}
set.seed(7231)

# Simulation parameters
B <- 200       # Number of Monte Carlo replications
alpha <- 0.05  # Significance level for CIs

# Storage for results
estimates <- matrix(NA, nrow = B, ncol = m)
se_estimates <- matrix(NA, nrow = B, ncol = m)
ci_lower <- matrix(NA, nrow = B, ncol = m)
ci_upper <- matrix(NA, nrow = B, ncol = m)
converged <- logical(B)
```

```{r monte-carlo-run, cache=TRUE, warning=FALSE, eval=run_long}
for (b in 1:B) {
    # Generate component lifetimes
    comp_times_b <- matrix(nrow = n, ncol = m)
    for (j in 1:m) {
        comp_times_b[, j] <- rexp(n, theta[j])
    }
    comp_times_b <- md_encode_matrix(comp_times_b, "t")

    # Apply masking pipeline
    data_b <- comp_times_b |>
        md_series_lifetime_right_censoring(tau) |>
        md_bernoulli_cand_c1_c2_c3(p) |>
        md_cand_sampler()
    data_b$omega <- ifelse(data_b$delta, "exact", "right")

    # Fit model
    tryCatch({
        result_b <- solver(data_b, par = theta0, method = "Nelder-Mead")
        estimates[b, ] <- result_b$par
        se_estimates[b, ] <- sqrt(diag(result_b$vcov))

        # Asymptotic Wald CIs
        z <- qnorm(1 - alpha/2)
        ci_lower[b, ] <- result_b$par - z * sqrt(diag(result_b$vcov))
        ci_upper[b, ] <- result_b$par + z * sqrt(diag(result_b$vcov))
        converged[b] <- result_b$converged
    }, error = function(e) {
        converged[b] <<- FALSE
    })
}

cat("Convergence rate:", mean(converged, na.rm = TRUE), "\n")
```

### Bias, Variance, and MSE

```{r monte-carlo-results}
# Compute summary statistics (only for converged runs)
valid <- converged & !is.na(estimates[, 1])
est_valid <- estimates[valid, , drop = FALSE]

bias <- colMeans(est_valid) - theta
variance <- apply(est_valid, 2, var)
mse <- bias^2 + variance
rmse <- sqrt(mse)

# Create results table
results_df <- data.frame(
    Component = 1:m,
    True = theta,
    Mean_Est = colMeans(est_valid),
    Bias = bias,
    Variance = variance,
    MSE = mse,
    RMSE = rmse,
    Rel_Bias_Pct = 100 * bias / theta
)

knitr::kable(results_df, digits = 4,
             caption = "Monte Carlo Results: Bias, Variance, and MSE",
             col.names = c("Component", "True θ", "Mean θ̂", "Bias",
                          "Variance", "MSE", "RMSE", "Rel. Bias %"))
```

### Confidence Interval Coverage

The asymptotic $(1-\alpha)$% Wald confidence interval is:
$$
\hat{\theta}_j \pm z_{1-\alpha/2} \cdot \text{SE}(\hat{\theta}_j)
$$

We assess coverage probability: the proportion of intervals that contain
the true parameter value.

```{r ci-coverage}
# Compute coverage for each component
coverage <- numeric(m)
for (j in 1:m) {
    valid_j <- valid & !is.na(ci_lower[, j]) & !is.na(ci_upper[, j])
    covered <- (ci_lower[valid_j, j] <= theta[j]) & (theta[j] <= ci_upper[valid_j, j])
    coverage[j] <- mean(covered)
}

# Mean CI width
mean_width <- colMeans(ci_upper[valid, ] - ci_lower[valid, ], na.rm = TRUE)

coverage_df <- data.frame(
    Component = 1:m,
    True = theta,
    Coverage = coverage,
    Nominal = 1 - alpha,
    Mean_Width = mean_width
)

knitr::kable(coverage_df, digits = 4,
             caption = paste0("Coverage Probability of ", 100*(1-alpha), "% Confidence Intervals"),
             col.names = c("Component", "True θ", "Coverage", "Nominal", "Mean Width"))
```

### Sampling Distribution Visualization

```{r sampling-dist-plot, fig.width=8, fig.height=4}
# Plot sampling distributions
oldpar <- par(mfrow = c(1, min(m, 5)), mar = c(4, 4, 2, 1))
on.exit(par(oldpar))
for (j in 1:min(m, 5)) {
    hist(est_valid[, j], breaks = 20, probability = TRUE,
         main = paste0("Component ", j),
         xlab = expression(hat(theta)[j]),
         col = "lightblue", border = "white")
    abline(v = theta[j], col = "red", lwd = 2, lty = 2)
    abline(v = mean(est_valid[, j]), col = "blue", lwd = 2)
    legend("topright", legend = c("True", "Mean Est."),
           col = c("red", "blue"), lty = c(2, 1), lwd = 2, cex = 0.7)
}
```

### Summary

The Monte Carlo simulation ($n = 7500$, $B = 200$, $p = 0.3$, $\sim 25\%$
censoring) demonstrates that the MLE for the exponential series system with
masked data satisfying conditions C1--C3 has desirable large-sample properties:

- **Essentially unbiased.** All relative biases are below 0.7%, with the largest
  being component 5 at 0.67%. The squared-bias contribution to MSE is negligible
  (e.g., component 1: $\text{bias}^2 \approx 2 \times 10^{-7}$ vs
  $\text{Var} = 0.0019$). The MLE is consistent and approximately unbiased at
  this sample size.

- **Uniform relative precision.** RMSE ranges from 0.042 to 0.046 across
  components --- roughly 4--5% of the true rates. This uniformity is consistent
  with the asymptotic efficiency of the MLE: components with higher failure rates
  (and thus higher variance) also have proportionally wider CIs.

- **Coverage near nominal.** Coverage ranges from 93.4% to 96.5% against a
  nominal 95%. Component 4 (93.4%) slightly undercovers, suggesting the Wald
  interval may be mildly liberal for some components at this sample size. This
  is a known finite-sample property of Wald intervals for rate parameters.

- **CI width as a design metric.** Mean CI widths of 0.170--0.180 (about
  $\pm 8.5\%$ of the true rate) provide a concrete basis for sample size
  planning. Since CI width scales as $1/\sqrt{n}$, halving the width requires
  roughly $4\times$ the sample size.

Sensitivity Analysis {#sensitivity}
===================================

The MLE properties depend on several factors: sample size, masking probability,
and right-censoring rate. In this section, we study how these factors affect
estimation accuracy.

## Effect of Masking Probability

The masking probability $p$ controls how much information about the failed
component is obscured. When $p = 0$, only the failed component is in the
candidate set (perfect information). As $p$ increases, more non-failed
components are included, making estimation more difficult.

```{r masking-sensitivity-setup, cache=TRUE, eval=run_long}
set.seed(7231)

# Use smaller sample for sensitivity study
n_sens <- 500
B_sens <- 100

# Masking probabilities to test
p_values <- seq(0, 0.5, by = 0.1)

# Fixed right-censoring (25% censored)
q_sens <- 0.25
tau_sens <- rep(-(1/sum(theta))*log(q_sens), n_sens)

# Storage
mask_results <- list()
```

```{r masking-sensitivity-run, cache=TRUE, warning=FALSE, message=FALSE, eval=run_long}
for (p_idx in seq_along(p_values)) {
    p_curr <- p_values[p_idx]
    est_p <- matrix(NA, nrow = B_sens, ncol = m)

    for (b in 1:B_sens) {
        # Generate data
        comp_b <- matrix(nrow = n_sens, ncol = m)
        for (j in 1:m) comp_b[, j] <- rexp(n_sens, theta[j])
        comp_b <- md_encode_matrix(comp_b, "t")

        data_b <- comp_b |>
            md_series_lifetime_right_censoring(tau_sens) |>
            md_bernoulli_cand_c1_c2_c3(p_curr) |>
            md_cand_sampler()
        data_b$omega <- ifelse(data_b$delta, "exact", "right")

        tryCatch({
            fit_b <- solver(data_b, par = theta0, method = "Nelder-Mead")
            if (fit_b$converged) est_p[b, ] <- fit_b$par
        }, error = function(e) NULL)
    }

    # Compute statistics
    valid_p <- !is.na(est_p[, 1])
    mask_results[[p_idx]] <- list(
        p = p_curr,
        bias = colMeans(est_p[valid_p, , drop = FALSE]) - theta,
        variance = apply(est_p[valid_p, , drop = FALSE], 2, var),
        mse = (colMeans(est_p[valid_p, , drop = FALSE]) - theta)^2 +
              apply(est_p[valid_p, , drop = FALSE], 2, var)
    )
}
```

```{r masking-sensitivity-plot, fig.width=8, fig.height=5}
# Extract bias and MSE for plotting
bias_mat <- sapply(mask_results, function(x) mean(abs(x$bias)))
mse_mat <- sapply(mask_results, function(x) mean(x$mse))

oldpar <- par(mfrow = c(1, 2), mar = c(4, 4, 2, 1))
on.exit(par(oldpar))

# Mean absolute bias vs masking probability
plot(p_values, bias_mat, type = "b", pch = 19, col = "blue",
     xlab = "Masking Probability (p)",
     ylab = "Mean Absolute Bias",
     main = "Bias vs Masking Probability")
grid()

# Mean MSE vs masking probability
plot(p_values, mse_mat, type = "b", pch = 19, col = "red",
     xlab = "Masking Probability (p)",
     ylab = "Mean MSE",
     main = "MSE vs Masking Probability")
grid()
```

```{r masking-sensitivity-table}
mask_df <- data.frame(
    p = p_values,
    Mean_Abs_Bias = bias_mat,
    Mean_MSE = mse_mat,
    Mean_RMSE = sqrt(mse_mat)
)

knitr::kable(mask_df, digits = 4,
             caption = "Effect of Masking Probability on Estimation Accuracy",
             col.names = c("Masking Prob.", "Mean |Bias|", "Mean MSE", "Mean RMSE"))
```

The simulation reveals several patterns in how masking degrades estimation
quality:

- **MSE increases monotonically** from 0.015 ($p = 0$) to 0.062 ($p = 0.5$)
  --- a $4.2\times$ degradation. RMSE roughly doubles from 0.12 to 0.25,
  meaning the typical estimation error doubles when moving from unmasked to
  maximally masked candidate sets.

- **Variance drives the degradation, not bias.** Mean $|\text{bias}|$
  fluctuates between 0.008 and 0.026 with no systematic trend, while variance
  grows steadily with $p$. This is consistent with the theoretical result that
  conditions C1--C3 preserve the unbiasedness of the MLE --- masking erodes
  *precision* without introducing systematic error.

- **The marginal cost of masking accelerates.** MSE roughly doubles from $p = 0$
  to $p = 0.3$ ($0.015 \to 0.033$), then roughly doubles again from $p = 0.3$
  to $p = 0.5$ ($0.033 \to 0.062$). For practitioners, this means moderate
  masking ($p \leq 0.3$) is far less costly than heavy masking.

- **Residual MSE at $p = 0$.** Even with no masking of non-failed components,
  MSE of 0.015 remains due to right-censoring alone --- a useful baseline for
  isolating the censoring contribution.

## Effect of Right-Censoring Rate

Right-censoring reduces the number of exact failure times observed. When a
system is right-censored, we know the system survived beyond the censoring
time, but we don't observe the actual failure time or failed component.

```{r censoring-sensitivity-setup, cache=TRUE, eval=run_long}
set.seed(7231)

# Censoring quantiles (proportion surviving past tau)
# q = 0.1 means 10% survive past tau (light censoring, ~10% censored)
# q = 0.9 means 90% survive past tau (heavy censoring, ~90% censored)
q_values <- seq(0.1, 0.9, by = 0.1)  # Survival probabilities (matches simulation framework)

# Fixed masking probability
p_cens <- 0.2

# Storage
cens_results <- list()
```

```{r censoring-sensitivity-run, cache=TRUE, warning=FALSE, message=FALSE, eval=run_long}
for (q_idx in seq_along(q_values)) {
    q_curr <- q_values[q_idx]
    tau_curr <- rep(-(1/sum(theta))*log(q_curr), n_sens)
    est_q <- matrix(NA, nrow = B_sens, ncol = m)

    for (b in 1:B_sens) {
        # Generate data
        comp_b <- matrix(nrow = n_sens, ncol = m)
        for (j in 1:m) comp_b[, j] <- rexp(n_sens, theta[j])
        comp_b <- md_encode_matrix(comp_b, "t")

        data_b <- comp_b |>
            md_series_lifetime_right_censoring(tau_curr) |>
            md_bernoulli_cand_c1_c2_c3(p_cens) |>
            md_cand_sampler()
        data_b$omega <- ifelse(data_b$delta, "exact", "right")

        tryCatch({
            fit_b <- solver(data_b, par = theta0, method = "Nelder-Mead")
            if (fit_b$converged) est_q[b, ] <- fit_b$par
        }, error = function(e) NULL)
    }

    # Compute statistics
    valid_q <- !is.na(est_q[, 1])
    cens_rate <- mean(data_b$omega == "right")  # Actual censoring rate
    cens_results[[q_idx]] <- list(
        q = q_curr,
        cens_rate = cens_rate,
        bias = colMeans(est_q[valid_q, , drop = FALSE]) - theta,
        variance = apply(est_q[valid_q, , drop = FALSE], 2, var),
        mse = (colMeans(est_q[valid_q, , drop = FALSE]) - theta)^2 +
              apply(est_q[valid_q, , drop = FALSE], 2, var)
    )
}
```

```{r censoring-sensitivity-plot, fig.width=8, fig.height=5}
# Extract statistics
cens_rates <- sapply(cens_results, function(x) x$cens_rate)
bias_cens <- sapply(cens_results, function(x) mean(abs(x$bias)))
mse_cens <- sapply(cens_results, function(x) mean(x$mse))

oldpar <- par(mfrow = c(1, 2), mar = c(4, 4, 2, 1))
on.exit(par(oldpar))

# Mean absolute bias vs censoring rate
plot(cens_rates * 100, bias_cens, type = "b", pch = 19, col = "blue",
     xlab = "Censoring Rate (%)",
     ylab = "Mean Absolute Bias",
     main = "Bias vs Censoring Rate")
grid()

# Mean MSE vs censoring rate
plot(cens_rates * 100, mse_cens, type = "b", pch = 19, col = "red",
     xlab = "Censoring Rate (%)",
     ylab = "Mean MSE",
     main = "MSE vs Censoring Rate")
grid()
```

```{r censoring-sensitivity-table}
cens_df <- data.frame(
    Cens_Rate_Pct = round(cens_rates * 100, 1),
    Mean_Abs_Bias = bias_cens,
    Mean_MSE = mse_cens,
    Mean_RMSE = sqrt(mse_cens)
)

knitr::kable(cens_df, digits = 4,
             caption = "Effect of Right-Censoring Rate on Estimation Accuracy",
             col.names = c("Censoring %", "Mean |Bias|", "Mean MSE", "Mean RMSE"))
```

The simulation reveals that censoring is the dominant source of information
loss in this model:

- **MSE grows from 0.019 (10% censored) to 0.183 (90% censored)** --- a
  $9.5\times$ degradation, far exceeding the $4.2\times$ range from masking.
  Censoring is more damaging because it eliminates *both* the failure time and
  the candidate set, whereas masking only dilutes the candidate set.

- **Convex MSE growth, accelerating sharply above 70%.** MSE grows gradually
  from 10% to 70% censoring ($0.019 \to 0.059$, roughly linear), then
  accelerates: $0.088$ at 80% and $0.183$ at 90%. The inflection point is
  near 70--80% censoring, beyond which each additional percentage point of
  censoring is dramatically more damaging.

- **Bias remains modest even under extreme censoring.** Mean $|\text{bias}|$
  grows from 0.005 (10% censoring) to 0.034 (90% censoring) but stays small
  relative to the true rates ($\sim 3\%$ at worst). The MLE remains
  approximately unbiased --- it is consistent but increasingly inefficient.

- **Robustness under extreme censoring.** At 90% censoring, only $\sim 50$ of
  500 observations are exact failures, yet the MLE still converges. The
  resulting estimates are imprecise (RMSE $\approx 0.43$, or $\sim 40\%$ of the
  true rates), but the likelihood-based approach remains functional.

```{r save-precomputed, include=FALSE, eval=run_long}
# Save simulation results for future fast vignette builds
saveRDS(list(
    estimates = estimates, se_estimates = se_estimates,
    ci_lower = ci_lower, ci_upper = ci_upper, converged = converged,
    B = B, alpha = alpha,
    mask_results = mask_results, p_values = p_values,
    cens_results = cens_results, q_values = q_values,
    theta = theta, m = m, n_sens = n_sens, B_sens = B_sens, p_cens = p_cens
), "precomputed_results.rds")
```

## Practical Recommendations

Based on these sensitivity analyses:

1. **Prioritize reducing censoring over masking.** Censoring degrades MSE by
   $9.5\times$ across its range versus $4.2\times$ for masking. In experimental
   design, extending the observation window to reduce censoring below 50% yields
   larger gains than improving diagnostic resolution to reduce masking.

2. **Moderate masking is tolerable.** MSE only doubles from $p = 0$ to
   $p = 0.3$. If reducing the masking probability below 0.3 requires expensive
   diagnostic equipment, the incremental benefit may not justify the cost.

3. **The 70% censoring threshold.** Below $\sim 70\%$ censoring, MSE growth is
   roughly linear and moderate. Above 70%, MSE accelerates sharply --- avoid
   experimental designs where more than two-thirds of systems survive past the
   observation window.

4. **CI widths for sample size planning.** At $n = 7500$ with $p = 0.3$ and
   25% censoring, 95% CI widths are $\sim 0.17$ for rates near 1.0. Width
   scales as $1/\sqrt{n}$, so $n = 1875$ gives widths of $\sim 0.34$ (double),
   and $n = 30000$ gives $\sim 0.085$ (half).

```{r cleanup, include=FALSE}
options(old_opts)
```

