Tujuan

Analisis ini bertujuan untuk menyelesaikan/ mencari koefisien beta bagi permasalahan data trees dengan menggunakan metode iterasi Jacobi.

Library

library(corrplot)
library(cmna)
library(Rlinsolve)

Persiapan Data

head(trees)

##   Girth Height Volume
## 1   8.3     70   10.3
## 2   8.6     65   10.3
## 3   8.8     63   10.2
## 4  10.5     72   16.4
## 5  10.7     81   18.8
## 6  10.8     83   19.7

Eksplorasi Data

summary(trees)

##      Girth           Height       Volume     
##  Min.   : 8.30   Min.   :63   Min.   :10.20  
##  1st Qu.:11.05   1st Qu.:72   1st Qu.:19.40  
##  Median :12.90   Median :76   Median :24.20  
##  Mean   :13.25   Mean   :76   Mean   :30.17  
##  3rd Qu.:15.25   3rd Qu.:80   3rd Qu.:37.30  
##  Max.   :20.60   Max.   :87   Max.   :77.00

boxplot(trees)

korelasi <- cor(trees)
corrplot(korelasi, method = "number",type = "upper")

Matriks

Langkah selanjutnya adalah bentuk dalam sebuah matriks

pred <- cbind(intercept=1, Girth=trees$Girth, Height=trees$Height)
head(pred)

##      intercept Girth Height
## [1,]         1   8.3     70
## [2,]         1   8.6     65
## [3,]         1   8.8     63
## [4,]         1  10.5     72
## [5,]         1  10.7     81
## [6,]         1  10.8     83

Kemudian, membentuk vektor b

resp<- trees$Volume
head(resp)

## [1] 10.3 10.3 10.2 16.4 18.8 19.7

Transformasi

A <- t(pred) %*% pred #(X'X)
b <- t(pred) %*% resp #(X'y)
Ab <- cbind(A,b)

Metode Iteratif Jacobi

Dengan menggunakan metode iterative Jacobi ?

x0 <- c(0,0,0)
result <- lsolve.jacobi(A,b,x0,maxiter=10000)

x <- result$x
iter <- result$iter
cat("The solution is x =", x, " after ", iter, " iterations\n")

## The solution is x = -57.89436 4.709369 0.3378165  after  3052  iterations

Metode Eliminasi Gauss/Row Echelon Form

replace_row <- function(m, row1, row2, k){
  m[row2, ] <- m[row2, ] + m[row1, ]*k
  return(m)
}

ref_matrix <- function(a){
  m <- nrow(a)
  n <- ncol(a)
  piv <- 1
  
  # cek elemen diagonal apakah bernilai nol
  for(row_curr in 1:m){
    if(piv <= n){
      i <- row_curr
      while(a[i, piv] == 0 && i < m){
        i <- i+1
        if(i > m){
          i <- row_curr
          piv <- piv+1
          if(piv > n)
            return(a)
        }
      }
      
      # jika diagonal bernilai nol, lakukan row swapping
      if(i != row_curr)
        a <- swap_row(a, i, row_curr)
      
      # proses triangulasi untuk membentuk matriks segitiga atas
      for(j in row_curr:m)
        if(j != row_curr){
          c <- a[j, piv]/a[row_curr, piv]
          a <- replace_row(a, row_curr, j, -c)
        }
      piv <- piv+1
    }
  }
  return(a)
}

Dengan menggunakan metode eliminasi Gauss/ Row echelon form:

ref_matrix(Ab)[]

##           intercept    Girth    Height          
## intercept        31 410.7000 2356.0000  935.3000
## Girth             0 295.4374  311.5000 1496.6435
## Height            0   0.0000  889.5641  301.7857

Menggunakan substitusi balik didapatkan koefficient variable height: 301.7857/889.5641 = 0.3392512

GRHeight <- 301.7857/889.5641
GRHeight

## [1] 0.3392512

GRGirth <- 1496.6435/311.5000

Metode Eliminasi Gauss Jordan

Dengan menggunakan metode eliminasi Gauss Jordan:

nrow <- nrow(A)
Ugmt.mtx <- cbind(A,b)
Ugmt.mtx[1,] <- Ugmt.mtx[1,]/Ugmt.mtx[1,1]
for (i in 2:nrow){
  for (j in i:nrow) {
    Ugmt.mtx[j, ] <- Ugmt.mtx[j, ] - Ugmt.mtx[i-1, ] * Ugmt.mtx[j, i-1]
  }
  Ugmt.mtx[i,] <- Ugmt.mtx[i,]/Ugmt.mtx[i,i]
}
for (i in j:2){
  for (j in i:2-1) {
    Ugmt.mtx[j, ] <- Ugmt.mtx[j, ] - Ugmt.mtx[i, ] * Ugmt.mtx[j, i]
  }
}
Ugmt.mtx

##           intercept Girth Height            
## intercept         1     0      0 -57.9876589
## Girth             0     1      0   4.7081605
## Height            0     0      1   0.3392512

Metode Iteratif Gauss Siedel

gaussSeidel <- function(A, b, epsilon, maxIterations) {
  # Extract the size of the system
  n <- length(b)
  
  # Initialize the solution vector x
  x <- numeric(n)
  
  # Iterate until either epsilon is reached or maxIterations is reached
  for (iteration in 1:maxIterations) {
    x_old <- x
    
    # Update each variable using the Gauss-Seidel formula
    for (i in 1:n) {
      x[i] <- (b[i] - sum(A[i, -i] * x[-i])) / A[i, i]
    }
    
    # Check if epsilon is reached
    if (max(abs(x - x_old)) < epsilon) {
      break
    }
  }
  
  # Return the solution and the number of iterations
  list(x = x, iterations = iteration)
}

result <- gaussSeidel(A, b, 1e-6, 1000)
x <- result$x
iterations <- result$iterations
cat("Number of iterations:", iterations, "\n")

## Number of iterations: 1000

cat("Solution:", x, "\n")

## Solution: -57.77524 4.71234 0.3357443

Metode Base R

Dengan menggunakan fungsi base R

## fungsi solve
baser <- data.frame(solve(A,b))
baser

##           solve.A..b.
## intercept -57.9876589
## Girth       4.7081605
## Height      0.3392512

Perbandingan

hasil <- data.frame("Hasil"=c("Intersep","Girth","Height"),
                    "Iteratif Jacobi"=c(x[1],x[2],x[3]),
                    "Eliminasi Gauss/Row Echelon Form"=c("",GRGirth,GRHeight),
                    "Eliminasi Gauss Jordan"=c(Ugmt.mtx[1,4],Ugmt.mtx[2,4],Ugmt.mtx[3,4]),
                    "Iteratif Gauss Siedel"=c(x[1],x[2],x[3]),
                    "Base R"=c(baser[1,1],baser[2,1],baser[3,1]),check.names = FALSE)
knitr::kable(hasil,"pipe",align = "lccccc")

Apabila dilihat dilihat dari perbandingan metode iteratif Jacobi dengan metode lainnya, hasil yang didapatkan tidak jauh berbeda. Maka dari itu, dapat dikatakan bahwa nilai koefisien beta yang dicari menggunakan metode iterasi Jacobi sudah sesuai dan berhasil dilakukan.