STA 506 Final Examination

Working Model for the Final Exam

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

The Kumaraswamy distribution is a two-parameter continuous probability distribution defined on the interval (0, 1). It is often used as an alternative to the Beta distribution due to its simple closed-form expressions for the cumulative distribution function (CDF) and quantile function. It is commonly used in

• Hydrology: Modeling rainfall, streamflow, or other bounded natural phenomena

• Economics: Income shares, proportions, or bounded indices

• Monte Carlo simulation: Efficient random variate generation (via inverse transform)

• Machine learning: Output layer for bounded targets, prior distributions in Bayesian models

• Reliability engineering: Modeling failure rates of systems with bounded lifetimes

Let \(X\) be the Kumaraswamy random variable with Cumulative Distribution Function (CDF)

\[ F(x; a, b) = 1 - (1 - x^a)^b \]

where \(a > 0\) and \(b > 0\) unknown parameters and \(0 < x < 1\).

The following are two special case of the Kumaraswamy distribution:

1. Uniform Distribution: When \(a = 1\) and \(b = 1\), the Kumaraswamy distribution becomes a uniform distribution over \([0, 1]\) with CDF \(F(x) = x\).

2. Power Distribution: when \(b = 1\) and \(a > 0\), the Kumaraswamy distribution becomes a power distribution over \([0, 1]\) with CDF \(F(x) = x^a\).

This final exam focuses on inferences of Kumaraswamy distribution and related data analysis.

Part A: Methodological Derivations

Problem A1:

Show that the density function of the Kumaraswamy distribution is:

\[ f(x; a, b) = ab \, x^{a-1} (1 - x^a)^{b-1}. \] Answer:

The CDF of the Kumaraswamy distribution is \[ F(x; a, b) = 1 - (1 - x^a)^b \quad, 0<x<1,a>0,b>0. \] To find the PDF \(f(x; a, b)\), we compute: \[ f(x; a, b) = \frac{d}{dx} F(x; a, b) = \frac{d}{dx} \left[ 1 - (1 - x^a)^b \right]. \] Apply the sum rule to simplifies the exmpression to: \[ f(x; a, b) = - \frac{d}{dx} (1 - x^a)^b \] Then using the chain rule, treat \((1 - x^a)\) as the inner function: \[ f(x; a, b) = - \left[ b(1 - x^a)^{b-1} \cdot \frac{d}{dx}(1 - x^a) \right] \] Therefore: \[ f(x; a, b) = ab \, x^{a-1} (1 - x^a)^{b-1}, 0 < x < 1. \] So the density is: \[ f(x; a, b) = ab \, x^{a-1} (1 - x^a)^{b-1}. \]

Problem A2:

Let \(\{x_1, x_2, \cdots, x_n \}\) be an i.i.d. random sample taken from a population that follows the aove 2-parameter Kumaraswamy distribution. Write out the loglikelihood function of \(a\) and \(b\), denoted by \(\ell(a,b)\), based on the above random sample and derive the gradient vector \([\ell_a^\prime(a,b), \ell_b^\prime(a,b)]\), the first order partial derivative of the log-likelihood with respect to parameters \(a\) and \(b\).

Answer:

Given the probability density function (PDF): \[ f(x; a, b) = ab \, x^{a-1} (1 - x^a)^{b-1} \] Taking logs gives he Log-Likelihood Function: \[ \ell(a, b) = n \ln(a) + n \ln(b) + (a-1) \sum_{i=1}^{n} \ln(x_i) + (b-1) \sum_{i=1}^{n} \ln(1 - x_i^a) \] The Gradient Vector with respect to \(a\): \[ \ell_a^\prime(a, b) = \frac{n}{a} + \sum_{i=1}^{n} \ln(x_i) - (b-1) \sum_{i=1}^{n} \frac{x_i^a \ln(x_i)}{1 - x_i^a} \] Partial derivative with respect to \(b\): \[ \ell_b^\prime(a, b) = \frac{n}{b} + \sum_{i=1}^{n} \ln(1 - x_i^a) \]

Problem A3:

Based on the gradients functions obtained in the above problem A2, derive the observed Fisher Information matrix (i.e, the negative Hessian Matrix).

Answer:

