FurtherMultiAnalysis.Rmd

---
title: "MultiFurtherAnalysis"
author: "Xiaohui Li Yuxin Ma"
date: "6/25/2017"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(mvtnorm)
library(reshape2)
source("PvalueFunction.R")
source("RatioPower.R")
```

# Check Y_max
```{r}
eps_level = 1
Y_i_exp <- NULL
Y_max <- NULL
for (i in 2:100){ # Treatment level
  K_i = i
  for (exp in 1:500){
    Y_exp = multi_y(eps_level, K_i)
    Y_i_exp[exp] = max(Y_exp)
  }
  Y_max[i] = mean(Y_i_exp)
}
plot(Y_max, type = "l", ylim = c(0,1), ylab = "Y_max", xlab = "K")
```

# Fitting
```{r}
K_max <- c(2:100)
lm_sqrt <- lm(Y_max[2:100] ~ 1/sqrt(K_max))
summary(lm_sqrt)
summary(lm_sqrt)[c("r.squared", "adj.r.squared")]
K_max_ln <- log(K_max, base = exp(1))
lm_max <- lm(Y_max[2:100] ~ K_max_ln) # epsilon = 1
summary(lm_max)
```

# Scaling (epsilon)
```{r}
N <- 10000
K <- 10
d <- 100
sigk <- rep(1, K)
Yk <- rep(0, K)
p_cv1 <- quantile(compute_pcombine(Yk, sigk, N, K, n = 10000 ) ,0.05)
p_cv2 <- quantile(compute_pcombine(Yk, 1.5*sigk, N, K, n = 10000 ) ,0.05)
p_cv3 <- quantile(compute_pcombine(Yk, 2*sigk, N, K, n = 10000 ) ,0.05)
eps1 <- 0
eps2 <- 0
eps3 <- 0
power_1 <- NULL
power_2 <- NULL
power_3 <- NULL
for (i in 1:d){ # Treatment level
  eps1 <- eps1 + 1/d
  eps2 <- eps2 + 1.5/d
  eps3 <- eps3 + 2/d
  power_exp1 <- NULL
  power_exp2 <- NULL
  power_exp3 <- NULL
  for (j in 1:100){
    Y_exp1 <- multi_y(eps1, K)
    Y_exp2 <- multi_y(eps2, K)
    Y_exp3 <- multi_y(eps3, K)
    power_exp1[j] <- compute_power(Y_exp1, sigk, N, K, p_cv1, R = 2000)
    power_exp2[j] <- compute_power(Y_exp2, 1.5*sigk, N, K, p_cv2, R = 2000)
    power_exp3[j] <- compute_power(Y_exp3, 2*sigk, N, K, p_cv3, R = 2000)
  }
  power_1[i] <- mean(power_exp1) 
  power_2[i] <- mean(power_exp2)
  power_3[i] <- mean(power_exp3)
}
plot(seq(0.01,1,0.01), power_1, type = "l", xlab = "Epsilon", ylab = "Power")
lines(seq(0.015,1.5,0.015), power_2, type = "l", col = "blue")
lines(seq(0.02,2,0.02), power_3, type = "l", col = "red")
legend("bottomright", c("Epsilon", "1.5Epsilon", "2Epsilon"), col = c("black","blue","red"), lty = 1)
```