SC_HW5_01

Exercise J-4.1

(a)

# EM Function
my_EM <- function(theta, maxit, tolerr){
  header = paste0("Iteration", "      Theta", "                     Mre")
  print(header)
  print("---------------------------------------------------")
  print(sprintf('%02.0f            %.12f             --', 0, theta))
  for (it in 1:maxit){
    x2star = 1997*theta/(2+theta)
    theta.n=(x2star+32)/(x2star+907+904+32)
    Mre = abs(theta.n - theta)/max(1, abs(theta.n))
    print(sprintf('%02.0f            %.12f           %2.1e',
                  it, theta.n, Mre))
    if (Mre < tolerr){ break }
    theta=theta.n
  }
}

my_EM(theta=0.02, maxit=200, tolerr=1e-6)

## [1] "Iteration      Theta                     Mre"
## [1] "---------------------------------------------------"
## [1] "00            0.020000000000             --"
## [1] "01            0.027793132773           7.8e-03"
## [1] "02            0.031742941132           3.9e-03"
## [1] "03            0.033721123764           2.0e-03"
## [1] "04            0.034705946241           9.8e-04"
## [1] "05            0.035194771667           4.9e-04"
## [1] "06            0.035437045211           2.4e-04"
## [1] "07            0.035557033560           1.2e-04"
## [1] "08            0.035616437341           5.9e-05"
## [1] "09            0.035645841640           2.9e-05"
## [1] "10            0.035660395186           1.5e-05"
## [1] "11            0.035667598091           7.2e-06"
## [1] "12            0.035671162906           3.6e-06"
## [1] "13            0.035672927162           1.8e-06"
## [1] "14            0.035673800302           8.7e-07"

(b)

# Check Convergence
Check_con <- function(theta, maxit, tolerr){
  header = paste0("It", "     Theta", "           Mre       Ratio_2      Ratio_G      Ratio_1")
  print(header)
  print("----------------------------------------------------------------------")
  print(sprintf('%02.0f   %.12f     --          --          --          --', 0, theta))
  for (it in 1:maxit){
    x2star = 1997*theta/(2+theta)
    theta.n=(x2star+32)/(x2star+907+904+32)
    Mre = abs(theta.n - theta)/max(1, abs(theta.n))
    tstar = (-1657+sqrt(3728689))/7680
    Ratio_2 = abs(theta.n-tstar)/(abs(theta-tstar)^(2))
    Ratio_G = abs(theta.n-tstar)/(abs(theta-tstar)^((1+sqrt(5))/2))
    Ratio_1 = abs(theta.n-tstar)/abs(theta-tstar)
    print(sprintf('%02.0f   %.12f   %2.1e   %.3e    %.3e    %.3e',
                  it, theta.n, Mre, Ratio_2, Ratio_G, Ratio_1))
    if (Mre < tolerr){ break }
    theta=theta.n
  }
}

Check_con(theta=0.02, maxit=200, tolerr=1e-6)

## [1] "It     Theta           Mre       Ratio_2      Ratio_G      Ratio_1"
## [1] "----------------------------------------------------------------------"
## [1] "00   0.020000000000     --          --          --          --"
## [1] "01   0.027793132773   7.8e-03   3.208e+01    6.559e+00    5.028e-01"
## [1] "02   0.031742941132   3.9e-03   6.329e+01    9.953e+00    4.989e-01"
## [1] "03   0.033721123764   2.0e-03   1.264e+02    1.524e+01    4.969e-01"
## [1] "04   0.034705946241   9.8e-04   2.538e+02    2.343e+01    4.959e-01"
## [1] "05   0.035194771667   4.9e-04   5.114e+02    3.611e+01    4.954e-01"
## [1] "06   0.035437045211   2.4e-04   1.032e+03    5.571e+01    4.951e-01"
## [1] "07   0.035557033560   1.2e-04   2.083e+03    8.599e+01    4.950e-01"
## [1] "08   0.035616437341   5.9e-05   4.208e+03    1.328e+02    4.950e-01"
## [1] "09   0.035645841640   2.9e-05   8.501e+03    2.051e+02    4.949e-01"
## [1] "10   0.035660395186   1.5e-05   1.718e+04    3.167e+02    4.949e-01"
## [1] "11   0.035667598091   7.2e-06   3.470e+04    4.891e+02    4.949e-01"
## [1] "12   0.035671162906   3.6e-06   7.012e+04    7.555e+02    4.949e-01"
## [1] "13   0.035672927162   1.8e-06   1.417e+05    1.167e+03    4.949e-01"
## [1] "14   0.035673800302   8.7e-07   2.863e+05    1.802e+03    4.949e-01"

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 algorithm, when we put \(\alpha\) = 2, (1+\(\sqrt{5}\))/2, 1, we can obtain the Convergence Ratios like the above result. When \(\alpha\) = 2, (1+\(\sqrt{5}\))/2, the ratio doesn’t converge. However, when \(\alpha\) = 1, the ratio converges to a constant(4.949e-01). Therefore, this algorithm is linearly convergent.