m1 <- 12 
m2 <- 18 # change to evaluate the contribution of effect size

mean_11 = m1
mean_21 = m2

var1 <- 40
iterations <- 1000
types1 <- 5:40

p_val <- array(0,iterations)
percent1 <- array(0,length(types1))

for (j in 1:length(types1) ){
var1 = 40

samp_size <- types1[j]

for (i in 1:iterations){
x1 <- sample(rnorm(samp_size,mean_11,sqrt(var1)))
x2 <- sample(rnorm(samp_size,mean_21,sqrt(var1)))
rslt <- t.test(x1, x2, alternative = "two.sided", var.equal = TRUE)
p_val[i] <- rslt$p.value
}

#hist(p_val,breaks = c(seq(-0.01,0.1,0.01),seq(0.2,1,0.1) ) )

percent1[j] <- sum( (p_val < 0.05)*1 )/iterations
}
    

###################################################################

mean_12 = m1
mean_22 = m2

samp_size <- 25
iterations <- 1000
types2 <- 5:60

var <- rep(0,length(types2))
p_val <- array(0,iterations)
percent2 <- array(0,length(types2))

for (j in 1:length(types2) ){
  
  var[j] <-types2[j]
  for (i in 1:iterations){
    x1 <- sample(rnorm(samp_size,mean_12,sqrt(var[j])))
    x2 <- sample(rnorm(samp_size,mean_22,sqrt(var[j])))
    rslt <- t.test(x1, x2, alternative = "two.sided", var.equal = TRUE)
    p_val[i] <- rslt$p.value
  }
  
  #hist(p_val,breaks = c(seq(-0.01,0.1,0.01),seq(0.2,1,0.1) ) )
  
  percent2[j] <- sum( (p_val <0.05)*1 )/iterations
}

###################################################################

mean_12 = 12
mean_22 = 18

samp_size <- 40
iterations <- 1000
types2 <- 5:60

var <- rep(0,length(types2))
p_val <- array(0,iterations)
percent3 <- array(0,length(types2))

for (j in 1:length(types2) ){
  
  var[j] <-types2[j]
  for (i in 1:iterations){
    x1 <- sample(rnorm(samp_size,mean_12,sqrt(var[j])))
    x2 <- sample(rnorm(samp_size,mean_22,sqrt(var[j])))
    rslt <- t.test(x1, x2, alternative = "two.sided", var.equal = TRUE)
    p_val[i] <- rslt$p.value
  }
  
  #hist(p_val,breaks = c(seq(-0.01,0.1,0.01),seq(0.2,1,0.1) ) )
  
  percent3[j] <- sum( (p_val <0.05)*1 )/iterations
}




par(mfrow=c(2,3))
hist(rnorm(10000,mean_11,sqrt(var1)),col=rgb(0,0,1,1/4),xlim=c(0,30),breaks = 100,
main = paste("means A: ",m1, "B =", m2, " var= ",var1), ylab ="Frequency",
xlab = "x")
hist(rnorm(10000,mean_21,sqrt(var1)),add=T,col=rgb(1,0,0,1/4),xlim=c(0,30),breaks = 100)
hist(rnorm(10000,mean_12,sqrt(var[1])),col=rgb(0,0,1,1/4),xlim=c(0,30),breaks = 100,
main = paste("means A: ",m1, "B =", m2, " var= ",var[1]), ylab ="Frequency",
xlab = "x")
hist(rnorm(10000,mean_22,sqrt(var[1])),add=T,col=rgb(1,0,0,1/4),xlim=c(0,30),breaks = 100)
hist(rnorm(10000,mean_12,sqrt(var[16])),col=rgb(0,0,1,1/4),xlim=c(0,30),breaks = 100,
main = paste("means A: ",m1, "B =", m2, " var= ",var[16]), ylab ="Frequency",
xlab = "x")
hist(rnorm(10000,mean_22,sqrt(var[16])),add=T,col=rgb(1,0,0,1/4),xlim=c(0,30),breaks = 100)

plot(types1,percent1,ylim=c(0,1),pch=16, ylab= "Frequency P < 0.05", xlab= "Sample size",
     main = paste("variance =", var1))
points(25,percent1[which(types1==25)],cex=3)
points(40,percent1[which(types1==40)],cex=3) 
plot(types2,percent2,ylim=c(0,1),pch=16, ylab= "Frequency P < 0.05", xlab= "Variance",
     main = paste("sample size =", 25)) 
points(40,percent2[which(types2==40)],cex=3)
plot(types2,percent3,ylim=c(0,1),pch=16, ylab= "Frequency P < 0.05", xlab= "Variance",
     main = paste("sample size =", 40)) 
points(40,percent3[which(types2==40)],cex=3)