StatComp_2017Q1Q2

#2017-Q1------------------------------------------------------------------------------------
# X1~exp(5); X2~trunc exp; X3~Weibull
exact <- 5+5*exp(-5)*(exp(5)-1)/(1-6*exp(-5))+2*gamma(1+1/2)

f2 <- function(x){
  exp(-5)*(5^x)/(factorial(x)*(1-6*exp(-5)))
}

X <- vector()
set.seed(1234)
for(i in 1:10000){
  x <- 0
  while(TRUE){
    u <- runif(1)
    if(u<1/3){
      x <- x+rexp(1,rate=1/5)
      break
    }
    else{
      if(u<2/3){
        # simulation x2
        u1 <- runif(1)
        x2 <- 2
        p <- f2(2)
        while(u1>p){
          x2 <- x2+1
          p <- p+f2(x2)
        }
        x <- x+x2
      }
      else{
        x <- x + rweibull(1,shape=2, scale = 2)
      }
    }
  }
  X[i] <- x
}
exact # 11.948
mean(X) # 12.05248

#2017-Q2------------------------------------------------------------------------------------
X <- c(9.853093,12.993558,11.294747,11.171975,8.913325,10.976033,12.353118)
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