StatComp_2021Q1Q2Q3

## 2021-Q1 ##############################
sim_n <- 10000; n1 <- 10; n2 <- 100

X_n1 <- matrix(rpois(sim_n*n1, lambda = 5), nrow = sim_n, ncol = n1)

theta <- exp(-5)
t1_n1 <- apply(X_n1, 1, function(x){length(x[x==0])/n1})
t2_n1 <- apply(X_n1, 1, function(x){sum((-1)^x)/n1})
t3_n1 <- apply(X_n1, 1, function(x){exp(-sum(x)/n1)})
t4_n1 <- apply(X_n1, 1, function(x){(1-1/n1)^sum(x)})

t1_bar_n1 <- mean(t1_n1)
t2_bar_n1 <- mean(t2_n1)
t3_bar_n1 <- mean(t3_n1)
t4_bar_n1 <- mean(t4_n1)

bias1_n1 <- t1_bar_n1 - theta
bias2_n1 <- t2_bar_n1 - theta
bias3_n1 <- t3_bar_n1 - theta
bias4_n1 <- t4_bar_n1 - theta

SE1_n1<-sqrt(sum((t1_n1 - t1_bar_n1)^2)/sim_n)
SE2_n1<-sqrt(sum((t2_n1 - t2_bar_n1)^2)/sim_n)
SE3_n1<-sqrt(sum((t3_n1 - t3_bar_n1)^2)/sim_n)
SE4_n1<-sqrt(sum((t4_n1 - t4_bar_n1)^2)/sim_n)

R_mse1_n1<-sqrt(sum((t1_n1 - theta)^2)/sim_n)
R_mse2_n1<-sqrt(sum((t2_n1 - theta)^2)/sim_n)
R_mse3_n1<-sqrt(sum((t3_n1 - theta)^2)/sim_n)
R_mse4_n1<-sqrt(sum((t4_n1 - theta)^2)/sim_n)


X_n2 <- matrix(rpois(sim_n*n2, lambda = 5), nrow = sim_n, ncol = n2)

theta <- exp(-5)
t1_n2 <- apply(X_n2, 1, function(x){length(x[x==0])/n2})
t2_n2 <- apply(X_n2, 1, function(x){sum((-1)^x)/n2})
t3_n2 <- apply(X_n2, 1, function(x){exp(-sum(x)/n2)})
t4_n2 <- apply(X_n2, 1, function(x){(1-1/n2)^sum(x)})

t1_bar_n2 <- mean(t1_n2)
t2_bar_n2 <- mean(t2_n2)
t3_bar_n2 <- mean(t3_n2)
t4_bar_n2 <- mean(t4_n2)

bias1_n2 <- t1_bar_n2 - theta
bias2_n2 <- t2_bar_n2 - theta
bias3_n2 <- t3_bar_n2 - theta
bias4_n2 <- t4_bar_n2 - theta

SE1_n2<-sqrt(sum((t1_n2 - t1_bar_n2)^2)/sim_n)
SE2_n2<-sqrt(sum((t2_n2 - t2_bar_n2)^2)/sim_n)
SE3_n2<-sqrt(sum((t3_n2 - t3_bar_n2)^2)/sim_n)
SE4_n2<-sqrt(sum((t4_n2 - t4_bar_n2)^2)/sim_n)

R_mse1_n2<-sqrt(sum((t1_n2 - theta)^2)/sim_n)
R_mse2_n2<-sqrt(sum((t2_n2 - theta)^2)/sim_n)
R_mse3_n2<-sqrt(sum((t3_n2 - theta)^2)/sim_n)
R_mse4_n2<-sqrt(sum((t4_n2 - theta)^2)/sim_n)

result<-{data.frame("n" = c(n1," "," "," "," ",n2," "," "," "," "),
                    "Mean estimate" = c(" ",t1_bar_n1,t2_bar_n1,t3_bar_n1,t4_bar_n1," ",t1_bar_n2,t2_bar_n2,t3_bar_n2,t4_bar_n2),
                    "Bias" = c(" ",bias1_n1,bias2_n1,bias3_n1,bias4_n1," ",bias1_n2,bias2_n2,bias3_n2,bias4_n2),
                    "Sample SE" = c(" ",SE1_n1,SE2_n1,SE3_n1,SE4_n1," ",SE1_n2,SE2_n2,SE3_n2,SE4_n2),
                    "Root MSE" = c(" ",R_mse1_n1,R_mse2_n1,R_mse3_n1,R_mse4_n1," ",R_mse1_n2,R_mse2_n2,R_mse3_n2,R_mse4_n2),
                    row.names = c(" ","theta1_hat","theta2_hat","theta3_hat","theta4_hat","  ","theta1-hat","theta2-hat","theta3-hat","theta4-hat"))
}
result
# write.csv(result,"C:/Users/user/Desktop/result.csv")

## 2 ########################

sim_n2 <- 5000

set.seed(123)
X <- rep(0, sim_n2)
for(i in 1:sim_n2)
{
  v <- vector()
  j <- 1
  k <- 0
  repeat{
    x <- sample(c(1:6), 1)
    k <- k + 1
    if(length(v[v == x]) == 0)
    {
      v[j] <- x
      j <- j+1
    }
    if(length(v) == 6)
      break
  }
  X[i] <- k
}
exact2 <- 6*sum(1/c(1:6))
exact2 # 14.7
mean(X) # 14.6988

## 3 ###########################

sim_n3 <- 5000
n <- 30
set.seed(123)
X <- matrix(rnorm(replicate*n), ncol=n, nrow=sim_n3)
Y <- X^2
ans <- matrix(0, nrow =7, ncol = 2)

# (1) #####
temp <- 20*X[,19] + 21*X[,20]
temp <- temp[temp<12/16]
ans[1,1] <- pnorm(12/16, mean = 0, sd = sqrt(20^2+21^2))
ans[1,2] <- length(temp)/sim_n3

# (2) #####
temp <- apply(Y[,1:12], 1,function(x){sqrt(sum(x))})
ans[2,1] <- sqrt(qchisq(.8, df = 12))
ans[2,2] <- quantile(temp, prob = .8)

# (3) #####
temp1 <- apply(Y[,1:12], 1, sum)
temp2 <- apply(Y[,1:16], 1, sum)
temp <- temp1/temp2
ans[3,1] <- NA
ans[3,2] <- median(temp)

# (4) #####
temp <- apply(X, 1, function(x) quantile(x, prob= .5+.2021))
ans[4,1] <- qnorm(.5+.2021)
ans[4,2] <- mean(temp)

# (5) #####
temp <- (Y[,1])/(Y[,2])
ans[5,1] <- qf(.2021,1,1)
ans[5,2] <- quantile(temp, prob = .2021)

temp <- (X[,1])/(X[,2])
ans[6,1] <- qcauchy(.2021)
ans[6,2] <- quantile(temp, prob = .2021)

# (6) #####
temp1 <- X[,1]
temp2 <- apply(Y[,7:15],1,sum)
temp2 <- sqrt(temp2)
temp <- temp1[temp1<temp2]
ans[7,1] <- pt(3,df=9)
ans[7,2] <- length(temp)/sim_n3

ans <- data.frame(ans)
# write.csv(ans,"C:/Users/user/Desktop/ans.csv")