SC_HW2_01
Myung
2/23/2022
Exercise GH-2.2
(a) Graph the Log-likelihood Function (-pi<theta<pi)
x <- c(3.91, 4.85, 2.28, 4.06, 3.7, 4.04, 5.46, 3.53, 2.28, 1.96, 2.53, 3.88,
2.22, 3.47, 4.82, 2.46, 2.99, 2.54, .52, 2.5)
theta <- seq(-pi,pi, length.out=1000)
L <- function(x, theta) {
l <- matrix(NA, length(x),length(theta))
for (i in 1:length(x)){
l[i,] <- log((1-cos(x[i]-theta))/(2*pi))
}
colSums(l)
}
plot(theta, L(x, theta), type="l", lwd=3, main="The Log-likelihood function",
xlab="theta", ylab="Log-likelihood")
(b) The method-of-moment estimator of theta
\[\mu_1^\prime=E(X_1)=\int_{0}^{2\pi}x\frac{[1-cos(x-\theta)]}{2\pi}dx=\pi+sin\theta\] \[m_1 = \frac{\sum_{i=1}^{n}x_i}{n}\] \[\mu_1^\prime=m_1\] Therefore, \[\tilde{\theta}=sin^{-1}(\frac{\sum_{i=1}^{n}x_i}{n}-\pi)\]
MME_theta <- asin(mean(x)-pi)
MME_theta## [1] 0.05844061(c) MLE using the Newton-Raphson method
### The Gradient of the Log-likelihood
Grad <- function(x, theta){
dl <- matrix(NA, length(x),length(theta))
for (i in 1:length(x)){
dl[i,] <- (-sin(x[i]-theta))/(1-cos(x[i]-theta))
}
colSums(dl)
}
### The Hessian of The Log-likelihood
Hess <- function(x, theta){
ddl <- matrix(NA, length(x),length(theta))
for (i in 1:length(x)){
ddl[i,] <- cos(x[i]-theta)-1/(1-cos(x[i]-theta))^2
}
colSums(ddl)
}
### Run Newton, using MME (Convergence Criteria)
maxit <- 1000
x <- c(3.91, 4.85, 2.28, 4.06, 3.7, 4.04, 5.46, 3.53, 2.28, 1.96, 2.53, 3.88,
2.22, 3.47, 4.82, 2.46, 2.99, 2.54, .52, 2.5)
tolgrad <- 1e-9
tolerr <- 1e-6
Newton <- function(maxit, x, theta){
header = paste0("Iteration"," Theta"," Ratio"," Gradient" )
print(header)
print(sprintf("%03d %.12f -- %2.1e", 0, theta, Grad(x,
theta)))
for (it in 1:maxit){
dl <- Grad(x, theta)
ddl <- Hess(x, theta)
tnew <- theta - dl/ddl
Mre <- abs(tnew- theta)/max(1, abs(tnew))
if((abs(Grad(x, tnew)) < tolgrad) && (Mre < tolerr)) break
print(sprintf("%03d %.12f %2.1e %2.1e", it, tnew, Mre, Grad(x,
tnew)))
theta <- tnew
}
}
theta <- MME_theta; Newton(maxit, x, theta) ## [1] "Iteration Theta Ratio Gradient"
## [1] "000 0.058440606140 -- -1.6e+00"
## [1] "001 0.044706428198 1.4e-02 -1.3e+00"
## [1] "002 0.032775338694 1.2e-02 -1.0e+00"
## [1] "003 0.022755660475 1.0e-02 -7.8e-01"
## [1] "004 0.014591631009 8.2e-03 -5.9e-01"
## [1] "005 0.008108080713 6.5e-03 -4.5e-01"
## [1] "006 0.003065655631 5.0e-03 -3.3e-01"
## [1] "007 -0.000791752455 3.9e-03 -2.5e-01"
## [1] "008 -0.003705304767 2.9e-03 -1.8e-01"
## [1] "009 -0.005884815879 2.2e-03 -1.3e-01"
## [1] "010 -0.007503479774 1.6e-03 -9.8e-02"
## [1] "011 -0.008699183102 1.2e-03 -7.2e-02"
## [1] "012 -0.009578953697 8.8e-04 -5.3e-02"
## [1] "013 -0.010224385414 6.5e-04 -3.8e-02"
## [1] "014 -0.010696887445 4.7e-04 -2.8e-02"
## [1] "015 -0.011042252581 3.5e-04 -2.0e-02"
## [1] "016 -0.011294401772 2.5e-04 -1.5e-02"
## [1] "017 -0.011478341222 1.8e-04 -1.1e-02"
## [1] "018 -0.011612441091 1.3e-04 -7.9e-03"
## [1] "019 -0.011710162454 9.8e-05 -5.8e-03"
## [1] "020 -0.011781351073 7.1e-05 -4.2e-03"
## [1] "021 -0.011833198782 5.2e-05 -3.1e-03"
## [1] "022 -0.011870953762 3.8e-05 -2.2e-03"
## [1] "023 -0.011898443133 2.7e-05 -1.6e-03"
## [1] "024 -0.011918456307 2.0e-05 -1.2e-03"
## [1] "025 -0.011933025599 1.5e-05 -8.6e-04"
## [1] "026 -0.011943631316 1.1e-05 -6.2e-04"
## [1] "027 -0.011951351478 7.7e-06 -4.5e-04"
## [1] "028 -0.011956971031 5.6e-06 -3.3e-04"
## [1] "029 -0.011961061462 4.1e-06 -2.4e-04"
## [1] "030 -0.011964038816 3.0e-06 -1.8e-04"
## [1] "031 -0.011966205959 2.2e-06 -1.3e-04"
## [1] "032 -0.011967783358 1.6e-06 -9.3e-05"
## [1] "033 -0.011968931495 1.1e-06 -6.8e-05"
## [1] "034 -0.011969767182 8.4e-07 -4.9e-05"
## [1] "035 -0.011970375446 6.1e-07 -3.6e-05"
## [1] "036 -0.011970818178 4.4e-07 -2.6e-05"
## [1] "037 -0.011971140424 3.2e-07 -1.9e-05"
## [1] "038 -0.011971374974 2.3e-07 -1.4e-05"
## [1] "039 -0.011971545693 1.7e-07 -1.0e-05"
## [1] "040 -0.011971669952 1.2e-07 -7.3e-06"
## [1] "041 -0.011971760395 9.0e-08 -5.3e-06"
## [1] "042 -0.011971826224 6.6e-08 -3.9e-06"
## [1] "043 -0.011971874139 4.8e-08 -2.8e-06"
## [1] "044 -0.011971909014 3.5e-08 -2.1e-06"
## [1] "045 -0.011971934397 2.5e-08 -1.5e-06"
## [1] "046 -0.011971952873 1.8e-08 -1.1e-06"
## [1] "047 -0.011971966321 1.3e-08 -7.9e-07"
## [1] "048 -0.011971976109 9.8e-09 -5.8e-07"
## [1] "049 -0.011971983233 7.1e-09 -4.2e-07"
## [1] "050 -0.011971988419 5.2e-09 -3.1e-07"
## [1] "051 -0.011971992193 3.8e-09 -2.2e-07"
## [1] "052 -0.011971994940 2.7e-09 -1.6e-07"
## [1] "053 -0.011971996940 2.0e-09 -1.2e-07"
## [1] "054 -0.011971998395 1.5e-09 -8.6e-08"
## [1] "055 -0.011971999454 1.1e-09 -6.2e-08"
## [1] "056 -0.011972000225 7.7e-10 -4.5e-08"
## [1] "057 -0.011972000787 5.6e-10 -3.3e-08"
## [1] "058 -0.011972001195 4.1e-10 -2.4e-08"
## [1] "059 -0.011972001492 3.0e-10 -1.7e-08"
## [1] "060 -0.011972001709 2.2e-10 -1.3e-08"
## [1] "061 -0.011972001866 1.6e-10 -9.3e-09"
## [1] "062 -0.011972001981 1.1e-10 -6.7e-09"
## [1] "063 -0.011972002064 8.3e-11 -4.9e-09"
## [1] "064 -0.011972002125 6.1e-11 -3.6e-09"
## [1] "065 -0.011972002169 4.4e-11 -2.6e-09"
## [1] "066 -0.011972002201 3.2e-11 -1.9e-09"
## [1] "067 -0.011972002225 2.3e-11 -1.4e-09"
## [1] "068 -0.011972002242 1.7e-11 -1.0e-09"=> Therefore, \(\hat{\theta}\)(MLE) = -0.011972
### Solutions (starting at -2.7, 2.7)
Newton <- function(maxit, x, theta){
for (it in 1:maxit){
dl <- Grad(x, theta)
ddl <- Hess(x, theta)
tnew <- theta - dl/ddl
Mre <- abs(tnew- theta)/max(1, abs(tnew))
if(abs(Grad(x, tnew) < tolgrad) && (Mre < tolerr)) break
theta <- tnew
}
return(theta)
}
### Start at -2.7
theta <- -2.7 ; Newton(maxit, x, theta) ## [1] -2.667439### Start at 2.7
theta <- 2.7 ; Newton(maxit, x, theta) ## [1] 2.873094(d) Graph of Solutions
Newton <- function(maxit, x, theta){
for (it in 1:maxit){
dl <- Grad(x, theta)
ddl <- Hess(x, theta)
tnew <- theta - dl/ddl
Mre <- abs(tnew- theta)/max(1, abs(tnew))
if(abs(Grad(x, tnew) < tolgrad) && (Mre < tolerr)) break
theta <- tnew
}
return(theta)
}
maxit <- 1000
x <- c(3.91, 4.85, 2.28, 4.06, 3.7, 4.04, 5.46, 3.53, 2.28, 1.96, 2.53, 3.88,
2.22, 3.47, 4.82, 2.46, 2.99, 2.54, .52, 2.5)
tolgrad <- 1e-9
tolerr <- 1e-6
theta_ini <- seq(-pi,pi, length.out=200)
theta <- theta_ini
theta_con <- Newton(maxit, x, theta)
My_theta <- cbind(theta_ini, theta_con)
#My_theta
plot(My_theta[,1], My_theta[,2], type="p", xlim=c(-pi,pi), ylim=c(-pi, pi),
xlab="Initial Value", ylab="Converged Value", main="The Graph of Solutions",
cex=0.7)