## 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 value
n <- 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 tosses
x <- 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 function
Lp_fct <- function(theta, h = h, n = n){
(h/theta) - (n - h)/(1 - theta)
}
## Second derivative of the log-likelihood function
Lpp_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-loop
i <- 0 # count iterations
theta <- 0.4 # starting value
Lp <- 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)
#> theta Lp
#> [1,] 0.4000000 -4.166667e+00 2.400000e-01
#> [2,] 0.1600000 1.488095e+00 -3.326733e-02
#> [3,] 0.1932673 2.159084e-01 -6.558924e-03
#> [4,] 0.1998263 5.433195e-03 -1.736356e-04
#> [5,] 0.1999999 3.539786e-06 -1.132731e-07
#> [6,] 0.2000000 1.504574e-12 -4.814638e-14