Accept_reject_sampling.Rmd

---
title: "A3q4"
author: "Tianyi Fang"
date: "March 29, 2018"
output:
  word_document: default
  pdf_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

#Q4.
(1) How to get the best corn/fertilizer?
```{r}
dataq4 = read.csv("C:/Users/Tianyi Fang/Desktop/stat529/corn.csv", header = TRUE)
colnames(dataq4) = c('type', 'fer1', 'fer2', 'fer3')
#calculate thetan for each cell
theta_df = aggregate(dataq4[,2:4], list(dataq4$type), mean)
row.names(theta_df) = theta_df[,1]
theta_df = theta_df[,2:4]
#calculate lambdai.zetaj for each corn/fertilizer
theta_df['yi'] = rowMeans(theta_df)
theta_df = rbind(theta_df, unlist(colMeans(theta_df)))
row.names(theta_df) = c('corn1','corn2','corn3','corn4','wj')
```
## Case1: tao^2 >>sigma^2
generate 5000 data based on priors for lambda i for each corn, find the P(lambda(i)|data) the largest.
```{r}
get_largest = function(q4n, df, type){
  #initialize df
  df['beta'] = rnorm(q4n, mean = 125, sd = 15)
  df['tao'] = runif(q4n, min = 2, max = 16)
#lambda_df['theta'] = rnorm(q4n, mean = lambda_df$beta, sd = lambda_df$tao)
  if(type =='yi'){
    type_get = theta_df$yi[1:4]
  }
  if(type =='wj'){
    type_get = unlist(theta_df['wj',1:3])
  }
  
    for(i in 1:length(type_get)){
      #print(head(theta_df[type]))
      df[,i+3] = rnorm(q4n, mean = type_get[i], sd = sqrt(df$tao/3))
}
    # 0/1 coding for largest lambda
    df_new = t(apply(df[,-(1:3)], 1, function(x) {1*(x==max(x))}))
    #get count(theta = 1 in the last row) and precentage
    #remember to drop the initial theta
    new_count = colSums(df_new)
    new_prec = new_count/q4n
    return(list(df, df_new, new_count, new_prec))
}

#get result of lambda
q4n = 5000
lambda_df = data.frame('index' = seq(1:q4n))
lambda_result = get_largest(q4n, lambda_df, 'yi')
print(paste('P(lambda2 = lambda[4]|data) is ', lambda_result[[4]][2]))
#get result of zeta
zeta_df = data.frame('index' = seq(1:q4n))
zeta_result = get_largest(q4n, zeta_df, 'wj')
print(paste('P(W1 = W[3]|data) is ', zeta_result[[4]][1]))
```
P(theta12>=theta4|data):
```{r}
x_bar12 = theta_df$fer1[4]
x_bar4 = theta_df$fer1[2]
sigma124 = 40
w = rnorm(50000, mean = x_bar12-x_bar4, sd = sqrt(sigma124))
print(paste('P(theta12>=theta4|data) in case 1 is', sum(w<=0)/length(w)))
```
we can see that theta2 is not always the largest. Because there are interaction between the corn and fertilizer. Using HB to calculate
## Case2. HB
```{r}
q4n = 5000
#get the result of HB method
get_hb = function(q4n, df, type){
    if(type =='yi'){
        type_get = theta_df$yi[1:4]
   }
    if(type =='wj'){
        type_get = unlist(theta_df['wj',1:3])
   }
   n = length(type_get)
    for(i in 1:n){
  #generate lambda1..lambda4
       df[i+3] = rnorm(5000, mean = df$beta, sd = sqrt(df$tao))
    }
    for(i in 1:n){
  #generate squre_sum
       df[i+n+3] = (df[i+3]- type_get[i])^2
    }
    for(i in 1:n){
  #generate ratio
       df[i+2*n+3] = exp(-0.5*12/80*df[i+n+3])
    }

     keeper = list()
  #find keeper
     for(i in 1:n){
         keeper[[i]] = df[i+3][df[i+2*n+3]>=runif(5000,0,1)] 
    }
#calculate the mean of keeper
     print(lapply(keeper, function(x) {mean(x)}))
     return(list(df, keeper))
}
#compare the result 
compare_result = function(keeper, type, w){
  m = min(unlist(lapply(keeper, function(x){length(x)})))
  print(m)
  compare_df = data.frame(matrix(NA, nrow = m, ncol = length(keeper)))
  for(i in 1:length(keeper)){
    compare_df[,i] = keeper[[i]][1:m]
  }
  compare_df['max'] = unlist(apply(compare_df,1, which.max))
  print(paste('P(',type,w, '=', type,'[4]|data) is:', sum(compare_df['max'] == w)/nrow(compare_df)))
  return(list(m,compare_df))
}
```
since after sampling, the min number of lambda is 741, truncate all other into 741, then get the precentage of lambda2 = lmabda[4]
```{r}
#get the HB result for corn
hb_ldf = data.frame('index' = seq(1,q4n))
hb_ldf['beta'] = rnorm(q4n, mean = 125, sd = 15)
hb_ldf['tao'] = runif(q4n,2,16)
corn_keeper = get_hb(q4n, hb_ldf, 'yi')[[2]]
corn_result = compare_result(corn_keeper, 'lambda',2)

#get the HB result for fertilizer
hb_fdf = data.frame('index' = seq(1,q4n))
hb_fdf['beta'] = rnorm(q4n, mean = 125, sd = 15)
hb_fdf['tao'] = runif(q4n,2,16)
fer_keeper = get_hb(q4n, hb_fdf, 'wj')[[2]]
fer_result = compare_result(fer_keeper, 'W', 1)
```
P(theta12>=theta4|data):
theta12 is (corn4, fer1), theta4 is (corn2, fer1), since we have samplers from above results, we can directly compare them.
x_bar_4 = 135, x_bar_12 = 131
```{r}
x_bar_4 = 135
x_bar_12 = 131
cc = data.frame('index' = seq(1,5000))
cc['beta'] = rnorm(5000, 125, 15)
cc['tao'] = runif(5000, 2,16)

get_hb = function(q4n, df){
   n = 4
   df['theta4'] = rnorm(5000, mean = df$beta, sd = sqrt(df$tao))
   df['sse4'] = df['theta4']- x_bar_4
   df['ratio4'] = exp(-0.5*df['sse4']/20)
   df['theta12'] = rnorm(5000, mean= df$beta, sd = sqrt(df$tao))
   df['sse12'] = df['theta12']- x_bar_12
   df['ratio12'] = exp(-0.5*df['sse12']/20)

    keeper = list()
  #find keeper
    keeper[[1]] = df['theta4'][df['ratio4']>=runif(5000,0,1)] 
    keeper[[2]] = df['theta12'][df['ratio12']>=runif(5000,0,1)] 
#calculate the mean of keeper
     print(lapply(keeper, function(x) {mean(x)}))
     return(list(df, keeper))
}
compare412 = get_hb(q4n, cc)[[2]]
c = min(unlist(lapply(compare412, function(x){length(x)})))
compare412_hb = data.frame(matrix(NA, nrow = c, ncol = length(compare412)))
for(i in 1:length(compare412)){
    compare412_hb[,i] = compare412[[i]][1:c]
  }
compare412_hb['max'] = unlist(apply(compare412_hb,1, which.max))
print(paste('P(theta12>theta4|data) in case 1 is:', sum(compare412_hb['max'] == 2)/nrow(compare412_hb)))
```
#Extra
tao~U(2,3)
```{r}
x_bar_4 = 135
x_bar_12 = 131
cc2 = data.frame('index' = seq(1,5000))
cc2['beta'] = rnorm(5000, 125, 15)
cc2['tao'] = runif(5000, 2,3)

get_hb = function(q4n, df){
   n = 4
   df['theta4'] = rnorm(5000, mean = df$beta, sd = sqrt(df$tao))
   df['sse4'] = df['theta4']- x_bar_4
   df['ratio4'] = exp(-0.5*df['sse4']/20)
   df['theta12'] = rnorm(5000, mean= df$beta, sd = sqrt(df$tao))
   df['sse12'] = df['theta12']- x_bar_12
   df['ratio12'] = exp(-0.5*df['sse12']/20)

    keeper = list()
  #find keeper
    keeper[[1]] = df['theta4'][df['ratio4']>=runif(5000,0,1)] 
    keeper[[2]] = df['theta12'][df['ratio12']>=runif(5000,0,1)] 
#calculate the mean of keeper
     print(lapply(keeper, function(x) {mean(x)}))
     return(list(df, keeper))
}
compare4122 = get_hb(q4n, cc2)[[2]]
c2 = min(unlist(lapply(compare4122, function(x){length(x)})))
compare4122_hb = data.frame(matrix(NA, nrow = c2, ncol = length(compare4122)))
for(i in 1:2){
    compare4122_hb[,i] = compare4122[[i]][1:c2]
  }
compare4122_hb['max'] = unlist(apply(compare4122_hb,1, which.max))
print(paste('P(theta12>theta4|data) in case 2 is:', sum(compare4122_hb['max'] == 2)/nrow(compare4122_hb)))

```