js_estimator.R

# JS estimator for Normal scale

Atrue <- .2
n <- 1000
res <- t(sapply(1:100, function(i) {
  set.seed(i)
  theta <- rnorm(n, 0, sqrt(Atrue))
  D <- rexp(n,rate=.1)
  X <- rnorm(n, theta, sqrt(D))
  (rough <- var(X) - mean(D))
  S <- X^2
  I <- function(A) 1/(2*(A + D)^2)
  Ahat <- function(A) sum((S - D) * I(A)) / sum(I(A))
  objective <- function(A) Ahat(A) - A
  minA <- .01
  if (objective(minA) < 0) {
    zero <- minA
  } else {
    zero <- uniroot(objective, interval=c(.01, 20))$root
  }
  #s <- seq(from=.01, to=2, by=.01)
  #plot(s,sapply(s, objective))
  c(rough, zero)
}))
head(res)
boxplot(res)
abline(h=Atrue, col="blue")