Q5:

Part b) Box plot and P_values

Let \(\theta_2\) = 0.1 and \(\theta_1\) = 0. Assume n = 10,100,1000, and 10,000. For each n value, generate a data set of D = {\(x_1,..., x_n, y_1,...,y_n\)} by random number generator in your favorite computer program. Make two box plots for {\(x_1,...,x_n\)} and {\(y_1,...,y_n\)}, respectively, for each n value. Calculate the p-value for each n.

par(mfrow = c(2,4))
set.seed(456)
n = c(10,100,1000,10000)
X1 = rnorm(n[1], mean = 0, sd = 1) 
Y1 = rnorm(n[1], mean = 0.1, sd = 1) 
boxplot(X1, xlab = "Distribution of X", main = paste("Sample Size: n = ", n[1]))
boxplot(Y1, xlab = "Distribution of Y", main = paste("Sample Size: n = ", n[1]))
p1 = t.test(X1,Y1, var.equal = T)[3]

X2 = rnorm(n[2], mean = 0, sd = 1) 
Y2 = rnorm(n[2], mean = 0.1, sd = 1) 
boxplot(X2, xlab = "Distribution of X", main = paste("Sample Size: n = ", n[2]))
boxplot(Y2, xlab = "Distribution of Y", main = paste("Sample Size: n = ", n[2]))
p2 = t.test(X2,Y2, var.equal = T)[3]


X3 = rnorm(n[3], mean = 0, sd = 1) 
Y3 = rnorm(n[3], mean = 0.1, sd = 1) 
boxplot(X3, xlab = "Distribution of X", main = paste("Sample Size: n = ", n[3]))
boxplot(Y3, xlab = "Distribution of Y", main = paste("Sample Size: n = ", n[3]))
p3 = t.test(X3,Y3, var.equal = T)[3]

X4 = rnorm(n[4], mean = 0, sd = 1) 
Y4 = rnorm(n[4], mean = 0.1, sd = 1) 
boxplot(X3, xlab = "Distribution of X", main = paste("Sample Size: n = ", n[4]))
boxplot(Y3, xlab = "Distribution of Y", main = paste("Sample Size: n = ", n[4]))

p4 =t.test(X4,Y4, var.equal = T)[3]


p_values = cbind(p1,p2,p3,p4)
p_values
##         p1         p2        p3        p4          
## p.value 0.04735338 0.6400347 0.2943447 1.461134e-12

Part c) Bayes Factor

get_BF = function(X,Y,n){
  D = append(X,Y)
  theta_mle = mean(D);theta_1 = mean(X); theta_2 = mean(Y)
  BIC_H0 = -2 * sum(log(dnorm(D, theta_mle))) + log(n)
  BIC_H1 = -2 * sum(log(dnorm(X, theta_1))+ log(dnorm(Y, theta_2))) + 2 * log(n)
  delta_BIC = BIC_H1 - BIC_H0
  BF = exp(delta_BIC/2)
  BF
}

BF_D1 = get_BF(X1,Y1,10)
BF_D2 = get_BF(X2,Y2,100)
BF_D3 = get_BF(X3,Y3,1000)
BF_D4 = get_BF(X4,Y4,10000)

(BF = cbind(BF_D1,BF_D2,BF_D3,BF_D4))
##          BF_D1    BF_D2    BF_D3        BF_D4
## [1,] 0.2214016 9.076242 18.55832 8.137164e-10
posterior_prob_D1 = BF_D1 / (1+BF_D1)
posterior_prob_D2 = BF_D2 / (1+BF_D2)
posterior_prob_D3 = BF_D3 / (1+BF_D3)
posterior_prob_D4 = BF_D4 / (1+BF_D4)


round((posterior_prob = cbind(posterior_prob_D1,posterior_prob_D2,posterior_prob_D3,posterior_prob_D4)),3)
##      posterior_prob_D1 posterior_prob_D2 posterior_prob_D3 posterior_prob_D4
## [1,]             0.181             0.901             0.949                 0

Comment:

The P(\(H_{o}|D\)) for the four data sets are: 0.612, 0.82, 0.609 and 0.0 respectively.

Part d) Repeat part b) and c) with different \(\theta\)

