plot sample mean and variance

asymp_meanvar.R

moments_4 <- function(truedens){
  k1 <- integrate(function(x){x*truedens(x)}, lower=-Inf, upper=Inf)
  k2 <- integrate(function(x){(x-k1$value)^2*truedens(x)}, lower=-Inf, upper=Inf)
  k3 <- integrate(function(x){(x-k1$value)^3*truedens(x)}, lower=-Inf, upper=Inf)
  k4 <- integrate(function(x){(x-k1$value)^4*truedens(x)}, lower=-Inf, upper=Inf)
  c(k1$value, k2$value, 
    k3$value, k4$value)  
}

Mu4 <- moments_4(truedens = function(x){dgamma(x,2,1)})
Cov <- matrix( c(Mu4[2],Mu4[3],Mu4[3],Mu4[4]), 2 ,2)

sim_meanvar <- function(rand, n_s, mu, sigma2, iter=10000){
  mhat <- numeric(iter)
  shat <- numeric(iter)
  for(i in seq_len(iter)){
    x <- rand(n_s)
    mhat[i] <- sqrt(n_s)*(mean(x)-mu)
    shat[i] <- sqrt(n_s)*(var(x)-sigma2)
  }
  data.frame(mean=mhat,
             variance=shat,
             n=n_s)  
}

n_s <- c(10,100,1000)

rand <- function(n){rgamma(n,2,1)}
ressim <- lapply(n_s, function(ns){sim_meanvar(rand = rand, ns, mu = 2, sigma2 = 2)})
ressim <- bind_rows(ressim)

ran_x <- range(ressim$mean)
ran_y <- range(ressim$variance)
df_dens <- expand.grid(
  m=seq(ran_x[1], ran_x[2], length.out=201),
  s=seq(ran_y[1], ran_y[2], length.out=201)) %>% 
  rowwise() %>% 
  mutate(dens = dmvnorm(c(m,s), sigma = Cov))

p1 <- ggplot()+
  geom_point(data = ressim, aes(x=mean, y=variance), alpha=0.05) +
  geom_contour(data=df_dens, aes(x=m, y=s, z=dens, colour=after_stat(level))) +
  scale_color_viridis_c()+
  facet_grid(cols=vars(n), labeller = label_both)+
  theme_bw(15)+labs(x="sample mean", y="sample variance")

print(p1)
ggsave("meanvar.png", width = 10, height = 7)