From Problem A2, the gradient functions are: \[ \ell_a^\prime(a, b) = \frac{n}{a} + \sum_{i=1}^{n} \ln(x_i) - (b-1) \sum_{i=1}^{n} \frac{x_i^a \ln(x_i)}{1 - x_i^a} \] and \[ \ell_b^\prime(a, b) = \frac{n}{b} + \sum_{i=1}^{n} \ln(1 - x_i^a). \] Second derivative with respect to \(a\) \[ \frac{\partial^2 \ell}{\partial a^2} = -\frac{n}{a^2} - (b-1) \sum_{i=1}^{n} \frac{x_i^a (\ln x_i)^2}{(1-x_i^a)^2} \] Second derivative with respect to \(b\) \[ \frac{\partial^2 \ell}{\partial b^2} = -\frac{n}{b^2} \] Mixed Observed Information: \[ \frac{\partial^2 \ell}{\partial a \partial b} = -\sum_{i=1}^{n} \frac{x_i^a \ln x_i}{1-x_i^a} \] Therefore, the Hseesian matrix is: \[ H(a,b) = \begin{bmatrix} -\frac{n}{a^2} - (b-1) \sum_{i=1}^{n} \frac{x_i^a (\ln x_i)^2}{(1-x_i^a)^2} & -\sum_{i=1}^{n} \frac{x_i^a \ln x_i}{1-x_i^a} \\ -\sum_{i=1}^{n} \frac{x_i^a \ln x_i}{1-x_i^a} & -\frac{n}{b^2} \end{bmatrix} \] The observed Fisher Information is defined as the negative Hessian matrix: \[ J(a,b) = -H(a,b) \]

Problem A4:

