## Program the Newton Raphson algorithm for a numerical computation of the
## ML estimate theta_hat of the parameter θ = P(Coin = HEAD)
## in our coin toss example of Chapter 6 of our script.
## Replicate the results in Table 6.1. of our script.
## Below I use the same data that was used to produce the
## results in Table 6.1 of our script. However, you
## can produce new data by setting another seed-value
theta_true <- 0.2# unknown true theta valuen<-5 # sample size## Use a common Random Number Generator:RNGkind(sample.kind = “Rounding”) # RNGkind()是random number generation## Warning in RNGkind(sample.kind = “Rounding”): non-uniform ‘Rounding’ sampler used set.seed(1)# ?sample# simulate data: n many (unfair) coin tossesx <- sample(x = c(0,1),size = n,replace = TRUE,prob = c(1-theta_true, theta_true))## number of heads (i.e., the number of “1”s in x)h <- sum(x)## First derivative of the log-likelihood functionLp_fct <- function(theta, h = h, n = n){(h/theta) – (n – h)/(1 – theta)}## Second derivative of the log-likelihood functionLpp_fct <- function(theta, h = h, n = n){- (h/theta^2) – (n – h)/(1 – theta)^2}t <- 1e-10 # convergence criterion 是1*10的-10次方check <- TRUE # for stopping the while-loopi <- 0 # count iterationstheta <- 0.4 # starting valueLp<- Lp_fct(theta, h=h, n=n)Lpp <- Lpp_fct(theta, h=h, n=n)while(check){i <- i + 1##theta_new <- theta[i] – (Lp_fct(theta[i], h=h, n=n) / Lpp_fct(theta[i], h=h, n=n))Lp_new <- Lp_fct(theta_new, h = h, n = n)Lpp_new <- Lpp_fct(theta_new, h = h, n = n)##theta <- c(theta, theta_new)Lp <- c(Lp, Lp_new)Lpp <- c(Lpp, Lpp_new)##if(abs(Lp_fct(theta_new, h=h, n=n)) < t ){check <- FALSE}}cbind(theta, Lp, Lp/Lpp)#>thetaLp
#> [1,] 0.4000000 -4.166667e+002.400000e-01
#> [2,] 0.16000001.488095e+00 -3.326733e-02
#> [3,] 0.19326732.159084e-01 -6.558924e-03
#> [4,] 0.19982635.433195e-03 -1.736356e-04
#> [5,] 0.19999993.539786e-06 -1.132731e-07
#> [6,] 0.20000001.504574e-12 -4.814638e-14
Reviews
There are no reviews yet.