Finding P-value and power of the Most Powerful Test

Question

You observe a sample $X_1, \quad, X_{20}$ with the density $$ f(x, \vartheta)=2\left(x / \vartheta^2\right) I_{[0 \leq x<\vartheta]} $$ with an unknown parameter $\vartheta>0$, yielding $$ \min \left(X_i, i=1, ..., 20\right)=0.035, \quad \max \left(X_i, 1 \leq i \leq 20\right)=3.66 $$ Based on these data, give numerical expressions for

(a) the p-value of the most powerful hypothesis test for $H_0: \vartheta=4.0$ versus $H_A: \vartheta=3.66$

(b) the power of the most powerful hypothesis test of significance level $0.05$ of $H_0: \vartheta=4.0$ versis $H_B: \vartheta \leq 3.75$, at the point alternative $\vartheta=3.70$

$\textbf{My attempt:}$

$\textbf{(a)}$ Using NP-Lemma, the best test is given by $\frac{L_0(\theta \mid X)}{L_1(\theta\mid X)} < k$, so $$\frac{L_0(\theta \mid X)}{L_1(\theta\mid X)} = \frac{(x) I_{[X_{(n)} < 4]}}{(x) I_{[X_{(n)} < 3.66]}} < k \left(\frac{4}{3.66}\right)^2 = k_1 $$

Now I can see how the best test would be reject $H_0$ if $X_{(n)} < k_2$ where $k_2$ is a constant satisfying $P(X_{(n)} < k_2\mid \theta = 4) = \alpha$, but I am struggling to see how I got from that ratio of indicator functions to this. When the sufficient statistic is something outside of the indicator function, such as the sample mean, it is clear how we get the best test since we can try to isolate the sample mean on the left-hand side. In this scenario, I'm not exactly sure how we can isolate $X_{(n)}$ within the indicator function.

Moving on to find the p-value we plug in the maximum for $k_2$: $$P(X_{(n)} < 3.66\mid\theta = 4) = (P(X_1 < 3.66\mid \theta = 4))^n = (3.66/4)^{40} = 0.029$$

$\textbf{b)}$ I first show that the sufficient statistic is $X_{(n)}$, and then that this is a monotone likelihood ratio family. This allows us to use Karlin-Rubin Theorem which states that the best test will be: Reject $H_0$ if $X_{(n)} < c$ where $c$ is a constant satisfying $$P(X_{(n)} < c\mid \theta = 4) = 0.05 \implies c = 0.05^{1/40}\times 4$$

$$\textrm{Power} = P(\textrm{Reject }~H_0\mid \textrm{point alternative}) = P(X_{(n)} < c\mid \theta = 3.7) = \frac{0.05\times 4^{40}}{3.7^{40}} = 1.13$$

Now this clearly does not make sense because power is between 0 and 1. However, if we did

$$\textrm{Power} = P(\textrm{Reject }~H_0\mid H_1) = P(X_{(n)} < c\mid \theta = 3.75) = \frac{0.05\times 4^{40}}{3.75^{40}} = 0.66 $$

this does make sense and in general we do use $H_1$, but then I would ask what is the point of the point alternative? Unless this is a mistake in the question or I am conceptually not understanding something?

Answer 1

The indicator function sets the limits of integration, but does not have to be carried through since the alternate hypothesis fits within the range of null hypothesis and the maximum sample fits within the range of both.

$$f(x|\theta) = \frac{2x}{\theta^2}$$

$$\int_0^{\theta} \frac{2x}{\theta^2} dx = 1$$

$$F(x|\theta) = \frac{x^2}{\theta^2}$$

Let $X_{(1)}$ and $X_{(n)}$ be the first and n-th order statistics.

$$f_{X_{(1)}}(x|\theta, n) = n\left(1-\frac{x^2}{\theta^2}\right)^{n-1}\left(\frac{2x}{\theta^2}\right)$$

$$F_{X_{(1)}}(x|\theta, n) = 1 - \left(1 - \frac{x^2}{\theta^2}\right)^n$$

$$f_{X_{(n)}}(x|\theta, n) = n\left(\frac{x^2}{\theta^2}\right)^{n-1}\left(\frac{2x}{\theta^2}\right)$$

Joint density ref

$$F_{X_{(1)}, X_{(n)}}(x_{(1)}, x_{(n)}) = [F(x_{(n)})]^n - [F(x_{(n)}) - F(x_{(1)})]^n$$

$$f_{X_{(1)}, X_{(n)}}(x_{(1)}, x_{(n)}) = n(n-1)[F(x_{(n)}) - F(x_{(1)})]^{(n-2)}f(x_{(1)})f(x_{(n)})$$