Consider power distribution \(F(x) = x^a, (a >0 \quad \text{ and }\quad x \in (0,1))\), a special case of the Kumaraswamy distribution with \(b = 1\), and a random sample from this distribution \(\{ x_1, x_2, \cdots, x_n\}\). Derive the MLE and MME of \(a\) respectively. [Hint: To find the MME, you need to compute the moment of the power distribution; that is, \(E[X^k] = \int_0^1 x^k F'(x) dx\). Note that both the MLE and the MME have closed-form expressions.]

Answer:

Maximum Likelihood Estimator (MLE) differenitate with respect to \(a\): \[ \frac{d\ell}{da} = \frac{n}{a} + \sum_{i=1}^n \ln x_i = 0 \] The MLE \(\hat{a}\) is: \[ \hat{a} = \frac{-n}{\sum_{i=1}^n \ln x_i} \] The MME \(a\):

\[ E[X] = \int_{0}^{1} x f(x) dx = \int_{0}^{1} x(ax^{a-1}) dx = a \int_{0}^{1} x^a dx=\frac{a}{a+1} \] we set the theoretical mean equal to the sample mean \(\bar{x}\): \[ \bar{x} = \frac{a}{a+1} \] Solving for \(a\): \[ \bar{x}(a + 1) = a \] \[ \bar{x}a + \bar{x} = a \] \[ \bar{x} = a - \bar{x}a \] \[ \bar{x} = a(1 - \bar{x}) \] The MME is therefore: \[ \tilde{a} = \frac{\bar{x}}{1 - \bar{x}} \]

Problem A5:

Using the same setting as in Problem A4, find the asymptotic (Wald) confidence interval for \(a\). [Hint: Compute the Fisher information for \(a\), then take its reciprocal to obtain the variance.]

Answer:

Find the second derivative of the log-likelihood derived in Problem A4: \[ \ell(a) = n \ln(a) + (a-1) \sum_{i=1}^{n} \ln(x_i) \] The second derivative of the log-likelihood \(\ell(a)\) with respect to \(a\) is: \[ \frac{d^2 \ell}{d a^2} = -\frac{n}{a^2} \] So, the Fisher information for the sample is the negative expectation of the second derivative: \[ I_n(a) = -E\left[ \frac{d^2 \ell}{d a^2} \right] = \frac{n}{a^2} \] Therefore, \[ \text{Var}(\hat{a}) \approx \frac{a^2}{n} \] Using the plug-in principle, we estimate this variance using the MLE \(\hat{a}\): \[ \widehat{\text{Var}}(\hat{a}) = \frac{\hat{a}^2}{n} \] Thus the asymptotic \((1 - \alpha) \times 100\%\) Wald confidence interval for \(a\) is: \[ \hat{a} \pm z_{1-\alpha/2} \frac{\hat{a}}{\sqrt{n}} \] For a 95% confidence interval (\(z_{0.975} \approx 1.96\)): \[ \left[ \hat{a} \left( 1 - \frac{1.96}{\sqrt{n}} \right), \quad \hat{a} \left( 1 + \frac{1.96}{\sqrt{n}} \right) \right] \]

Problem A6:

Using the same setting as in Problem A4, perform a likelihood ratio test for the hypothesis \(H_0 :a=1\) (i.e., the power distribution reduces to a uniform distribution). [Hint: Evaluate the log-likelihood function at the maximum likelihood estimate \(\hat{a}\) and at \(a=1\), then use these values to construct the LRT test statistic.]

Answer:

We test: \[ H_0: a = 1 \quad \text{vs.} \quad H_a: a \neq 1 \] The likelihood ratio test statistic is: \[ LR = 2 \left[ \ell(\hat{a}) - \ell(1) \right] \] Under the power model: \[ \ell(a) = n \ln(a) + (a-1) \sum_{i=1}^{n} \ln(x_i) \]

Under the Alternative Hypothesis (\(H_a: a = \hat{a}\)): \[\ell(H_a) = \ell(\hat{a}) = n \ln(\hat{a}) + (\hat{a}-1) \sum_{i=1}^{n} \ln(x_i)\] \[\ell(1) = n \ln(1) + (1-1) \sum_{i=1}^{n} \ln(x_i)=0\] Hence, \[ LR = 2 [ n \ln \hat{a} + (\hat{a}-1) \sum_{i=1}^{n} \ln(x_i)] \] Under \(H_0\), for large \(n\): \[LR \sim \chi^2_1\] Reject \(H_0\) at the 5% significance level if: \[LR > 3.841\]

Part B: Numerical Analysis

All code must be well commented and adhere to best coding practices

Working Dataset: A small reservoir supplies water to a town. During the dry season (50 days), engineers record the fraction of usable storage filled each morning. Values near 0 mean the reservoir is nearly empty; values near 1 mean it’s full. The distribution tends to be right‑skewed (mostly low levels due to drought) but with occasional replenishment.

The following 50 data points (ordered for clarity) represent the daily proportion of usable storage:

0.12, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.20, 0.21, 0.22,
0.23, 0.24, 0.25, 0.26, 0.27, 0.28, 0.29, 0.30, 0.31, 0.32,
0.33, 0.34, 0.35, 0.36, 0.37, 0.38, 0.39, 0.40, 0.41, 0.42,
0.43, 0.44, 0.45, 0.46, 0.47, 0.48, 0.49, 0.50, 0.51, 0.52,
0.53, 0.54, 0.55, 0.56, 0.57, 0.58, 0.59, 0.60, 0.61, 0.78

Problem B1:

Fit the Kumaraswamy distribution to the above data. Use the derivations in Problem A2 to find the MLE of \(a\) and \(b\). Please copy the key formulas before coding.

Key formulas:

The log-likelihood is \[\ell(a,b) = n \ln a + n \ln b + (a-1) \sum_{i=1}^{n} \ln x_i + (b-1) \sum_{i=1}^{n} \ln(1 - x_i^a)\] The gradient is \[\ell_a^\prime(a,b) = \frac{n}{a} + \sum_{i=1}^{n} \ln x_i - (b-1) \sum_{i=1}^{n} \frac{x_i^a \ln x_i}{1 - x_i^a}\] \[\ell_b^\prime(a,b) = \frac{n}{b} + \sum_{i=1}^{n} \ln(1 - x_i^a)\]

# data
x <- c(
  0.12, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.20, 0.21, 0.22,
  0.23, 0.24, 0.25, 0.26, 0.27, 0.28, 0.29, 0.30, 0.31, 0.32,
  0.33, 0.34, 0.35, 0.36, 0.37, 0.38, 0.39, 0.40, 0.41, 0.42,
  0.43, 0.44, 0.45, 0.46, 0.47, 0.48, 0.49, 0.50, 0.51, 0.52,
  0.53, 0.54, 0.55, 0.56, 0.57, 0.58, 0.59, 0.60, 0.61, 0.78
)

n <- length(x)
xbar <- mean(x)
sumlogx <- sum(log(x))


# Log-likelihood function for Kumaraswamy distribution
kuma.loglik <- function(params, data){
  a <- params[1]
  b <- params[2]
  
  # Keep parameters positive
  if(a <= 0 || b <= 0){
    return(-1e10)
  }
  
  n <- length(data)
  
  # Log-likelihood
  ll <- n * log(a) + n * log(b) +
        (a - 1) * sum(log(data)) +
        (b - 1) * sum(log(1 - data^a))
  
  return(ll)
}

# Gradient (score vector) of the log-likelihood
kuma.grad <- function(params, data){
  a <- params[1]
  b <- params[2]
  
  n <- length(data)
  
  # Partial derivative with respect to a
  da <- n / a + sum(log(data)) -
        (b - 1) * sum((data^a * log(data)) / (1 - data^a))
  
  # Partial derivative with respect to b
  db <- n / b + sum(log(1 - data^a))
  
  return(c(da, db))
}

# MLE using optim() with gradient, following lecture style
mle_result <- optim(
  par = c(1, 1),
  fn = kuma.loglik,
  gr = kuma.grad,
  data = x,
  method = "L-BFGS-B",
  hessian = TRUE,
  lower = c(1e-8, 1e-8),
  control = list(
    maxit = 1000,
    fnscale = -1,
    trace = FALSE,
    abstol = 1e-8
  )
)

# Extract MLEs
a_hat <- mle_result$par[1]
b_hat <- mle_result$par[2]

# One clean output box
data.frame(
  a_hat = a_hat,
  b_hat = b_hat
)

     a_hat    b_hat
1 2.529601 7.883389

Conclusion:

The fitted Kumaraswamy distribution gives MLEs approximately \(\hat{a} \approx 2.53 \quad \text{and} \quad \hat{b} \approx 7.88\).

Problem B2:

Fit the power distribution to the above data using the derived of \(a\) obtained in Problem A4 to test the following hypothesis using likelihood ratio procedure ar significance level \(\alpha = 0.05\):

\[ H_0: b = 1 \quad \text{ versus } \quad H_a: b \ne 1. \]

State the statistical decision clearly. What is the practical implication of the testing result?

# Power distribution MLE of a from Part A4
a_hat_power <- -n / sum(log(x))

# Log-likelihood under the power model (b = 1)
power.loglik <- function(a, data){
  n <- length(data)
  ll <- n * log(a) + (a - 1) * sum(log(data))
  return(ll)
}

# Restricted and unrestricted log-likelihood values
ll_full <- kuma.loglik(c(a_hat, b_hat), x)
ll_restricted <- power.loglik(a_hat_power, x)

# Likelihood ratio test statistic
LR_b1 <- 2 * (ll_full - ll_restricted)

# p-value from chi-square(1)
p_value_b1 <- 1 - pchisq(LR_b1, df = 1)

data.frame(
  a_hat_power = a_hat_power,
  ll_full = ll_full,
  ll_restricted = ll_restricted,
  LR_stat = LR_b1,
  p_value = p_value_b1
)

  a_hat_power  ll_full ll_restricted  LR_stat      p_value
1   0.9393197 24.56271     0.1000435 48.92533 2.658984e-12

Conclusion:

Since the p-value is extremely small, we reject \(H_0: b = 1\). This means the power distribution is too restrictive for these data. In practical terms, the reservoir storage proportions are not adequately described by the one-parameter power model, and the second Kumaraswamy parameter b is needed for a better fit.

Problem B3:

Use the procedure and code from Problem B1 to estimate the MLEs of \(a\) and \(b\), and then complete the following analyses:

(1). Obtain the bootstrap sampling distributions of \(\hat{a}\) and \(\hat{b}\) and plot each distribution using Gaussian kernel density curves.

# Bootstrap sampling distributions of a_hat and b_hat
set.seed(123)

# Number of bootstrap samples
B <- 2000

# Store bootstrap estimates
boot_a <- NULL
boot_b <- NULL

for(i in 1:B){
  
  # Bootstrap sample
  boot_sample <- sample(x, size = n, replace = TRUE)
  
  # Refit Kumaraswamy model to bootstrap sample
  boot_fit <- optim(
    par = c(a_hat, b_hat),
    fn = kuma.loglik,
    gr = kuma.grad,
    data = boot_sample,
    method = "L-BFGS-B",
    lower = c(1e-8, 1e-8),
    control = list(
      maxit = 1000,
      fnscale = -1,
      trace = FALSE,
      abstol = 1e-8
    )
  )
  
  # Save bootstrap estimates
  boot_a[i] <- boot_fit$par[1]
  boot_b[i] <- boot_fit$par[2]
}

# Plot bootstrap distribution of a_hat
hist(boot_a,
     probability = TRUE,
     breaks = 20,
     main = "Bootstrap Distribution of a_hat",
     xlab = expression(hat(a)))
lines(density(boot_a), col = "blue", lwd = 2)

# Plot bootstrap distribution of b_hat
hist(boot_b,
     probability = TRUE,
     breaks = 20,
     main = "Bootstrap Distribution of b_hat",
     xlab = expression(hat(b)))
lines(density(boot_b), col = "blue", lwd = 2)

Conclusion:

The bootstrap sampling distribution of \(\hat{a}\) appears fairly concentrated and roughly symmetric around the MLE. In contrast, the bootstrap sampling distribution of \(\hat{b}\) is more spread out and shows noticeable right-skewness. The Gaussian kernel density curves confirm that the normal approximation is more reasonable for \(\hat{a}\) than for \(\hat{b}\).

(2). Construct both the \(95\%\) bootstrap confidence interval and the Wald confidence interval for \(b\). Do these intervals agree with the results obtained in Problem B2? [Compute the standard error of \(\hat{b}\) using the observed Fisher information matrix, i.e., the inverse of the negative Hessian obtained from optim()]

# 95% percentile bootstrap CI for b
ci_boot_b <- quantile(boot_b, probs = c(0.025, 0.975))

# Wald CI for b using observed Fisher information
Hess <- mle_result$hessian
covar <- solve(-Hess)
se_b_hat <- sqrt(covar[2, 2])

z <- qnorm(0.975)
ci_wald_b <- c(
  b_hat - z * se_b_hat,
  b_hat + z * se_b_hat
)

data.frame(
  lower_boot = ci_boot_b[1],
  upper_boot = ci_boot_b[2],
  lower_wald = ci_wald_b[1],
  upper_wald = ci_wald_b[2]
)

     lower_boot upper_boot lower_wald upper_wald
2.5%    5.03667   16.34823   3.484028   12.28275

Answer:

Yes. Both the 95% percentile bootstrap confidence interval and the Wald confidence interval for b lie entirely above 1. Therefore, both intervals agree with the likelihood ratio test result in Problem B2, which rejects \(H_0: b=1\).

(3). Based on the bootstrap sampling distributions from part (1) of this problem, assess whether the validity of the Wald confidence interval is supported.

Answer:

The Wald confidence interval is based on asymptotic normality, so it is better supported when the sampling distribution of the estimator is approximately symmetric. Here, the bootstrap distribution of \(\hat{b}\) shows noticeable right-skewness. This suggests that the normal approximation is not perfect, so the Wald confidence interval is less strongly supported than the bootstrap percentile interval. Therefore, the bootstrap interval is more appropriate for b in this analysis.

Problem B4:

In the introduction to the working model for this exam, the Kumaraswamy distribution reduces to the uniform distribution on (0,1). In this problem, we perform a likelihood ratio test for the following hypothesis to assess whether the data come from the uniform distribution on (0,1):

\[ H_0: a = 1\quad \& \quad b = 1\quad \text{ versus } \quad H_a: a \ne 1 \quad \text{or} \quad b \ne 1 \quad \text{or}\quad (a \ne 1 \quad \& \quad b \ne 1). \]

Provide a practical interpretation of the above test result. [Hint: \(H_a\) basically says that there is no constraints for \(a\) and \(b\). Please review the lecture note for module 11 on the likelihood ratio test before coding.]

# Log-likelihood under the null model a = 1, b = 1 (uniform distribution)
ll_null <- kuma.loglik(c(1, 1), x)

# Likelihood ratio statistic
LR_uniform <- 2 * (ll_full - ll_null)

# p-value from chi-square(2)
p_value_uniform <- 1 - pchisq(LR_uniform, df = 2)

data.frame(
  ll_full = ll_full,
  ll_null = ll_null,
  LR_stat = LR_uniform,
  p_value = p_value_uniform
)

   ll_full ll_null  LR_stat      p_value
1 24.56271       0 49.12542 2.150558e-11

Conclusion:

The p-value is extremely small, so we reject \(H_0: a=1, b = 1\). This means the reservoir storage proportions are not consistent with a uniform distribution on (0,1). In practical terms, the data show a clear distributional shape instead of being spread evenly across all possible storage levels. Therefore, the full Kumaraswamy model is more appropriate than the uniform model.