2

I want to get a 95% confidence interval for the following question. enter image description here

I have written function f_n in my R code. I first randomly sample 100 with Normal and then I define function h for lambda. Then I can get f_n. My question is that how to define a function of f_n-chi-square and use uniroot` to find Confidence interval.

# I first get 100 samples 
set.seed(201111)
x=rlnorm(100,0,2)

Based on the answer by @RuiBarradas, I try the following code.

set.seed(2011111)
# I define function h, and use uniroot function to find lambda
h <- function(lam, n)
{
  sum((x - theta)/(1 + lam*(x - theta)))
}

# sample size
n <- 100
# the parameter of interest must be a value in [1, 12],
#true_theta<-1
#true_sd<- exp(2)
#x <- rnorm(n, mean = true_theta, sd = true_sd)
x=rlnorm(100,0,2)

xmax <- max(x)
xmin <- min(x)
theta_seq = seq(from = 1, to = 12, by = 0.01)

f_n <- rep(NA, length(theta_seq))

for (i in seq_along(theta_seq))
{
  theta <- theta_seq[i]
  lambdamin <- (1/n-1)/(xmax - theta)
  lambdamax <- (1/n-1)/(xmin - theta)
  lambda = uniroot(h, interval = c(lambdamin, lambdamax), n = n)$root
  f_n[i] = -sum(log(1 + lambda*(x - theta)))
}
j <- which.max(f_n)
max_fn <- f_n[j]
mle_theta <- theta_seq[j]

plot(theta_seq, f_n, type = "l", 
     main = expression(Estimated ~ theta),
     xlab = expression(Theta),
     ylab = expression(f[n]))
points(mle_theta, f_n[j], pch = 19, col = "red")
segments(
  x0 = c(mle_theta, xmin),
  y0 = c(min(f_n)*2, max_fn),
  x1 = c(mle_theta, mle_theta),
  y1 = c(max_fn, max_fn),
  col = "red",
  lty = "dashed"
)

I got the following plot of f_n.

enter image description here

For 95% CI, I try


LR <- function(theta, lambda)
{
  2*sum(log(1 + lambda*(x - theta))) - qchisq(0.95, df = 1)
}

lambdamin <- (1/n-1)/(xmax - mle_theta)
lambdamax <- (1/n-1)/(xmin - mle_theta)
lambda <- uniroot(h, interval = c(lambdamin, lambdamax), n = n)$root

uniroot(LR, c(xmin, mle_theta), lambda = lambda)$root

The result is 0.07198144. Then the logarithm is log(0.07198144)=-2.631347.

But there is NA in the following code.

uniroot(LR, c(mle_theta, xmax), lambda = lambda)$root

So the 95% CI is theta >= -2.631347.

But the question is that the 95% CI should be a closed interval...

Hermi
  • 357
  • 1
  • 11

1 Answers1

3

Here is a solution.

First of all, the data generation code is wrong, the parameter theta is in the interval [1, 12], and the data is generated with rnorm(., mean = 0, .). I change this to a true_theta = 5.

set.seed(2011111)
# I define function h, and use uniroot function to find lambda
h <- function(lam, n)
{
  sum((x - theta)/(1 + lam*(x - theta)))
}

# sample size
n <- 100
# the parameter of interest must be a value in [1, 12],
true_theta <- 5
true_sd <- 2
x <- rnorm(n, mean = true_theta, sd = true_sd)
xmax <- max(x)
xmin <- min(x)
theta_seq <- seq(from = xmin + .Machine$double.eps^0.5, 
                 to = xmax - .Machine$double.eps^0.5, by = 0.01)
f_n <- rep(NA, length(theta_seq))

for (i in seq_along(theta_seq))
{
  theta <- theta_seq[i]
  lambdamin <- (1/n-1)/(xmax - theta)
  lambdamax <- (1/n-1)/(xmin - theta)
  lambda = uniroot(h, interval = c(lambdamin, lambdamax), n = n)$root
  f_n[i] = -sum(log(1 + lambda*(x - theta)))
}
j <- which.max(f_n)
max_fn <- f_n[j]
mle_theta <- theta_seq[j]

plot(theta_seq, f_n, type = "l", 
     main = expression(Estimated ~ theta),
     xlab = expression(Theta),
     ylab = expression(f[n]))
points(mle_theta, f_n[j], pch = 19, col = "red")
segments(
  x0 = c(mle_theta, xmin),
  y0 = c(min(f_n)*2, max_fn),
  x1 = c(mle_theta, mle_theta),
  y1 = c(max_fn, max_fn),
  col = "red",
  lty = "dashed"
)


LR <- function(theta, lambda)
{
  2*sum(log(1 + lambda*(x - theta))) - qchisq(0.95, df = 1)
}

lambdamin <- (1/n-1)/(xmax - mle_theta)
lambdamax <- (1/n-1)/(xmin - mle_theta)
lambda <- uniroot(h, interval = c(lambdamin, lambdamax), n = n)$root

uniroot(LR, c(xmin, mle_theta), lambda = lambda)$root
#> [1] 4.774609

Created on 2022-03-25 by the reprex package (v2.0.1)

The one-sided CI95 is theta >= 4.774609.

Rui Barradas
  • 70,273
  • 8
  • 34
  • 66
  • @quasAliki The interval is not closed by construction. I have tried a two-sided interval and on the right, either the function `LR` doesn't change sign or it becomes `NaN` due to `log` of negative arguments. (That's why there was a bug in the 1st code I posted, `qchisq(0.05/2, 1)`is wrong, corrected to `qhisq(0.95, 1)`. See the edit.) – Rui Barradas Mar 25 '22 at 11:40
  • 1
    @quasAliki Why not consider the [normal distribution and its transformation to log-normal](https://en.wikipedia.org/wiki/Log-normal_distribution#Generation_and_parameters)? In this case, `exp(0)` gives `1` which is the minimum of your interval for `theta` and `exp(2)` gives the normal parameter `sigma`, the standard deviation. Apply the code above but generating `rnorm(100, 1, exp(2))` normal variates and compute the logarithm of the CI95 value. – Rui Barradas Mar 26 '22 at 12:40
  • Thanks! I apply your code and just modify `rnorm(100, 1, exp(2))`. But it still shows " LR doesn't change sign or it becomes NaN"? – Hermi Mar 27 '22 at 01:08
  • Also, why not use my original code for `theta_seq` as `theta_seq = seq(from = 1, to = 12, by = 0.01)`? – Hermi Mar 27 '22 at 01:14
  • @quasAliki Because the parameter cannot be smaller than `min(x)` nor greater than `max(x)`, but try your sequence and see how it goes. – Rui Barradas Mar 27 '22 at 06:37
  • For `x=rlnorm(100,0,2)`, is the expectation exp(2)=7.39? Why do you write `rnorm(100, 1, exp(2))`? – Hermi Mar 27 '22 at 20:30
  • So we can only find one root in `uniroot(LR, c(xmin, mle_theta), lambda = lambda)$root `? But there is no solution in `uniroot(LR, c(mle_theta, xmax), lambda = lambda)$root`? – Hermi Mar 27 '22 at 20:36
  • Can you please see my update? I use `uniroot(LR, c(xmin, mle_theta), lambda = lambda)$root` and got `0.07198144`. So I still need to take logarithm that is `log(0.07198144)=-2.631347`, right? So the 95% CI is `theta >= -2.631347`? – Hermi Mar 27 '22 at 20:50
  • @quasAliki Sorry, exp(2) is wrong, I was in a hurry when I wrote that. As for the CI95, the image you posted in the question asks for a one-sided interval like said above. And yes, that's the CI95. – Rui Barradas Mar 27 '22 at 21:36
  • That is ok. I think the 95% CI is a closed interval. Is there an error in the code? – Hermi Mar 27 '22 at 21:54
  • So if the expectation is `exp(2)`. how to transform `x=rlnorm(100,0,2)` into `rnorm`? Why not just used `x=rlnorm(100,0,2)`? – Hermi Mar 27 '22 at 21:57
  • @quasAliki I have already posted a link to the [Wikipedia](https://en.wikipedia.org/wiki/Log-normal_distribution) where it is explained how to transform from normal to log-normal and vice-versa, please read it. Why not just use `x=rlnorm(100,0,2)`? Because the normal distribution is very well behaved and the transformation to log-normal is a simple one. – Rui Barradas Mar 28 '22 at 11:35
  • I used `x=rlnorm(100,0,2)` in my question. Can you please see it? Do you think that is right? Because the 95% CI should be closed. I just would like to get a closed interval. – Hermi Mar 28 '22 at 13:05