The objective of this assignment is to build and evaluate the most powerful likelihood ratio test for deciding whether one observation comes from the null mixture distribution or from a Poisson(20) alternative. It also aims to compare that test to a simpler tail-based test to show how choosing the rejection region by likelihood ratios can give much higher power.
# Assignment 6: Likelihood Ratio Test, 1 sample
# H0 = 0.5*Poisson(10) + 0.5*Poisson(40)
# HA = Poisson(20)
# Step 1: Set up observed values and probabilities
x <- 0:100
# HA probabilities using Poisson(20)
p_alt <- dpois(x, 20)
# H0 probabilities using mixture of Poisson(10) and Poisson(40)
p_null <- (dpois(x, 10) + dpois(x, 40)) / 2
# Compute likelihood ratios
LR <- p_alt / p_null
# Store everything in a data frame for easy viewing
results <- data.frame(
x = x,
p_alt = p_alt,
p_null = p_null,
LR = LR
)
# Step 2: Order likelihood ratios from largest to smallest
ordered_results <- results[order(-results$LR), ]
# View the ordered table
ordered_results## x p_alt p_null LR
## 23 22 7.691369e-02 5.344203e-04 1.439199e+02
## 22 21 8.460506e-02 6.271600e-04 1.349019e+02
## 24 23 6.688147e-02 6.660047e-04 1.004219e+02
## 21 20 8.883532e-02 1.029039e-03 8.632839e+01
## 25 24 5.573456e-02 1.000249e-03 5.572071e+01
## 20 19 8.883532e-02 1.914081e-03 4.641148e+01
## 26 25 4.458765e-02 1.556494e-03 2.864620e+01
## 19 18 8.439355e-02 3.568354e-03 2.365055e+01
## 27 26 3.429819e-02 2.377720e-03 1.442482e+01
## 18 17 7.595420e-02 6.392258e-03 1.188222e+01
## 28 27 2.540607e-02 3.516295e-03 7.225239e+00
## 17 16 6.456107e-02 1.085376e-02 5.948269e+00
## 29 28 1.814719e-02 5.021044e-03 3.614226e+00
## 16 15 5.164885e-02 1.736078e-02 2.975031e+00
## 30 29 1.251530e-02 6.924808e-03 1.807314e+00
## 15 14 3.873664e-02 2.603921e-02 1.487628e+00
## 31 30 8.343536e-03 9.232821e-03 9.036822e-01
## 14 13 2.711565e-02 3.645420e-02 7.438278e-01
## 32 31 5.382927e-03 1.191323e-02 4.518443e-01
## 13 12 1.762517e-02 4.739024e-02 3.719156e-01
## 33 32 3.364329e-03 1.489152e-02 2.259225e-01
## 12 11 1.057510e-02 5.686822e-02 1.859580e-01
## 34 33 2.038987e-03 1.805032e-02 1.129613e-01
## 11 10 5.816307e-03 6.255502e-02 9.297905e-02
## 35 34 1.199404e-03 2.123566e-02 5.648066e-02
## 10 9 2.908153e-03 6.255502e-02 4.648953e-02
## 36 35 6.853739e-04 2.426933e-02 2.824033e-02
## 9 8 1.308669e-03 5.629952e-02 2.324476e-02
## 37 36 3.807633e-04 2.696592e-02 1.412017e-02
## 8 7 5.234676e-04 4.503961e-02 1.162238e-02
## 38 37 2.058180e-04 2.915235e-02 7.060083e-03
## 7 6 1.832137e-04 3.152773e-02 5.811191e-03
## 39 38 1.083253e-04 3.068668e-02 3.530041e-03
## 6 5 5.496410e-05 1.891664e-02 2.905596e-03
## 40 39 5.555141e-05 3.147352e-02 1.765021e-03
## 5 4 1.374102e-05 9.458319e-03 1.452798e-03
## 41 40 2.777571e-05 3.147352e-02 8.825103e-04
## 4 3 2.748205e-06 3.783327e-03 7.263989e-04
## 42 41 1.354913e-05 3.070587e-02 4.412552e-04
## 3 2 4.122307e-07 1.134998e-03 3.631994e-04
## 43 42 6.451964e-06 2.924369e-02 2.206276e-04
## 2 1 4.122307e-08 2.269996e-04 1.815997e-04
## 44 43 3.000914e-06 2.720343e-02 1.103138e-04
## 1 0 2.061154e-09 2.269996e-05 9.079986e-05
## 45 44 1.364052e-06 2.473039e-02 5.515690e-05
## 46 45 6.062452e-07 2.198257e-02 2.757845e-05
## 47 46 2.635849e-07 1.911528e-02 1.378922e-05
## 48 47 1.121638e-07 1.626832e-02 6.894612e-06
## 49 48 4.673491e-08 1.355694e-02 3.447306e-06
## 50 49 1.907547e-08 1.106689e-02 1.723653e-06
## 51 50 7.630189e-09 8.853509e-03 8.618265e-07
## 52 51 2.992231e-09 6.943928e-03 4.309133e-07
## 53 52 1.150858e-09 5.341483e-03 2.154566e-07
## 54 53 4.342860e-10 4.031308e-03 1.077283e-07
## 55 54 1.608467e-10 2.986154e-03 5.386416e-08
## 56 55 5.848970e-11 2.171749e-03 2.693208e-08
## 57 56 2.088918e-11 1.551249e-03 1.346604e-08
## 58 57 7.329537e-12 1.088596e-03 6.733020e-09
## 59 58 2.527426e-12 7.507557e-04 3.366510e-09
## 60 59 8.567547e-13 5.089869e-04 1.683255e-09
## 61 60 2.855849e-13 3.393246e-04 8.416275e-10
## 62 61 9.363440e-14 2.225079e-04 4.208137e-10
## 63 62 3.020464e-14 1.435535e-04 2.104069e-10
## 64 63 9.588776e-15 9.114509e-05 1.052034e-10
## 65 64 2.996492e-15 5.696568e-05 5.260172e-11
## 66 65 9.219977e-16 3.505580e-05 2.630086e-11
## 67 66 2.793932e-16 2.124594e-05 1.315043e-11
## 68 67 8.340097e-17 1.268414e-05 6.575214e-12
## 69 68 2.452970e-17 7.461261e-06 3.287607e-12
## 70 69 7.110057e-18 4.325369e-06 1.643804e-12
## 71 70 2.031445e-18 2.471639e-06 8.219018e-13
## 72 71 5.722380e-19 1.392473e-06 4.109509e-13
## 73 72 1.589550e-19 7.735960e-07 2.054755e-13
## 74 73 4.354931e-20 4.238882e-07 1.027377e-13
## 75 74 1.177008e-20 2.291288e-07 5.136886e-14
## 76 75 3.138689e-21 1.222020e-07 2.568443e-14
## 77 76 8.259708e-22 6.431685e-08 1.284222e-14
## 78 77 2.145379e-22 3.341135e-08 6.421108e-15
## 79 78 5.500971e-23 1.713403e-08 3.210554e-15
## 80 79 1.392651e-23 8.675456e-09 1.605277e-15
## 81 80 3.481627e-24 4.337728e-09 8.026385e-16
## 82 81 8.596611e-25 2.142088e-09 4.013192e-16
## 83 82 2.096734e-25 1.044921e-09 2.006596e-16
## 84 83 5.052372e-26 5.035763e-10 1.003298e-16
## 85 84 1.202946e-26 2.397983e-10 5.016491e-17
## 86 85 2.830460e-27 1.128462e-10 2.508245e-17
## 87 86 6.582466e-28 5.248662e-11 1.254123e-17
## 88 87 1.513211e-28 2.413178e-11 6.270613e-18
## 89 88 3.439115e-29 1.096899e-11 3.135307e-18
## 90 89 7.728348e-30 4.929884e-12 1.567653e-18
## 91 90 1.717411e-30 2.191059e-12 7.838266e-19
## 92 91 3.774529e-31 9.631030e-13 3.919133e-19
## 93 92 8.205498e-32 4.187405e-13 1.959567e-19
## 94 93 1.764623e-32 1.801034e-13 9.797833e-20
## 95 94 3.754518e-33 7.663975e-14 4.898917e-20
## 96 95 7.904248e-34 3.226937e-14 2.449458e-20
## 97 96 1.646718e-34 1.344557e-14 1.224729e-20
## 98 97 3.395295e-35 5.544565e-15 6.123646e-21
## 99 98 6.929174e-36 2.263088e-15 3.061823e-21
## 100 99 1.399833e-36 9.143789e-16 1.530911e-21
## 101 100 2.799666e-37 3.657516e-16 7.654557e-22# Step 3: Find rejection region for alpha about 0.05
# Add cumulative null probability as I move from largest LR downward
ordered_results$cum_null_prob <- cumsum(ordered_results$p_null)
# Show rows where cumulative null prob is still under 0.05
# Using < instead of <= to avoid floating point issues at the boundary
ordered_results[ordered_results$cum_null_prob < 0.05, ]## x p_alt p_null LR cum_null_prob
## 23 22 0.07691369 0.0005344203 143.919852 0.0005344203
## 22 21 0.08460506 0.0006271600 134.901884 0.0011615803
## 24 23 0.06688147 0.0006660047 100.421927 0.0018275850
## 21 20 0.08883532 0.0010290395 86.328385 0.0028566245
## 25 24 0.05573456 0.0010002486 55.720707 0.0038568731
## 20 19 0.08883532 0.0019140807 46.411479 0.0057709539
## 26 25 0.04458765 0.0015564942 28.646204 0.0073274480
## 19 18 0.08439355 0.0035683542 23.650553 0.0108958022
## 27 26 0.03429819 0.0023777204 14.424821 0.0132735226
## 18 17 0.07595420 0.0063922580 11.882217 0.0196657806
## 28 27 0.02540607 0.0035162946 7.225239 0.0231820752
## 17 16 0.06456107 0.0108537572 5.948269 0.0340358324
## 29 28 0.01814719 0.0050210445 3.614226 0.0390568769# Also show the next row that pushes it just over 0.05
ordered_results[1:(sum(ordered_results$cum_null_prob < 0.05) + 1), ]## x p_alt p_null LR cum_null_prob
## 23 22 0.07691369 0.0005344203 143.919852 0.0005344203
## 22 21 0.08460506 0.0006271600 134.901884 0.0011615803
## 24 23 0.06688147 0.0006660047 100.421927 0.0018275850
## 21 20 0.08883532 0.0010290395 86.328385 0.0028566245
## 25 24 0.05573456 0.0010002486 55.720707 0.0038568731
## 20 19 0.08883532 0.0019140807 46.411479 0.0057709539
## 26 25 0.04458765 0.0015564942 28.646204 0.0073274480
## 19 18 0.08439355 0.0035683542 23.650553 0.0108958022
## 27 26 0.03429819 0.0023777204 14.424821 0.0132735226
## 18 17 0.07595420 0.0063922580 11.882217 0.0196657806
## 28 27 0.02540607 0.0035162946 7.225239 0.0231820752
## 17 16 0.06456107 0.0108537572 5.948269 0.0340358324
## 29 28 0.01814719 0.0050210445 3.614226 0.0390568769
## 16 15 0.05164885 0.0173607790 2.975031 0.0564176559# These x values make up my rejection region
reject_x <- ordered_results$x[ordered_results$cum_null_prob < 0.05]
reject_x## [1] 22 21 23 20 24 19 25 18 26 17 27 16 28# Check the actual alpha for this rejection region
alpha_actual <- sum(p_null[x %in% reject_x])
alpha_actual## [1] 0.03905688# Also check what happens if I include the boundary-crossing point
reject_x2 <- ordered_results$x[1:(sum(ordered_results$cum_null_prob < 0.05) + 1)]
reject_x2## [1] 22 21 23 20 24 19 25 18 26 17 27 16 28 15alpha_actual2 <- sum(p_null[x %in% reject_x2])
alpha_actual2## [1] 0.05641766# Step 4: How much error from stopping at 100?
# The H0 mass above x=100 is whatever is left after summing p_null over 0:100
# Poisson(40) at x=100 is about 9.5 SDs above its mean, so I expect this to be tiny
error_above_100 <- 1 - sum(p_null)
cat("Truncation error (H0 mass above x=100):", formatC(error_above_100, format = "e", digits = 3), "-- negligible\n")## Truncation error (H0 mass above x=100): 2.220e-16 -- negligible# Step 5: Power of the likelihood ratio test
# Power is the probability I reject under HA = Poisson(20)
power_LRT <- sum(p_alt[x %in% reject_x])
power_LRT## [1] 0.8091533# Power if I use the slightly larger rejection region
power_LRT2 <- sum(p_alt[x %in% reject_x2])
power_LRT2## [1] 0.8608022# Step 6: Compare to test with alpha = 0.025 in each tail
# Lower tail: largest c1 where P(X <= c1 | Poisson(10)) is still <= 0.025
c1 <- max(x[ppois(x, 10) <= 0.025])
c1## [1] 3# Upper tail: smallest c2 where P(X >= c2 | Poisson(40)) is still <= 0.025
# P(X >= c2) = 1 - P(X <= c2-1)
c2 <- min(x[(1 - ppois(x - 1, 40)) <= 0.025])
c2## [1] 54# Rejection region for this comparison test
reject_compare <- x[x <= c1 | x >= c2]
reject_compare## [1] 0 1 2 3 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68
## [20] 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87
## [39] 88 89 90 91 92 93 94 95 96 97 98 99 100# Actual size of the comparison test under my H0 mixture
# This ends up smaller than 0.05 because the mixture spreads probability
# away from the extreme tails -- so the LRT uses its alpha budget better,
# which is why it has higher power
alpha_compare <- sum(p_null[x %in% reject_compare])
alpha_compare## [1] 0.01516533# Power of the comparison test under HA = Poisson(20)
power_compare <- sum(p_alt[x %in% reject_compare])
power_compare## [1] 3.203971e-06# Step 7: Print summary
cat("Likelihood ratio test rejection region:\n")## Likelihood ratio test rejection region:print(sort(reject_x))## [1] 16 17 18 19 20 21 22 23 24 25 26 27 28cat("\nActual alpha of LRT:\n")##
## Actual alpha of LRT:print(alpha_actual)## [1] 0.03905688cat("\nPower of LRT:\n")##
## Power of LRT:print(power_LRT)## [1] 0.8091533cat("\nError from truncating at 100 under H0:\n")##
## Error from truncating at 100 under H0:print(error_above_100)## [1] 2.220446e-16cat("\nComparison test rejection region:\n")##
## Comparison test rejection region:print(reject_compare)## [1] 0 1 2 3 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68
## [20] 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87
## [39] 88 89 90 91 92 93 94 95 96 97 98 99 100cat("\nActual alpha of comparison test:\n")##
## Actual alpha of comparison test:print(alpha_compare)## [1] 0.01516533cat("\nPower of comparison test:\n")##
## Power of comparison test:print(power_compare)## [1] 3.203971e-06The likelihood ratio test rejects when 16 <= x <= 28. Its actual significance level is about 0.0391, which is below the required 0.05, so it is a valid level-0.05 test. Its power is about 0.8092, meaning it correctly rejects the null about 81% of the time when the true distribution is Poisson(20). This is a strong test because values in the middle are much more likely under Poisson(20) than under the null mixture of Poisson(10) and Poisson(40).
The comparison test rejects when x<= 3 or x >= 54. Its actual significance level is only about 0.0152, so it does not even use most of the available 0.05 error rate. Its power is essentially 0.0000032, which is almost zero. This means it almost never rejects when the alternative is true. The reason is that very small values are consistent with Poisson(10), and very large values are consistent with Poisson(40), so these tails do not provide strong evidence against the null mixture.
The probability under the null above X =100 is about 2.22 x 10-16, which is negligible. So stopping the calculations at 100 causes essentially no error.
Overall the LRT is clearly the better test. It places the rejection region where the alternative distribution is most favored relative to the null, giving high power while still keeping the significance level below 0.05. The tail-based comparison test is very poor because it rejects in regions that are actually more compatible with the null than with the alternative.