Now, the likelihood ratio:

The order statistics are independent samples in the limit as the sample gets large. If this sample is large enough, assuming independence:

$$\frac{L(\theta_0|X)}{L(\theta_1|X)} = \frac{f_{X_{(1)}}(x_{(1)}|\theta_0, n)f_{X_{(n)}}(x_{(n)}|\theta_0, n)}{f_{X_{(1)}}(x_{(1)}|\theta_1, n)f_{X_{(n)}}(x_{(n)}|\theta_1, n)} = \left(\frac{\theta_1}{\theta_0}\right)^{4n}\left(\frac{\theta_0^2-x_{(1)}^2}{\theta_1^2-x_{(1)}^2}\right)^{n-1} \le k$$

If we can't assume independence, then:

$$\frac{L(\theta_0|X)}{L(\theta_1|X)} = \frac{f_{X_{(1)}, X_{(n)}}(x_{(1)}, x_{(n)}|\theta_0)}{f_{X_{(1)}, X_{(n)}}(x_{(1)}, x_{(n)}|\theta_1)} = \left( \frac{\theta_1}{\theta_0} \right)^{2n}$$

The joint density formulation does not allow us to form a critical region, in the same way, but we can still answer the questions

(a) Under the independent min and max densities, the likelihood ratio is a function of $X_{(1)}$, and $k$ can be translated into a critical region on $X_{(1)}$, so the p-value is

$$P(X_{(1)} < x_{(1)}|\theta_0) = 1 - (1 - F(x|\theta_0))^n = 1 - \left(1 - \frac{0.035^2}{4^2}\right)^{20} = 0.00153$$

Under the joint density, we can think of a p-value as the region of the density with a more extreme min or max than the point estimates:

p-value = $P(X_{(1)} < x_{(1)}, X_{(n)} > x_{(n)} | \theta_0) = P(X_{(1)} < x_{(1)}) - P(X_{(1)} < x_{(1)}, X_{(n)} < x_{(n)})$

$$P(X_{(1)} < x_{(1)}, X_{(n)} > x_{(n)} | \theta_0) = F_{X_{(1)}}(0.035) - F_{X_{(1)}, X_{(n)}}(0.035, 3.66) = 0.00148$$

(b) Now, creating a critical region for a p-value of 0.05 under the approximation of independence of the min and max densities

$$P(X_{(1)} < x_{(1),0.05}|\theta_0) = 1 - (1 - F(x_{(1),0.05}|\theta_0))^n = 1 - \left(1 - \frac{x_{(1),0.05}^2}{4^2}\right)^{20} = 0.05$$

$$x_{(1),0.05} = 0.2024$$

At the boundary of the critical region,

$$\frac{L(\theta_0|X)}{L(\theta_1|X)} = \left(\frac{\theta_1}{\theta_0}\right)^{4n}\left(\frac{\theta_0^2-x_{(1),0.05}^2}{\theta_1^2-x_{(1),0.05}^2}\right)^{n-1} = \left(\frac{3.75}{4}\right)^{80}\left(\frac{4^2-0.2024^2}{3.75^2-0.2024^2}\right)^{19} = k = 0.06694$$

calculating power to reject under the alternate hypothesis:

$$P(X_{(1)} \le x_{(1),0.05} | \theta_1) = 1 - (1 - F(x_{(1),0.05}|\theta_1))^n = 1 - \left(1 - \frac{0.2024^2}{3.75^2}\right)^{20} = 0.0567$$

For the joint density, the simulation below will give an estimate of the power.

R code to show the calculations:

fn <- function(x, theta) 2*x/theta^2
Fn <- function(x, theta) x^2 / theta^2
Finv <- function(p, theta) sqrt(p*theta^2)

# check
integrate(fn, lower = 0, upper = 5, theta = 5)
#> 1 with absolute error < 1.1e-14
Fn(0, 5)
#> [1] 0
Fn(5, 5)
#> [1] 1
Finv(Fn(1.5, 5), 5)
#> [1] 1.5

# order statistics
fn_max <- function(x_max, theta, n) n*(x_max^2/theta^2)^(n-1)*2*x_max/theta^2
fn_min <- function(x_min, theta, n) n*(1-x_min^2/theta^2)^(n-1)*2*x_min/theta^2
Fn_min <- function(x_min, theta, n) 1 - (1 - x_min^2/theta^2)^n
joint_f <- function(x_min, x_max, theta, n) ifelse(x_max >= x_min, n*(n-1)*(Fn(x_max, theta) - Fn(x_min, theta))^(n-2)*fn(x_min, theta)*fn(x_max, theta), 0)
joint_F <- function(x_min, x_max, theta, n) ifelse(x_max >= x_min, (Fn(x_max, theta))^n - (Fn(x_max, theta) - Fn(x_min, theta))^n, 0)