set.seed(456)
generate_data = function(n,theta){
  x = rnorm(n, mean = 0, sd = 1) 
  y = rnorm(n, mean = theta, sd = 1)
  data = list("x" = x,"y" = y)
  data
}
get_p_val = function(X,Y){
  p = t.test(X,Y, var.equal = T)[3]
  p
}
n = c(10,100,1000,10000)
# n = 10
df_theta_0.1 =generate_data(n[1],0.1)
df_theta_0.05 =generate_data(n[1],0.05)
df_theta_0.01 =generate_data(n[1],0.01)
df_theta_0.001 =generate_data(n[1],0.001)
df_theta = list(df_theta_0.1, df_theta_0.05,df_theta_0.01,df_theta_0.001)


# p_val from t test
p_val_10_obs = rep(0,4)
for(i in 1:4){
  p_val_10_obs[i] = get_p_val(df_theta[[i]]$x,df_theta[[i]]$y)
}

p_val_10_obs = unlist(p_val_10_obs)

p_0.1 = get_p_val(df_theta_0.1$x,df_theta_0.1$y)
p_0.05 = get_p_val(df_theta_0.1$x,df_theta_0.1$y)
p_0.01 = get_p_val(df_theta_0.1$x,df_theta_0.1$y)
p_0.001 = get_p_val(df_theta_0.1$x,df_theta_0.1$y)

# posterior prob
theta = c(0.1,0.05,0.01,0.001)
posterior_10 = rep(0,4)
BF_theta_0.1 = get_BF(df_theta_0.1$x,df_theta_0.1$y,n[1]) 
posterior_10[1] = BF_theta_0.1/ (1+BF_theta_0.1)
  
BF_theta_0.05 = get_BF(df_theta_0.05$x,df_theta_0.05$y,n[1]) 
posterior_10[2] = BF_theta_0.05/ (1+BF_theta_0.05)
  
BF_theta_0.01 = get_BF(df_theta_0.01$x,df_theta_0.01$y,n[1])
posterior_10[3] = BF_theta_0.01/ (1+BF_theta_0.01)
  
BF_theta_0.001 = get_BF(df_theta_0.001$x,df_theta_0.001$y,n[1]) 
posterior_10[4] = BF_theta_0.001/ (1+BF_theta_0.001)
 
posterior_10
## [1] 0.1812685 0.7541259 0.7438201 0.7596865
# n = 100
df_theta_0.1 =generate_data(n[2],0.1)  
df_theta_0.05 =generate_data(n[2],0.05)
df_theta_0.01 =generate_data(n[2],0.01)
df_theta_0.001 =generate_data(n[2],0.001)

df_theta = list(df_theta_0.1, df_theta_0.05,df_theta_0.01,df_theta_0.001)

# p_val from t_test
p_val_100_obs = rep(0,4)
for(i in 1:4){
  p_val_100_obs[i] = get_p_val(df_theta[[i]]$x,df_theta[[i]]$y)
}

p_val_100_obs = unlist(p_val_100_obs)

# posterior prob
posterior_100 = rep(0,4)
BF_theta_0.1 = get_BF(df_theta_0.1$x,df_theta_0.1$y,n[2]) 
posterior_100[1] = BF_theta_0.1/ (1+BF_theta_0.1)
  
BF_theta_0.05 = get_BF(df_theta_0.05$x,df_theta_0.05$y,n[2]) 
posterior_100[2] = BF_theta_0.05/ (1+BF_theta_0.05)
  
BF_theta_0.01 = get_BF(df_theta_0.01$x,df_theta_0.01$y,n[2])
posterior_100[3] = BF_theta_0.01/ (1+BF_theta_0.01)
  
BF_theta_0.001 = get_BF(df_theta_0.001$x,df_theta_0.001$y,n[2]) 
posterior_100[4] = BF_theta_0.001/ (1+BF_theta_0.001)
# n = 1000
df_theta_0.1 =generate_data(n[3],0.1)  
df_theta_0.05 =generate_data(n[3],0.05)
df_theta_0.01 =generate_data(n[3],0.01)
df_theta_0.001 =generate_data(n[3],0.001)

# p_val from t test
df_theta = list(df_theta_0.1, df_theta_0.05,df_theta_0.01,df_theta_0.001)
p_val_1000_obs = rep(0,4)
for(i in 1:4){
  p_val_1000_obs[i] = get_p_val(df_theta[[i]]$x,df_theta[[i]]$y)
}

p_val_1000_obs = unlist(p_val_1000_obs)

#posterior prob
posterior_1000 = rep(0,4)
BF_theta_0.1 = get_BF(df_theta_0.1$x,df_theta_0.1$y,n[3]) 
posterior_1000[1] = BF_theta_0.1/ (1+BF_theta_0.1)
  
