theta <- seq(0.01, 0.99, 0.01)
x <- c(1997, 907, 904, 32)
### The Gradient of The Log-likelihood
Grad <- function (x, theta) {
dl <- x[1]/(2+theta)-(x[2]+x[3])/(1-theta)+x[4]/theta
}
### The Hessian of The Log-likelihood
Hess <- function (x, theta) {
ddl <- -(x[1]/(2+theta)^2+(x[2]+x[3])/(1-theta)^2+x[4]/theta^2)
}
### The Log-likelihood
Like <- function (x, theta) {
l <- x[1]*log(2+theta)+(x[2]+x[3])*log(1-theta)+x[4]*log(theta)
}
### Secant Method for The Multinomial
Secant <- function (theta0, theta1, maxit) {
header <- paste0("Iteration", " Theta", " MRE", " Gradient")
print(header)
for (it in 1:maxit){
theta2 <- theta1-Grad(x, theta1)*(theta1-theta0)/(Grad(x, theta1)-Grad(x,
theta0))
Mre <- abs(theta2- theta1)/max(1, abs(theta2))
if((Grad(x, theta2) < tolgrad) && (Mre < tolerr)) break
print(sprintf(" %2d %.12f %2.1e %2.1e", it, theta2, Mre,
Grad(x, theta2)))
theta0 <- theta1
theta1 <- theta2
}
}
### Run Secant
x <- c(1997, 907, 904, 32)
tolgrad <- 1e-9
tolerr <- 1e-6
Secant(theta0=0.02, theta1=0.01, maxit=20)## [1] "Iteration Theta MRE Gradient"
## [1] " 1 0.024561848011 1.5e-02 4.3e+02"
## [1] " 2 0.027823212918 3.3e-03 2.7e+02"
## [1] " 3 0.033351047002 5.5e-03 6.8e+01"
## [1] " 4 0.035197558267 1.8e-03 1.3e+01"
## [1] " 5 0.035646185180 4.5e-04 7.9e-01"
## [1] " 6 0.035674309037 2.8e-05 9.6e-03"
## [1] " 7 0.035674655571 3.5e-07 6.9e-06"Secant2 <- function (theta0, theta1, maxit) {
header <- paste0("Iteration", " Theta", " Ratio", " Digits")
print(header)
print(sprintf(' %2d %.12f -- --', 0, theta1))
for (it in 1:maxit){
theta2 <- theta1-Grad(x, theta1)*(theta1-theta0)/(Grad(x, theta1)-Grad(x,
theta0))
Mre <- abs(theta2- theta1)/max(1, abs(theta2))
ratio = abs(theta2-tstar)/(abs(theta1-tstar)^((1+sqrt(5))/2))
digits = -log10(abs(theta2-tstar)/abs(tstar))
if((Grad(x, theta2) < tolgrad) && (Mre < tolerr)) break
print(sprintf(" %2d %.12f %.3e %.3f", it, theta2, ratio,
digits))
theta0 <- theta1
theta1 <- theta2
}
}
x <- c(1997, 907, 904, 32)
tolgrad <- 1e-9
tolerr <- 1e-6
tstar <- -1657/7680 + sqrt(3728689)/7680
Secant2(theta0=0.02, theta1=0.01, maxit=20)## [1] "Iteration Theta Ratio Digits"
## [1] " 0 0.010000000000 -- --"
## [1] " 1 0.024561848011 4.162e+00 0.507"
## [1] " 2 0.027823212918 1.140e+01 0.657"
## [1] " 3 0.033351047002 5.918e+00 1.186"
## [1] " 4 0.035197558267 8.715e+00 1.874"
## [1] " 5 0.035646185180 6.738e+00 3.098"
## [1] " 6 0.035674309037 7.852e+00 5.012"
## [1] " 7 0.035674655571 7.135e+00 8.151"The convergence of a sequence \(\{x^{(n)}\}\) to \(x^{*}\) is of order \(\alpha\), if \(\displaystyle \lim_{n \to
\infty}\frac{|x^{(n+1)}-x^{*}|}{|x^{(n)}-x^{*}|^\alpha}=C\),
where C is a constant.
In this example, \(\alpha\) is (1+\(\sqrt{5}\))/2 (1<\(\alpha\)<2) and the convergence ratio
converges to
a constant (around 7.xx), so the secant method locally converges
super-linearly.
If we put \(\alpha\)=2 and run the
secant method, it does not converge to a constant.
This shows that this algorithm does not converge quadratically.
Secant2 <- function (theta0, theta1, maxit) {
header <- paste0("Iteration", " Theta", " Ratio", " Digits")
print(header)
print(sprintf(' %2d %.12f -- --', 0, theta1))
for (it in 1:maxit){
theta2 <- theta1-Grad(x, theta1)*(theta1-theta0)/(Grad(x, theta1)-Grad(x,
theta0))
Mre <- abs(theta2- theta1)/max(1, abs(theta2))
ratio = abs(theta2-tstar)/(abs(theta1-tstar)^((1+sqrt(5))/2))
digits = -log10(abs(theta2-tstar)/abs(tstar))
if((Grad(x, theta2) < tolgrad) && (Mre < tolerr)) break
print(sprintf(" %2d %.12f %.3e %.3f", it, theta2, ratio,
digits))
theta0 <- theta1
theta1 <- theta2
}
}
x <- c(1997, 907, 904, 32)
tolgrad <- 1e-9
tolerr <- 1e-6
tstar <- -1657/7680 + sqrt(3728689)/7680
Secant2(theta0=0.08, theta1=0.09, maxit=10)## [1] "Iteration Theta Ratio Digits"
## [1] " 0 0.090000000000 -- --"
## [1] " 1 -0.006087947626 4.652e+00 -0.068"
## [1] " 2 0.102137650875 1.133e+01 -0.270"
## [1] " 3 0.117525594044 6.579e+00 -0.361"
## [1] " 4 -0.037438751778 4.195e+00 -0.312"
## [1] " 5 0.291354660240 1.761e+01 -0.855"
## [1] " 6 60.370606350928 5.482e+02 -3.228"
## [1] " 7 58.056930549586 7.631e-02 -3.211"
## [1] " 8 118.921236872163 1.666e-01 -3.523"
## [1] " 9 177.489859174286 7.789e-02 -3.697"
## [1] " 10 296.948377158938 6.816e-02 -3.920"