Newton-Rapshon, Fisher-Scoring, and EM Algorithm

Author

Fatkhurokhman Fauzi

Newton-Rapshon

Metode Newton-rapshon adalah salah satu metode numerik untuk mencari akar suatu fungsi non-linear \(f(x)=0\). Metode ini termasuk dalam kategori metode iteratif yang berbasis pada pendekatan turunan pertama (derivatif) dari fungsi tersebut.

\[ x_{i+1} = x_i - \frac{f(x_i)}{f'(x_i)} \]

\(x_i\) adalah perkiraan nilai akar fungsi saat iterasi ke-i

\(x_{i+1}\) adalah updataing perkiraan nilai akar fungsi untuk iterasi berikutnya

\(f(x_i)\) adalah fungsi \(f(x)\) yang dievaluasi pada titik \(x_i\)

\(f'(x_i)\) adalah nilai turunan pertama fungsi \(f(x_i)\) yang dievaluasi pada titik \(x_i\)

Newton-Rapshon

Diketahui \(f(x)=x^2-2\) carilah niali \(x\) menggunakan newton-rapshon

Mendefinisikan fungsi \(f(x)\)
```
fungsi1<-function(x) x^2-2
```
Mendefinisikan turunan pertama \(f'(x)\)
```
fungsiturunan<-function (x) 2*x
```

Membuat fungsi newton-rapshon

newton_rapshon <- function(f, f_prime, x0, tol = 1e-6, max_iter = 100) {
  x <- x0
  for (i in 1:max_iter) {
    x_new <- x - f(x) / f_prime(x)
    cat("Iterasi ke:", i, "x baru:", x_new, "\n")
    if (abs(x_new - x) < tol) {
      return(x_new)
    }
    x <- x_new
  }
  return(x)
}

Menjalankan fungsi newton-rapshon

x0 <- 1
akar <- newton_rapshon(fungsi1, fungsiturunan, x0)

Iterasi ke: 1 x baru: 1.5 
Iterasi ke: 2 x baru: 1.416667 
Iterasi ke: 3 x baru: 1.414216 
Iterasi ke: 4 x baru: 1.414214 
Iterasi ke: 5 x baru: 1.414214

cat("Akar yang ditemukan:", akar, "\n")

Akar yang ditemukan: 1.414214

Newton-Rapshon untuk Optimisasi

Jika kita memiliki fungsi \(f(x)\) yang ingin kita minimalkan atau memaksimalkan, titik optimalnya dapat dicari dengan menyelesaikan \(f'(x)=0\) metode newton-rapshon memperbaharui nilai \(x\) sebagai berikut:

\[ x_{i+1}=x_i-\frac{f'(x_i)}{f''(x_i)} \]

\(x_i\) adalah perkiraan nilai akar fungsi saat iterasi ke-i

\(x_{i+1}\) adalah updataing perkiraan nilai akar fungsi untuk iterasi berikutnya

\(f'(x_i)\) adalahturunan pertama dari fungsi \(f(x)\) yang dievaluasi pada titik \(x_i\)

\(f''(x_i)\) adalah nilai turunan kedua dari fungsi \(f(x_i)\) yang dievaluasi pada titik \(x_i\)

Contoh Soal

Diketahui \(f(x)=x^2+3x+2\) carilah titik minimum menggunan newton-rapshon

Mendefinisikan turunan pertama \(f'(x)\)
```
f_prime <- function(x) 2 * x + 3
```
Mendefiniskan turunan kedua \(f''(x)\)
```
f_double_prime <- function(x) 2
```

Membuat fungsi newton-rapshon

newton_rapshon_min <- function(f_prime, f_double_prime, x0, tol = 1e-6, max_iter = 100) {
  x <- x0
  for (i in 1:max_iter) {
    x_new <- x - f_prime(x) / f_double_prime(x)
    cat("Iterasi ke:", i, "x baru:", x_new, "\n")
    if (abs(x_new - x) < tol) {
      return(x_new)
    }
    x <- x_new
  }
  return(x)
}

Menjalankan fungsi newton-rapshon

x0 <- 0
min_point <- newton_rapshon_min(f_prime, f_double_prime, x0)

Iterasi ke: 1 x baru: -1.5 
Iterasi ke: 2 x baru: -1.5

cat("Titik minimum ditemukan di x =", min_point, "\n")

Titik minimum ditemukan di x = -1.5

Newton-Rapshon Multidimensi

Optimasi multidimensi melibatkan fungsi multivariabel \(f(x_1, x_2, x_3, ...,x_p)\) sehingga nilai yang menyebabkan fungsi\(f\) maksimum atau minimum berasal dari \(p\) variabel

\[ x_{i+1}=x_i-\frac{f'(x_i)}{f''(x_i)} \]

\(x=[x_1, x_2, x_3, ...,x_p]^T\)

\(f'(x_i)\) adalah matriks khusus yang dinamakan Jacobian

\(f''(x_i)\) adalah matriks khusus yang dinamakan Hessian

Contoh

Diketahui \(f(x,y)=x^2+y^2+3x+2y+xy\) carilah titik minimul \(x,y\) menggunakan newton-rapshon

Mendefinisikan matriks Jacobian

grad_f <- function(x, y) {
  return(c(2*x + 3 + y, 2*y + 2 + x))
}

Mendefinisikan matriks Hessian

hessian_f <- function(x, y) {
  return(matrix(c(2, 1, 1, 2), nrow = 2, byrow = TRUE))
}

Membuat fungsi newton-rapshon

newton_rapshon_2d <- function(grad_f, hessian_f, x0, tol = 1e-6, max_iter = 100) {
  x <- x0
  for (i in 1:max_iter) {
    grad <- grad_f(x[1], x[2])          
    hessian <- hessian_f(x[1], x[2])    # Hitung Hessian
    hessian_inv <- solve(hessian)       # Hitung invers Hessian
    x_new <- x - hessian_inv %*% grad   # Update nilai x
    cat("Iterasi ke:", i, "x baru:", x_new, "\n")
    if (sum(abs(x_new - x)) < tol) {
      return(x_new)
    }
    x <- x_new
  }
  return(x)
}

Menjalankan fungsi newton-rapshon

x0 <- c(0, 0)

# Panggil fungsi Newton-rapshon untuk menemukan titik minimum
min_point <- newton_rapshon_2d(grad_f, hessian_f, x0)

Iterasi ke: 1 x baru: -1.333333 -0.3333333 
Iterasi ke: 2 x baru: -1.333333 -0.3333333

cat("Titik minimum ditemukan di (x, y) =", min_point[1], ",", min_point[2], "\n")

Titik minimum ditemukan di (x, y) = -1.333333 , -0.3333333

Newton-Rapshon untuk MLE

Formula newton-rapshon untuk MLE adalah sebagai berikut:

\[\begin{align*}\theta^{(i+1)} &= \theta^{(i)} - \left[ H_f(\theta^{(i)}) \right]^{-1} j_f(\theta^{(i)}) \\ \theta^{(i+1)} &= \theta^{(i)} - \left[ -I(\theta^{(i)}) \right]^{-1} S(\theta^{(i)}) \\ \theta^{(i+1)} &= \theta^{(i)} + \left[ I(\theta^{(i)}) \right]^{-1} S(\theta^{(i)}) \end{align*}\]

dengan

\(\theta=[\theta_1,\theta_2, \theta_3, ...,\theta_p]^T\)

\(J_f(\theta^{(i)})\) adalah matriks jacobian untuk perkiraan parameter dugaan \(\theta^{(i)}\)

\(H_f(\theta^{(i)})\) adalah matriks Hessian untuk perkiraan parameter dugaan \(\theta^{(i)}\)

\(S(\theta^{(i)})\) adalah Score Function untuk perkiraan parameter dugaan \(\theta^{(i)}\)

\(I(\theta^{(i)})\) adalah Observed Fisher Information untuk perkiraan parameter dugaan \(\theta^{(i)}\)

Contoh soal

Misalkan kita memiliki data sampel \(X_1, X_2, X_3,...,X_n\) yang berasal dari distribusi normal dengan parameter \(\mu\) dan \(\sigma^2\) . fungsi kepadatan peluang sebagai berikut:

\[f(x) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left( -\frac{(x - \mu)^2}{2\sigma^2} \right)\]

dengan data bangkitan distribusi normal, \(\mu\) = 10, \(\sigma^2\) =2. Estimasi menggunakan newton-rapshon dan fisher scoring

Membangkitkan data berdistribusi normal dengan \(\mu\) = 10, \(\sigma^2\) =2

set.seed(123)
mu_true <- 10
sigma2_true <- 2
n <- 100

data_sample <- rnorm(n, mean = mu_true, sd = sqrt(sigma2_true))

Membuat fungsi matriks Score function

score_function <- function(theta, x) {
  mu <- theta[1]
  sigma2 <- theta[2]
  n <- length(x)
  g1 <- sum(x - mu) / sigma2
  g2 <- -n / (2 * sigma2) + sum((x - mu)^2) / (2 * sigma2^2)
  return(c(-g1, -g2))
}

Membuat fungsi matriks Observed Fisher Information

fisher_information <- function(theta, x) {
  mu <- theta[1]
  sigma2 <- theta[2]
  n <- length(x)
  h11 <- n / sigma2
  h12 <- sum((x - mu)) / (sigma2^2)
  h22 <- - n / (2 * sigma2^2) + sum((x - mu)^2) / sigma2^3
  return(matrix(c(h11, h12, h12, h22), nrow = 2, byrow = TRUE))
}

Membuat fungsi Newton-Rapshon

solve_safe <- function(matrix) {
  if (det(matrix) == 0) {
    diag(matrix) <- diag(matrix) + 1e-6  # Tambahkan nilai kecil ke diagonal
  }
  return(solve(matrix))
}

newton_raphson <- function(init_theta, x, tol = 1e-6, max_iter = 100) {
  theta <- init_theta
  for (i in 1:max_iter) {
    grad <- score_function(theta, x)
    hess <- fisher_information(theta, x)

    theta_new <- theta - (solve(hess)%*%grad)
    cat("Iterasi ke:", i, "mu:", theta_new[1,1],"simga^2:", theta_new[2,1], "\n")
    if (sum(abs(theta_new - theta)) < tol) {
      break
    }
    theta <- theta_new
  }
  return(theta)
}

Mengestimasi dengan fungsi Newton-Raphson

# Inisialisasi estimasi awal
init_theta <- c(11, 2)

# Estimasi menggunakan Newton-Raphson
est_newton <- newton_raphson(init_theta, data_sample)

Iterasi ke: 1 mu: 9.382383 simga^2: 0.2904937 
Iterasi ke: 2 mu: 9.928776 simga^2: 0.3680699 
Iterasi ke: 3 mu: 10.0435 simga^2: 0.5240223 
Iterasi ke: 4 mu: 10.09386 simga^2: 0.7352231 
Iterasi ke: 5 mu: 10.11574 simga^2: 0.9971069 
Iterasi ke: 6 mu: 10.12442 simga^2: 1.279696 
Iterasi ke: 7 mu: 10.12722 simga^2: 1.514166 
Iterasi ke: 8 mu: 10.12781 simga^2: 1.629193 
Iterasi ke: 9 mu: 10.12785 simga^2: 1.649293 
Iterasi ke: 10 mu: 10.12785 simga^2: 1.649801 
Iterasi ke: 11 mu: 10.12785 simga^2: 1.649801

# Hasil estimasi
cat("Estimasi Newton-Raphson:\n")

Estimasi Newton-Raphson:

print(est_newton)

          [,1]
[1,] 10.127853
[2,]  1.649801

Fisher-Scoring

Formulasi Fisher Scoring sebagai berikut:

\[ \theta^{(i+1)} = \theta^{(i)} + J(\theta^{(i)})^{-1} S(\theta^{(i)}) \]

\[ J(\theta) = E_{\theta} (I(\theta)) = -E \left( H_f(\theta) \right) =\begin{bmatrix}E \left( \frac{\partial l(\theta | x)}{\partial^2 \theta_1} \right) & E \left( \frac{\partial l(\theta | x)}{\partial \theta_1 \partial \theta_2} \right) & \dots & E \left( \frac{\partial l(\theta | x)}{\partial \theta_1 \partial \theta_p} \right) \\E \left( \frac{\partial l(\theta | x)}{\partial \theta_2 \partial \theta_1} \right) & E \left( \frac{\partial l(\theta | x)}{\partial^2 \theta_2} \right) & \dots & E \left( \frac{\partial l(\theta | x)}{\partial \theta_2 \partial \theta_p} \right) \\\vdots & \vdots & \ddots & \vdots \\E \left( \frac{\partial l(\theta | x)}{\partial \theta_p \partial \theta_1} \right) & E \left( \frac{\partial l(\theta | x)}{\partial \theta_p \partial \theta_2} \right) & \dots & E \left( \frac{\partial l(\theta | x)}{\partial^2 \theta_p} \right)\end{bmatrix} \]

\(\theta=[\theta_1,\theta_2, \theta_3, ...,\theta_p]^T\)

\(S(\theta^{(i)})\) adalah Score Function untuk perkiraan parameter dugaan \(\theta^{(i)}\)

\(J(\theta^{(i)})\) adalah Expexted Fisher Information untuk perkiraan parameter dugaan \(\theta^{(i)}\)

Membuat Expexted Fisher Information

fisher_info <- function(theta, x) {
  sigma2 <- theta[2]
  n <- length(x)
  I11 <- n / sigma2
  I12 <- 0
  I22 <- n / (2 * sigma2^2)
  return(matrix(c(I11, I12, I12, I22), nrow = 2, byrow = TRUE))
}

Membuat fungsi Fisher-Scoring

fisher_scoring <- function(init_theta, x, tol = 1e-6, max_iter = 100) {
  theta <- init_theta
  for (i in 1:max_iter) {
    grad <- score_function(theta, x)
    fisher <- fisher_info(theta, x)

    theta_new <- theta - solve_safe(fisher) %*% grad
    cat("Iterasi ke:", i, "mu:", theta_new[1,1],"simga^2:", theta_new[2,1], "\n")
    if (sum(abs(theta_new - theta)) < tol) {
      break
    }
    theta <- theta_new
  }
  return(theta)
}

Menjalankan fungsi Fisher-Scoring

# Estimasi menggunakan Fisher Scoring
est_fisher <- fisher_scoring(init_theta, data_sample)

Iterasi ke: 1 mu: 10.12785 simga^2: 2.410441 
Iterasi ke: 2 mu: 10.12785 simga^2: 1.649801 
Iterasi ke: 3 mu: 10.12785 simga^2: 1.649801

cat("\nEstimasi Fisher Scoring:\n")


Estimasi Fisher Scoring:

print(est_fisher)

          [,1]
[1,] 10.127853
[2,]  1.649801

EM-Algorithm

Didapat nilai Expectation

\[ E_{\theta^{(0)}} \left[ z_2 \mid \theta^{(0)}, y \right] = y_1 \left( \frac{\frac{\theta^{(0)}}{4}}{\frac{1}{2} + \frac{\theta^{(0)}}{4}} \right) \]

Didapat nilai Maximization

\[ \theta^{(1)} = \frac{\hat{z}_2 + y_4}{\hat{z}_2 + y_4 + y_2 + y_3} \]

Fungsi Expectation

expectation <- function(theta){
  100*( ( (1/4)*theta ) / ( (1/2) + (1/4)*theta) )
}

Fungsi Maximization

maximization <- function(z2){
  (z2 + y4) / (z2 + y2 + y3 + y4)
}

Menjalankan fungsi

z2=0
y2<- 12
y3<- 18
y4<- 29

niter <- 100
theta0 <- 2
save_iter <- data.frame("iter"=0,"theta"=theta0,"z2"=z2)
for (i in 1:niter){
  x2 <- expectation(theta0)
  theta <- maximization(x2)
  criteria <- abs(theta-theta0)
  theta0 <- theta
  save_iter <- rbind(save_iter,c(i,theta0,x2))
  if (criteria<10^-7){
    break
  }
}

save_iter

   iter     theta       z2
1     0 2.0000000  0.00000
2     1 0.7247706 50.00000
3     2 0.6495300 26.59933
4     3 0.6407827 24.51491
5     4 0.6397040 24.26488
6     5 0.6395701 24.23393
7     6 0.6395534 24.23008
8     7 0.6395513 24.22961
9     8 0.6395511 24.22955
10    9 0.6395511 24.22954

Referensi

Hogg, Robert V., Joseph W. McKean, and Allen T,Craig.2013. Introduction to Mathematical Statistics. 7th ed. Boston: Pearson.

https://rpubs.com/alfanugraha/algoritma-em

https://gerrydito.github.io/posts/2023-9-01-Newton-Raphson-dan-Fisher-Scoring-untuk-MLE/

https://www.statisticshowto.com/em-algorithm-expectation-maximization/