BF_theta_0.05 = get_BF(df_theta_0.05$x,df_theta_0.05$y,n[3]) 
posterior_1000[2] = BF_theta_0.05/ (1+BF_theta_0.05)
  
BF_theta_0.01 = get_BF(df_theta_0.01$x,df_theta_0.01$y,n[3])
posterior_1000[3] = BF_theta_0.01/ (1+BF_theta_0.01)
  
BF_theta_0.001 = get_BF(df_theta_0.001$x,df_theta_0.001$y,n[3]) 
posterior_1000[4] = BF_theta_0.001/ (1+BF_theta_0.001)
# n = 10000
df_theta_0.1 =generate_data(n[4],0.1)  
df_theta_0.05 =generate_data(n[4],0.05)
df_theta_0.01 =generate_data(n[4],0.01)
df_theta_0.001 =generate_data(n[4],0.001)
df_theta = list(df_theta_0.1, df_theta_0.05,df_theta_0.01,df_theta_0.001)

# p_val from t-test
p_val_10000_obs = rep(0,4)
for(i in 1:4){
  p_val_10000_obs[i] = get_p_val(df_theta[[i]]$x,df_theta[[i]]$y)
}

p_val_10000_obs = unlist(p_val_10000_obs)

# posterior prob
posterior_10000 = rep(0,4)
BF_theta_0.1 = get_BF(df_theta_0.1$x,df_theta_0.1$y,n[4]) 
posterior_10000[1] = BF_theta_0.1/ (1+BF_theta_0.1)
  
BF_theta_0.05 = get_BF(df_theta_0.05$x,df_theta_0.05$y,n[4]) 
posterior_10000[2] = BF_theta_0.05/ (1+BF_theta_0.05)
  
BF_theta_0.01 = get_BF(df_theta_0.01$x,df_theta_0.01$y,n[4])
posterior_10000[3] = BF_theta_0.01/ (1+BF_theta_0.01)
  
BF_theta_0.001 = get_BF(df_theta_0.001$x,df_theta_0.001$y,n[4]) 
posterior_10000[4] = BF_theta_0.001/ (1+BF_theta_0.001)
 
posterior = rbind(posterior_10,posterior_100,posterior_1000,posterior_10000)
colnames(posterior) = c("theta = 0.1","theta = 0.05","theta = 0.01","theta = 0.001" )

# posterior prob
round(posterior,3)
##                 theta = 0.1 theta = 0.05 theta = 0.01 theta = 0.001
## posterior_10          0.181        0.754        0.744         0.760
## posterior_100         0.472        0.896        0.854         0.905
## posterior_1000        0.969        0.530        0.887         0.960
## posterior_10000       0.000        0.597        0.984         0.989
# p_values from t_test
p_val = rbind(p_val_10_obs,p_val_100_obs,p_val_1000_obs,p_val_10000_obs)
colnames(p_val) = c("theta = 0.1","theta = 0.05","theta = 0.01","theta = 0.001" )
round(p_val,3)
##                 theta = 0.1 theta = 0.05 theta = 0.01 theta = 0.001
## p_val_10_obs          0.047        0.806        0.589         0.981
## p_val_100_obs         0.023        0.565        0.303         0.776
## p_val_1000_obs        0.970        0.009        0.100         0.450
## p_val_10000_obs       0.000        0.004        0.327         0.698

Comment:

From comparing columns (different values of \(\theta\)), I noticed that when having relatively larger theta(left two columns vs. right two columns in both tables), in general, we are more likely to reject the null and conclude that X and Y have significant difference. When \(\theta\) is small, both perspectives fail to reject the null since X and Y barely have difference given really small \(\theta\).

From comparing rows, we can see that when we have relatively small sample size (first two row of both tables, n = 10, 100), the p_values and the posterior probability are roughly the same scale. As we move on to the larger sample size, the situation is different. Since increasing the sample size makes the hypothesis test more sensitive, thus more powerful. I also noticed that we are more likely to reject the null hypothesis if using p_value. From comparison, the p_value table has a smaller scale than the posterior probability. As we can see, for example, given sample size 10000 (\(\theta\) = 0.05), the p_value is 0.004. In this case, if using p_value, we would reject the null since 0.004 < 5% significance threshold. However, the \(H_{o}|D\) is 0.597, we fail to reject the null. This example shows us, using p_value could be misleading, while using posterior probability makes more sense. In Bayesian statistics, we transform them to direct measures the evidence against the null hypothesis(Bayes factor).