# check
integrate(fn_max, lower = 0, upper = 5, theta = 5, n = 20)
#> 1 with absolute error < 2.8e-09
integrate(fn_min, lower = 0, upper = 5, theta = 5, n = 20)
#> 1 with absolute error < 3.5e-06
Fn_min(0, 5, 20)
#> [1] 0
Fn_min(5, 5, 20)
#> [1] 1
inner_integral <- function(x_max) {
  sapply(x_max, function(z) { # integrate out x_min
    integrate(joint_f, lower = 0, upper = z, x_max = z, theta = 5, n = 20)$value
  }) 
}
integrate(inner_integral, lower = 0, upper = 5)
#> 1 with absolute error < 2.8e-09
joint_F(0, 0, 5, 20)
#> [1] 0
joint_F(0, 5, 5, 20)
#> [1] 0
joint_F(5, 5, 5, 20)
#> [1] 1

# likelihood ratio test with independent min and max density
lrt <- function(x_min, x_max, theta_0, theta_1, n) fn_max(x_max, theta_0, n)*fn_min(x_min, theta_0, n)/fn_max(x_max, theta_1, n)/fn_min(x_min, theta_1, n)
lrt(0.035, 3.66, 4, 3.66, 20)
#> [1] 0.02397777

# likelihood ratio test with joint min and max density
lrt_joint <- function(x_min, x_max, theta_0, theta_1, n) joint_f(x_min, x_max, theta_0, n) / joint_f(x_min, x_max, theta_1, n)
lrt_joint(0.035, 3.66, 4, 3.66, 20)
#> [1] 0.02863148

# a - p-value with independent min and max density
Fn_min(0.035, 4, 20)
#> [1] 0.001530137

# a - p-value with joint min and max density
Fn_min(0.035, 4, 20) - joint_F(0.035, 3.66, 4, 20)
#> [1] 0.001477816

# find critical region with independent min and max density
temp <- function(x_min) (Fn_min(x_min, 4, 20) - 0.05)^2
o <- optimize(f = temp, interval = c(0, 4))
o$minimum
#> [1] 0.2024429
Fn_min(o$minimum, 4, 20)
#> [1] 0.05000127

# b - power with independent min and max density
Fn_min(o$minimum, 3.75, 20)
#> [1] 0.05670122

# b - power with joint min and max density
#    need to do this one by simulation only due to the rejection region

# simulation study

set.seed(10430938)

n <- 20
sims <- 100000
x_1s <- numeric(sims)
x_ns <- numeric(sims)
for (i in 1:sims) {
  p <- runif(n, 0, 1)
  x <- Finv(p, 4)
  x_1s[i] <- min(x)
  x_ns[i] <- max(x)
}
# p-value for independent min and max densities
length(which(x_1s < o$minimum)) / sims # should be equal to 0.05
#> [1] 0.04972
length(which(x_1s < 0.035)) / sims # p-value of part A
#> [1] 0.00148

# critical values for independent min and max densities
# part A
lrts <- lrt(x_1s, x_ns, 4, 3.66, 20)
length(which(lrts < (3.66 / 4)^(4*20)*((4^2-o$minimum^2) / (3.66^2 - o$minimum^2))^(19))) / sims
#> [1] 0.04972
# part B
lrts <- lrt(x_1s, x_ns, 4, 3.75, 20)
length(which(lrts < (3.75 / 4)^(4*20)*((4^2-o$minimum^2) / (3.75^2 - o$minimum^2))^(19))) / sims
#> [1] 0.04972

## power

n <- 20
sims <- 100000
x_1s <- numeric(sims)
x_ns <- numeric(sims)
reject <- logical(sims)
reject_joint <- logical(sims)
for (i in 1:sims) {
  p <- runif(n, 0, 1)
  x <- Finv(p, 3.75) # distribution under the alternate hypothesis
  x_1s[i] <- min(x)
  x_ns[i] <- max(x)
  reject[i] <- x_1s[i] < o$minimum # reject under the null hypothesis
  reject_joint[i] <- Fn_min(x_1s[i], 4, 20) - joint_F(x_1s[i], x_ns[i], 4, 20) < 0.05 # reject under the null hypothesis
}
# power when we use the independent min and max densities
sum(reject) / sims
#> [1] 0.05511
# power when we use the joint density
sum(reject_joint) / sims
#> [1] 0.05787

^{Created on 2023-12-04 with reprex v2.0.2}