StatComp_2019Q1Q2Q3

#2019-Q1------------------------------------------------------------------------------------
z <- c(2:5)
exact<-10+sum(z*(((choose(5,z))*(0.5^5))/(1-6*(0.5^5))))+10

f2  <-function(x){                                 #pdf of X2
  ((choose(5,x))*(0.5^5))/(1-6*(0.5^5))
}
replicate <- 10000
X1 <- rep(0, replicate)
X2 <- rep(0, replicate)
for(i in 1:replicate){
  x1 <- 0
  x2 <- 0
  while(TRUE){
    u <- runif(1)
    x2 <- x2+1
    if(u<1/3){
      x1 <- x1+rgamma(1,2,scale = 5)
      break
    }
    else{
      if(u<2/3){
        u2 <- runif(1)
        x <- 2
        p <- f2(x)
        while(u2>p){
          x <- x+1
          p <- p+f2(x)
        }
        x1 <- x1+x
      }
      else{
        x1 <- x1 +rexp(1,rate = 1/10)
      }
    }
  }
  X1[i] <- x1
  X2[i] <- x2
}
exact
mean(X1)
mean(X2)

#2019-Q2------------------------------------------------------------------------------------
temp <- sum(1/c(1:10))
c <- 1/temp  # c = 0.3414172
p <- c/c(1:10)
theta <- 1/sum(p^2) # theta = 5.535575

sim_n <- 10000; n1 <- 5; n2 <- 20

set.seed(1234)
N5 <- rmultinom(n = sim_n, size = n1, prob = p)
# row-i means Xi & column-j means j-th simulation with n = 5

g1 <- function(x, n) (x/n)^2
g2 <- function(x, n) x*(x-1)/(n*(n-1))

theta1_n1 <- apply(N5, 2, function(x) 1/sum(g1(x,n1)))
theta2_n1 <- apply(N5, 2, function(x) 1/sum(g2(x,n1)))
theta1_n1 <- theta1_n1[is.finite(theta1_n1)]
theta2_n1 <- theta2_n1[is.finite(theta2_n1)]

theta1_bar_n1 <- mean(theta1_n1)
theta2_bar_n1 <- mean(theta2_n1)
bias1_n1 <- theta1_bar_n1 - theta
bias2_n1 <- theta2_bar_n1 - theta
sim1_n1 <- length(theta1_n1)
sim2_n1 <- length(theta2_n1)
SE1_n1<-sqrt(sum((theta1_n1 - theta1_bar_n1)^2)/sim1_n1)
SE2_n1<-sqrt(sum((theta2_n1 - theta2_bar_n1)^2)/sim2_n1)
R_mse1_n1<-sqrt(sum((theta1_n1 - theta)^2)/sim1_n1)
R_mse2_n1<-sqrt(sum((theta2_n1 - theta)^2)/sim2_n1)

N20 <- rmultinom(n = sim_n,size = n2,prob = p)
# row-i means Xi & column-j means jth simulation with n = 20

theta1_n2 <- apply(N20, 2, function(x) 1/sum(g1(x,n2)))
theta2_n2 <- apply(N20, 2, function(x) 1/sum(g2(x,n2)))
theta1_n2 <- theta1_n2[is.finite(theta1_n2)]
theta2_n2 <- theta2_n2[is.finite(theta2_n2)]

theta1_bar_n2 <- mean(theta1_n2)
theta2_bar_n2 <- mean(theta2_n2)
bias1_n2 <- theta1_bar_n2 - theta
bias2_n2 <- theta2_bar_n2 - theta
sim1_n2 <- length(theta1_n2)
sim2_n2 <- length(theta2_n2)
SE1_n2<-sqrt(sum((theta1_n2 - theta1_bar_n2)^2)/sim1_n2)
SE2_n2<-sqrt(sum((theta2_n2 - theta2_bar_n2)^2)/sim2_n2)
R_mse1_n2<-sqrt(sum((theta1_n2 - theta)^2)/sim1_n2)
R_mse2_n2<-sqrt(sum((theta2_n2 - theta)^2)/sim2_n2)

result<-{data.frame("n" = c(n1," "," ",n2," "," "),
                    "Mean estimate" = c(" ",theta1_bar_n1,theta2_bar_n1," ",theta1_bar_n2,theta2_bar_n2),
                    "Bias" = c(" ",bias1_n1,bias2_n1," ",bias1_n2,bias2_n2),
                    "Sample SE" = c(" ",SE1_n1,SE2_n1," ",SE1_n2,SE2_n2),
                    "Root MSE" = c(" ",R_mse1_n1,R_mse2_n1," ",R_mse1_n2,R_mse2_n2),
                    row.names = c(" ","theta1_hat","theta2_hat","  ","theta1-hat","theta2-hat"))
}
result
# write.csv(result,"C:/Users/**Username**/Desktop/result.csv")

#2019-Q3------------------------------------------------------------------------------------
X <- c(11.2,10.9,9.2, 10.5, 7.1, 6.9,11.5,5.4,7.8, 10.4)
as.vector(X)

L_f <- function(theta){
  f <- exp(-(X-theta))/(1+exp(-(X-theta)))^2
  lnL <- sum(log(f))
  return(lnL)
}
p <- optimize(f = L_f,c(-20,20), maximum = TRUE)
p