Pemrograman Statistika (Fungsi, OOP dan Optimasi Numerik) - Tugas Responsi 3

Pokok bahasan dalam praktikum STA561 pertemuan 6 dan 7 adalah sebagai berikut:

  1. Praktikum 6: Pemrograman Fungsi dan OOP di R
    • Fungsi
    • OOP (Object Oriented Programming) di R
    • Class System S3
    • Class System S4
  2. Praktikum 7: Optimasi secara Numerik
    • Diferensial
    • Integral
    • Optimasi secara NUmerik
      • Golden Section search
      • Newton Rapshon
      • Fungsi Optimasi Built-in > Optimasi Satu Variabel (Optimize/Optimise) > Optimasi Lebih dari Satu Variabel (Optim) > Regresi Linier
    • MLE: Maximum Likelihood Estimator
    • Latihan

Praktikum 6: Pemrograman Fungsi dan OOP

Fungsi

Fungsi digunakan untuk membuat sekumpulan mekanisme dasar yang dijalankan secara simultan.

function(argumen list){
  expression
  return (value)
}

Contoh 1

angka_acak1 <- function(n,pw){
  x <- runif(n)
  y <- runif(n)
  z <- (x+y)^pw
  return(z)
}

angka_acak1(10,2)
##  [1] 1.6190381 0.3360103 0.6215977 1.4363796 2.4243391 0.2859981 0.7040056
##  [8] 0.9756431 2.4239974 0.4982917

Contoh 2

angka_acak2 <- function(n,pw){
  x <- runif(n)
  y <- runif(n)
  z <- (x+y)^pw
  return(list(x = x, y = y, z = z))
}

angka_acak2(10,2)
## $x
##  [1] 0.47164979 0.80560026 0.16714127 0.61014590 0.08337948 0.81595572
##  [7] 0.60857190 0.94580033 0.40512597 0.10674867
## 
## $y
##  [1] 0.4444645 0.8200039 0.9170443 0.2374087 0.1866832 0.1818032 0.9577259
##  [8] 0.7121519 0.2608791 0.7532914
## 
## $z
##  [1] 0.83926541 2.64258895 1.17545844 0.71834885 0.07293387 0.99552287
##  [7] 2.45328888 2.74880565 0.44356275 0.73966891

Contoh 3

angka_acak3 <- function(n=10, pw=2){
  x <- runif(n)
  y <- runif(n)
  z <- (x+y)^pw
  return(z)
}

angka_acak3()
##  [1] 1.4345168 0.1563539 1.2025178 1.3182684 0.6968917 0.2911175 1.7571489
##  [8] 1.2241553 1.0113047 1.3851274

Contoh 4

angka_acak4 <- function(){
  x <- runif(n)
  y <- runif(n)
  z <- (x+y)^pw
  return(z)
}

n <- 5; pw <- 3
angka_acak4()
## [1] 1.4284204 0.1105117 3.6780662 1.1528365 2.3088689

Latihan 1

Membuat fungsi Median dari suatu vektor

med <- function(vect){
  n <- length(vect)
  vects <- sort(vect)
  if(n%%2 == 1){
    m <- vects[(n+1)/2]}
  else {
    m <- (vects[n/2] + vects[(n/2) + 1])/2}
  return(m)
}

x1 <- c(1,5,3,7,3,4,2,7)
med(x1)
## [1] 3.5

Latihan 2

Membuat fungsi Modus dari suatu vektor

modus <- function(vect){
  v <- unique(vect)
  f <- NULL
  for(i in v)
  { banyak <- sum(vect == i)
    f <- c(f, banyak)
  }
  fmax <- max(f)
  vf <- cbind(v,f)
  mode <- vf[f == fmax,]
  return(mode)
}
 
modus(x1) 
##      v f
## [1,] 3 2
## [2,] 7 2

Latihan 3

Membuat fungsi penduga parameter pada regresi berganda

# Membuat fungsi
p.est <- function(A) {
  if(!is.matrix(A))
    stop("input must be on matrix")
  x1 <- A[,-1]
  y  <- A[,1]
  one <- rep(1,nrow(A))
  x <- cbind(one,x1)
  colnames(x) <- paste("x",1:ncol(x),sep="")
  b.est <- as.vector(solve(t(x) %*% x) %*% (t(x) %*% y))
  names(b.est) <- paste("b",0:(length(b.est)-1),sep="")
  fitted.value <- as.vector(x%*%b.est)
  error <- as.vector(y-fitted.value)
  names(fitted.value) <- names(error) <- 1:nrow(A)
  list(beta.est=b.est,fit.val=fitted.value,error=error)
}

# Memasukkan data
Pendapatan <- c(3.5,3.2,3.0,2.9,4.0,2.5,2.3)
Biaya.Iklan <- c(3.1,3.4,3.0,3.2,3.9,2.8,2.2)
Jumlah.Warung <- c(30,25,20,30,40,25,30)
X <- cbind(Pendapatan,Biaya.Iklan,Jumlah.Warung)
p.est(X)
## $beta.est
##          b0          b1          b2 
## -0.21381852  0.89843390  0.01745279 
## 
## $fit.val
##        1        2        3        4        5        6        7 
## 3.094910 3.277176 2.830539 3.184754 3.988185 2.738116 2.286320 
## 
## $error
##           1           2           3           4           5           6 
##  0.40508982 -0.07717642  0.16946108 -0.28475357  0.01181483 -0.23811608 
##           7 
##  0.01368033
# Melalukan pengecekan dengan fungsi lm
model <- lm(Pendapatan~Biaya.Iklan+Jumlah.Warung)

# Mengeluarkan koefisien dari model regresi
model$coefficients
##   (Intercept)   Biaya.Iklan Jumlah.Warung 
##   -0.21381852    0.89843390    0.01745279
# Mengeluarkan fitted values
model$fitted.values
##        1        2        3        4        5        6        7 
## 3.094910 3.277176 2.830539 3.184754 3.988185 2.738116 2.286320
# Mengeluarkan residuals
model$residuals
##           1           2           3           4           5           6 
##  0.40508982 -0.07717642  0.16946108 -0.28475357  0.01181483 -0.23811608 
##           7 
##  0.01368033

Latihan 4

  • Membuat fungsi bernama three.M yang digunakan untuk menghitung mean, median, dan modus dari suatu vektor modus (tanpa mengunakan fungsi mean, quantile maupun fungsi instan yang telah tersedia di R).

  • Hitung mead, median, modus dari X di bawah ini dengan menggunakan fungsi three.M tersebut

# Membuat fungsi
three.M <- function(vect){
  
  # menghitung banyak data
  n <- length(vect)
  
  # menghitung rerata
  jumlah <- sum(vect)
  rerata <- jumlah/n
  
  # menghitung median
  vects <- sort(vect)     # mengurutkan data
  if(n%%2 == 1){me <- vects[(n+1)/2]}
  else {me <- (vects[n/2] + vects[(n/2) + 1])/2}
  
  # menghitung modus
  v <- unique(vect)
  f <- NULL
  for(i in v){
    banyak <- sum(vect == i)
    f <- c(f, banyak)
  }
  fmax <- max(f)
  vf <- cbind(v,f)
  mode <- vf[f == fmax,]
  
  # output
  list3 <- list("Mean:" = rerata, 
                "Median:" = me,
                "Modus:" = mode)
  return(list3)
}

# Menghitung Mean, Median, Modus suatu data berdistribusi Binomial
set.seed(123)
x2 <- rbinom(100, 10, 0.5)
three.M(x2)
## $`Mean:`
## [1] 4.99
## 
## $`Median:`
## [1] 5
## 
## $`Modus:`
##      v  f
## [1,] 6 24
## [2,] 5 24

Object Oriented Programming / Pemrograman Berbasis Objek

  • PBO menitikberatkan pada identifikasi objek-objek yang terlibat dalam sebuah program dan bagaimana objek tersebut saling berinteraksi.
  • Pada PBO, program terdiri dari beberapa objek.
  • Prinsip PBO yaitu: a) Abtraksi, b) Enkapsulasi, c) Inheritance, dan d) Polymorphism.
  • Pengembangan PBO di R menggunakan Class System S3 dan Class System S4.

Object - S3

Suatu _class_ dalam system S3 tidak didefinisikan dengan ketat. Fungsi class digunakan untuk menjadikan sebuah objek menjadi class yang diinginkan.

Contoh 1

A1 <- c(1:10)
A2 <- matrix(A1,2,5)
class(A1); class(A2)
## [1] "integer"
## [1] "matrix" "array"
A3 <- 1:12
A4 <- letters[1:12]
B1 <- data.frame(A3,A4)
B1
##    A3 A4
## 1   1  a
## 2   2  b
## 3   3  c
## 4   4  d
## 5   5  e
## 6   6  f
## 7   7  g
## 8   8  h
## 9   9  i
## 10 10  j
## 11 11  k
## 12 12  l
# menampilkan class data fram B1
class(B1)
## [1] "data.frame"
# mengambil kolom A4 dari data frame B1
B1$A4
##  [1] "a" "b" "c" "d" "e" "f" "g" "h" "i" "j" "k" "l"
A5 <- 10 + A3 + rnorm(12)

# membuat model linear
B2 <- lm(A5~A3)
B2
## 
## Call:
## lm(formula = A5 ~ A3)
## 
## Coefficients:
## (Intercept)           A3  
##     10.3946       0.9661
# melihat class B2
class(B2)
## [1] "lm"
# melihat daftar methos yang terdapat pada class B2 (lm: linear modelling)
methods(class = class(B2))
##  [1] add1           alias          anova          case.names     coerce        
##  [6] confint        cooks.distance deviance       dfbeta         dfbetas       
## [11] drop1          dummy.coef     effects        extractAIC     family        
## [16] formula        hatvalues      influence      initialize     kappa         
## [21] labels         logLik         model.frame    model.matrix   nobs          
## [26] plot           predict        print          proj           qr            
## [31] residuals      rstandard      rstudent       show           simulate      
## [36] slotsFromS3    summary        variable.names vcov          
## see '?methods' for accessing help and source code
# ringkasa model linier B2
summary(B2)
## 
## Call:
## lm(formula = A5 ~ A3)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.70634 -0.37935 -0.03672  0.38335  1.32502 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 10.39462    0.51345   20.25 1.91e-09 ***
## A3           0.96614    0.06976   13.85 7.51e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8343 on 10 degrees of freedom
## Multiple R-squared:  0.9504, Adjusted R-squared:  0.9455 
## F-statistic: 191.8 on 1 and 10 DF,  p-value: 7.514e-08
# melihat daftar kolom/output dari linier model
names(B2)
##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "xlevels"       "call"          "terms"         "model"
# mengambil "koefesien" lm B2
B2[1]
## $coefficients
## (Intercept)          A3 
##  10.3946172   0.9661381

Mengubah menjadi Class Baru

1) Pengubahan class secara Langsung
Ilustrasi 1:
Membuat object dengan nama Mobil1

# Membuat object Mobil1 dari class list 
Mobil1 <- list(Nama = "Toyot",
              Panjang = 3.5,
              Lebar = 2,
              Kecepatan = 180)

# melihat class object Mobil1
class(Mobil1)
## [1] "list"
# mengubah class object Mobil1 -> class "mobil"
class(Mobil1) <- "mobil"
class(Mobil1)
## [1] "mobil"

Membuat object dengan nama Mobil2

# Membuat object Mobil2 dari class list 
Mobil2 <- list(Nama = "Suzuki",
              Panjang = 1,
              Lebar = 1.8,
              Kecepatan = 150)

# mengubah class Mobil2 dari list -> mobil
class(Mobil2) <- "mobil"
class(Mobil2)
## [1] "mobil"

2) Pengubahan class menggunakan Fungsi Konstruktor
Fungsi konstruktor berfungsi untuk melakukan screening sebelum menambahkan class

Mobil <- function(Nama,Panjang,Lebar,Kecepatan){
  # screening ukuran
  if(Panjang<2 || Lebar<1.5 || Kecepatan<80)
    
    # program behenti jika tidak memenuhi kondisi di atas
    stop("atribut tidak sesuai")
  
  # assign ke dalam object "Mobil"
  Mobil <- list(Nama=Nama, Panjang=Panjang,
                Lebar=Lebar, Kecepatan=Kecepatan)
  
  # pengubahan class list -> "mobil"
  class(Mobil) <- "mobil"
  Mobil
}

Menambahkan data object menggunakan class “Mobil”

# Contoh data yang sesuai dengan kondisi screening
Mobil3 <- Mobil(Nama="Daihatsu",
                Panjang=3,
                Lebar=1.9,
                Kecepatan=120)
Mobil3
## $Nama
## [1] "Daihatsu"
## 
## $Panjang
## [1] 3
## 
## $Lebar
## [1] 1.9
## 
## $Kecepatan
## [1] 120
## 
## attr(,"class")
## [1] "mobil"
# Contoh data yang sesuai dengan kondisi screening
Mobil4 <- Mobil(Nama="Toyota",
                Panjang=1.8,
                Lebar=1.9,
                Kecepatan=120)

Fungsi Aksesor

Fungsi aksesor digunakan untuk mengakses data dalam class. Pengaksesan data dalam class secara langsung dengan menggunakan nama_object$data_yg_diakses tidak direkomendasikan.

1) Mengakses secara Langsung

Mobil2$Nama; Mobil3$Panjang
## [1] "Suzuki"
## [1] 3

2) Mengakses menggunakan Fungsi Aksesor 

# Membuat fungsi aksesor Nama & Kecepatan
nama <- function(objek)objek$Nama
panjang <- function(objek)objek$Panjang

# Mengakses data nama & kecepatan
nama(Mobil2); panjang(Mobil3)
## [1] "Suzuki"
## [1] 3

Fungsi Generik

Contoh pada method print
Dalam ilustrasi class mobil di atas, method print bawaan R menghasilkan output 4 data seperti di bawah ini.

print(Mobil1)
## $Nama
## [1] "Toyot"
## 
## $Panjang
## [1] 3.5
## 
## $Lebar
## [1] 2
## 
## $Kecepatan
## [1] 180
## 
## attr(,"class")
## [1] "mobil"

Method print akan disesuaikan kebutuhan yakni mengeluarkan output hanya berupa data nama dan kecepatan, sehingga perlu membuat fungsi print untuk class mobil

print.mobil <- function(objek){
  print(cat("Nama : ", nama(objek), "\n",
            "Panjang : ", panjang(objek),
            sep=" "))
  }
Mobil1
## Nama :  Toyot 
##  Panjang :  3.5NULL

Object - S4

  • Sistem objek S4 mengatasi kekurangan dalam sistem objek S3
  • Mendefenisikan class baru
setClass("class.name",
  representation(x="type"),
  prototype(x="..."))
  • Mendefinisikan objek baru dari class tertentu
new(class.name,…)

Membuat Object S4

1) Tanpa Fungsi Konstruktor

# membuat class dengan nama "car"
setClass("car", representation(Nama="character",
                               Panjang="numeric",
                               Lebar="numeric",
                               Kecepatan="numeric"))

# membuat objek bernama "Car1" dari class "car"
Car1 <- new("car", Nama="Toyota",
            Panjang=3.5, 
            Lebar=2,
            Kecepatan=180)
Car1
## An object of class "car"
## Slot "Nama":
## [1] "Toyota"
## 
## Slot "Panjang":
## [1] 3.5
## 
## Slot "Lebar":
## [1] 2
## 
## Slot "Kecepatan":
## [1] 180

2) Dengan Menggunakan Fungsi Konstruktor

Car <- function(Nama,Panjang,Lebar,Kecepatan){
  if(Panjang<2 || Lebar<1.5 || Kecepatan<80)
    stop("atribut tidak sesuai")
  new("car", Nama=Nama, 
      Panjang=Panjang,
      Lebar=Lebar, 
       Kecepatan=Kecepatan)
}
Car2 <- Car("Suzuki", 2.4, 1.8, 150)
Car2
## An object of class "car"
## Slot "Nama":
## [1] "Suzuki"
## 
## Slot "Panjang":
## [1] 2.4
## 
## Slot "Lebar":
## [1] 1.8
## 
## Slot "Kecepatan":
## [1] 150
class(Car2); class(Mobil1)
## [1] "car"
## attr(,"package")
## [1] ".GlobalEnv"
## [1] "mobil"

Akses terhadap Slot

1) Akses secara Langsung
Akses terhadap slot class secara langsung secara umum dengan syntax object_s4@data_yg_diakses

# mengakses data Nama dari object Car1
Car1@Nama
## [1] "Toyota"
# mengakses data Kecepatan dari object Car2
Car2@Kecepatan
## [1] 150

2) Akses menggunakan Fungsi Aksesor

# Fungsi aksesor
nama1 <- function(objek)objek@Nama
kecepatan1 <- function(objek)objek@Kecepatan
# mengakses data nama dari object Car1
nama1(Car1)
## [1] "Toyota"
# mengakses data kecepatan dari object Car2
kecepatan1(Car2)
## [1] 150

Class Method

Secara umum syntax dalam membuat class method adalah sebagai berikut:

setMethod("method", "class.name", function(...){...})

Contoh

setMethod(show, "car", function(object) {
  print(cat("Nama : ", nama1(object), "\n",
            "Kecepatan : ", kecepatan1(object),
            sep=" ")
        )}
  )
Car1
## Nama :  Toyota 
##  Kecepatan :  180NULL

Mencipatkan Fungsi Generik S4

Secara umum syntax untuk membuat fungsi generik baru adalah sebagai berikut.

setGeneric("fungsibaru",
           function(objek)
             standarGeneric("fungsibaru"))

Praktikum 7: Optimasi secara Numerik

Latihan 1: Fungsi Integral

Soal 1

  1. Carilah hasil integral tak tentu berikut menggunakan fungsi yac_str dari packages Ryacas \(\int x^2 + 4 dx\)
# install.packages("Ryacas")
library(Ryacas)
## 
## Attaching package: 'Ryacas'
## The following object is masked from 'package:stats':
## 
##     integrate
## The following objects are masked from 'package:base':
## 
##     %*%, diag, diag<-, lower.tri, upper.tri
yac_str("Integrate(x) x^2 + 4")
## [1] "x^3/3+4*x"
  1. Carilah hasil integral tentu berikut menggunakan fungsi integrate \(\int_{-10}^{10} x^2 + 4 dx\)
f1 <- function(x) x^2 + 4
integrate(f1, lower=-10, upper=10)
## 746.6667 with absolute error < 8.3e-12


Soal 2

  1. Carilah hasil integral berikut \(\int t^4 e^{-t} dt\)
yac_str("Integrate(t) t^4 * Exp(-t)")
## [1] "4*(3*((-2)*(t+1)*Exp(-t)-t^2*Exp(-t))-t^3*Exp(-t))-t^4*Exp(-t)"
  1. Carilah hasil integral berikut \(\int_{0}^{\infty} t^4 e^{-t} dt\)
f2 <- function(t) t^(4) * exp(-t)
integrate(f2, lower = 0, upper = Inf)
## 24 with absolute error < 2.2e-05
gamma(5)
## [1] 24

Latihan 2: Fungsi Diferensial

Cari diferensial dari fungsi \(e(x^2)\) menggunakan fungsi D atau deriv

xfs <- expression(exp(x^2))
D(xfs, "x")
## exp(x^2) * (2 * x)
xturunan <- deriv(~x^2, "x")
x <- 2
eval(xturunan)
## [1] 4
## attr(,"gradient")
##      x
## [1,] 4

Optimasi Numerik

Optimasi merupakan proses untuk mencari kondisi yang optimum atau dalam arti paling menguntungkan, sehingga dapat berupa maksimasi atau minimasi.

Beberapa metode statistik menggunakan metode pendugaan nilai optimum dari suatu fungsi tujuan, seperti Metode Kemungkinan Maksimum (likelihood) dan Metode Kuadrat Terkecil.

Metode untuk optimasi numerik,diantaranya:
* Golden Section Search
* Newton-Raphson
* Nelder-Mead

A. Golden Section Search (GSS)

Algoritma Golden Section digunakan untuk menyelesaikan NLP (Non-Linier Programming). GSS menggunakan prinsip “mengurangi daerah batas x yang mungkin menghasilkan nilai fungsi obyektif optimum (maksimum atau minimum) secara iteratif (berulang)”.
* Untuk mendapatkan sebuah titik baru yang simetris, dibutuhkan nilai r (Golden Ratio).
* Batas awal GSS sudah ditentukan.
* GSS merupakan teknik optimasi satu variabel.


Tahapan Optimasi
1. Menentukan range/interval \([a,b]\), \(a\): batas bawah, \(b\): batas atas.
2. Menentukan nilai fungsi di titik \(a\) atau \(f(a)\), dan nilai fungsi di tiitk \(b\) atau \(f(b)\).
3. Menentukan batas baru \(x_{1}\) yang berada dalam \([a,b]\), dimana \(x_{1} = \frac{b-(b-a)*2)}{\sqrt{5}+1}\).
4. Menentukan batas baru \(x_{2}\) yang berada dalam \([a,b]\), dimana \(x_{2} = \frac{a+(b-a)*2)}{\sqrt{5}+1}\).
5. Menentukan nilai fungsi dari \(f(x_{1})\) dan \(f(x_{2})\).
6. Pilihan:
Maksimasi
* Jika \(f(x_{1})>f(x_{2})\), titik optimal berada pada \([a,x_{2}]\), sehingga diperoleh batas baru \(a=a, b=x_{2}\).
* Jika \(f(x_{1})<f(x_{2})\), titik optimal berada pada \([x_{1},b]\), sehingga diperoleh batas baru \(a=x_{1}, b=b\).
Minimasi
* Jika \(f(x_{1})<f(x_{2})\), titik optimal berada pada \([a,x_{2}]\), sehingga diperoleh batas baru \(a=a, b=x_{2}\).
* Jika \(f(x_{1})>f(x_{2})\), titik optimal berada pada \([x_{1},b]\), sehingga diperoleh batas baru \(a=x_{1}, b=b\).
7. Melakukan iterasi sampai diperoleh \(\left | a-b \right | < toleransi\)


Contoh 1 Menghitung nilai maksimum fungsi \(f(x)=-x^2+6x+10\)

x.max <- function(f,a,b,tol=0.0000000001){
  ratio <- 2/(sqrt(5)+1)
  x1 <- b - ratio * (b-a)
  x2 <- a + ratio * (b-a)
  f1 <- f(x1)
  f2 <- f(x2)
  while(abs(b-a) > tol) {
    if (f1 > f2){
      b <- x2
      x2 <- x1
      f2 <- f1
      x1 <- b - ratio * (b-a)
      f1 <- f(x1)
    } else{
      a <- x1
      x1 <- x2
      f1 <- f2
      x2 <- a + ratio * (b-a)
      f2 <- f(x2)
    }
  } 
  return((a+b)/2)
}

f <- function(x){
  -1*x^2+6*x+10
}
# nilai maks jika a=1, b=2
x.max(f,1,2)
## [1] 2
# nilai maks jika a=3, b=5
x.max(f,3,5)
## [1] 3
x <- seq(-10, 10, by = 0.001)
f_x <- -1*x^2+6*x+10

plot(x,f_x, type = "l",
     main = "Grafik Fungsi 2x^2+6x+10",
     xlab = "nilai x",
     ylab = "nilai y",
     xlim = c(-10,10),
     ylim = c(0,25))
abline(v=2, lty=3, lwd=1, col="red")
abline(v=3, lty=3, lwd=1, col="red")
abline(v=7, lty=3, lwd=1, col="red") # nilai maks jika a=6, b=10

B. Newton Rapshon

Metode ini digunakan untuk mencari nilai optimum dari titik awal \(x_{0}\) dalam kasus univariat yang \(f(x)\) dioptimalkan. Tujuannya untuk konvergen ke titik stasioner yang turunannya = 0. Dalam kasus univariat, metode ini menghasilkan urutan iterasi dalam bentuk \[x_{n} = x_{n-1} - \frac{f'(x_{n-1})}{f''(x_{n-1})}\]

Iterasi pada metode ini terus berjalan sampai \(f'(x_{n-1})\approx 0\) atau lebih kecil dari nilai toleransi

Proses Iterasi

newtonr <- function(fx, x0=1){
  fx1 <- deriv(fx, "x")        #turunan pertama
  fx2 <- deriv(D(fx,"x"), "x") #turunan kedua
    e <- 1000
   while (e > 1e-6) {
      x <- x0
     f1 <- attr(eval(fx1), "gradient")[1]
     f2 <- attr(eval(fx2), "gradient")[1]
      e <- abs(f1) #bisa juga e <- abs(x1-x0)
     x1 <- x0-f1/f2
     x0 <- x1
   }
    return(x1)
}

Contoh 1 Menghitung nilai minimum dari fungsi \(f(x) = 6x^{2}-x-12\)

fx1 <- expression(6*x^2-x-12)
newtonr(fx1)
## [1] 0.08333333
# plot fungsi 6*x^2-x-12
plot(x = seq(-10,10, by=0.001), y = 6*x^2-x-12,
     type = "l",
     xlab = "nilai x",
     ylab = "nilai f(x)",
     xlim = c(-10,10),
     ylim = c(-15,50))
abline(v = 0.083333, lty =3, lwd = 1, col = "blue")


Contoh 2 Menghitung nilai minimum dari fungsi \(f(x) = e^{-x} + x^4\)

fx2 <- expression(exp(-x) + x^4)
newtonr(fx2)
## [1] 0.5282519
curve(exp(-x) + x^4, ylab = "nilai f(x)")
abline(v=0.5282519, lty=3, lwd=1, col="brown")


Contoh 3 Menghitung nilai minimum dari \(f(x) = 2x^{2} - x\)

# hitung f(x)
fx3 <- expression(x^2 - x)

# nilai minimum
newtonr(fx3)
## [1] 0.5
# plot fungsi
curve(x^2 - x, ylab = "fx3")
abline(v = newtonr(fx3), lty = 3, lwd = 1, col = "brown")

C. Fungsi Optimasi Built-in

Pada R algoritma Nelder-Mead (salah satu metode optimasi untuk fungsi lebih dari satu variabel) diterapkan pada fungsi:
* optimize() atau optimize() untuk fungsi satu variabel
* optim() untuk fungsi lebih dari satu variabel


Optimasi 1 Variabel: optimize()

Contoh 1
Menghitung nilai minimun dari \(f(x) = \left (x-\frac{1}{3} \right )^{2}\)

f <- function(x,a){
  (x-a)^2
}

# titik minimum (x,y) dg fungsi optimize(nama_fungsi, vektor_interval)
titikMin <- optimize(f, c(0,1), tol=0.0001, a=1/3)
titikMin
## $minimum
## [1] 0.3333333
## 
## $objective
## [1] 0


Contoh 2 Menghitung titik maksimum dari fungsi \(f(x) = -x^2+4x+5\)

Maksimum pada interval [-1,5]

f <- function(x){-x^2+4*x+5}

# titik maksimum (x,y) dg fungsi optimize(nama_fungsi, vektor_interval)


titikMaks <- optimise(f, c(-1,1), maximum = TRUE)
titikMaks
## $maximum
## [1] 0.999959
## 
## $objective
## [1] 7.999918
curve(f, from=-1, to=5)
abline(v=2, h=9, lty=3, lwd=1, col="brown")


Maksimum pada interval [0,3]

titikMaks <- optimise(f, c(0,3), maximum = TRUE)
titikMaks
## $maximum
## [1] 2
## 
## $objective
## [1] 9
curve(f, from=0, to=3)
abline(v=2, h=9, lty=3, lwd=1, col="brown")


Contoh 3 Menghitung titik maksimum dan minimum dari fungsi \(f(x) = sin(x) + sin(2x) + cos(3x)\)

f3 <- function(x) sin(x) + sin(2*x) + cos(3*x)

# plot fungsi f3
curve(f3, from=0, to=2*pi)

Minimum Lokal

optimize(f3, interval = c(0, 2*pi))
## $minimum
## [1] 3.033129
## 
## $objective
## [1] -1.054505

Minimum Global

optimize(f3, interval = c(4, 2*pi))
## $minimum
## [1] 5.273383
## 
## $objective
## [1] -2.741405

Maksimum Lokal

optimize(f3, interval = c(0,2*pi), maximum = TRUE)
## $maximum
## [1] 4.0598
## 
## $objective
## [1] 1.096473

Maksimum Global

optimize(f3, interval = c(0,1.5), maximum = TRUE)
## $maximum
## [1] 0.3323289
## 
## $objective
## [1] 1.485871


Contoh 4 ……..

Optimasi Lebih Dari 1 Variabel: optim()

Contoh 1

Carilah nilai \(x_{1}\) dan \(x_{2}\) agar fungsi berikut \(f(xf(x_{1}, x_{2}) = 100 \left ( x_{2} - x_{1}^{2} \right )^{2} + \left ( 1-x_{1} \right )^{2}_{1}, x_{2}) = 100 \left ( x_{2} - x_{1}^{2} \right )^{2} + \left ( 1-x_{1} \right )^{2}\) minimum!

f4 <- function(x){
  x1 <- x[1]
  x2 <- x[2]
  100 * (x2 - (x1)^2)^2 + (1-x1)^2
}

# mencari x1 dan x2 fungsi
optim(c(-1.2,1),f4)
## $par
## [1] 1.000260 1.000506
## 
## $value
## [1] 8.825241e-08
## 
## $counts
## function gradient 
##      195       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

Regresi Linier

Contoh Kasus

Lakukan pendugaan parameter regresi dengan memiminimumkan Jumlah Kuadrat Galat (residual sum of square) dari data berikut. Kemudian bandingkan hasilnya dengan output dari fungsi lm!

# data frame
data5 <- data.frame(x = c(1,2,3,4,5,6),
                    y = c(1,3,5,6,8,12))
data5
##   x  y
## 1 1  1
## 2 2  3
## 3 3  5
## 4 4  6
## 5 5  8
## 6 6 12
JKG <- function(data,b){
  with(data, sum((b[1] + b[2]*x - y)^2))
}

hasil1 <- optim(par = c(1,1), fn = JKG,  data = data5)
hasil1
## $par
## [1] -1.266302  2.028449
## 
## $value
## [1] 2.819048
## 
## $counts
## function gradient 
##       67       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL
hasil2 <- lm(y = x, data = data5)
## Warning in if (ret.y) z$y <- y: the condition has length > 1 and only the first
## element will be used
hasil2
## 
## Call:
## lm(data = data5, y = x)
## 
## Coefficients:
## (Intercept)            y  
##      0.7327       0.4744
# Plot regresi
plot(data5)
abline(hasil1$par, col=4)

hasil1$par
## [1] -1.266302  2.028449
hasil2$coefficients
## (Intercept)           y 
##   0.7327394   0.4743875
hasil1$value
## [1] 2.819048
# Jumlah Kuadrat Galat
sum(hasil2$residuals^2)
## [1] 0.6592428


D. MLE: Maximum Likelihood Estimator

Maximum Likelihood Estimator merupakan metode yang paling sering digunakan untuk menduga parameter sebaran.

# install.packages("bbmle")
library(bbmle)
## Loading required package: stats4
Contoh 1: Sebaran Normal

Jika \(x_{1}, x_{2}, ..., x_{n}\) berasal dari peubah acak \(X\sim N(\mu,\sigma)\) berikut penduga \(\mu\) dan \(\sigma\) menggunakan MLE.

  • menggunakan fungsi optim()
lnorm <- function (para ,xd){
  nilai <- -1*sum(dnorm(xd, mean = para[1], sd = para[2], log = TRUE))
  return (nilai)
  }
set.seed(2)
x <- rnorm(10,2,5)
hasil <- optim(c(1,1), lnorm, xd = x)
hasil
## $par
## [1] 3.056350 4.671425
## 
## $value
## [1] 29.6057
## 
## $counts
## function gradient 
##       63       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL
  • menggunakan fungsi mle()
lnorm <- function (mean, sd){
  nilai <- -1*sum(dnorm(x, mean = mean, sd = sd, log = TRUE))
  return (nilai)
  }
suppressWarnings(
  estnorm <- mle2(minuslogl = lnorm, start=list(mean = 1, sd = 2))
  )
summary(estnorm)
## Maximum likelihood estimation
## 
## Call:
## mle2(minuslogl = lnorm, start = list(mean = 1, sd = 2))
## 
## Coefficients:
##      Estimate Std. Error z value     Pr(z)    
## mean   3.0557     1.4775  2.0682   0.03862 *  
## sd     4.6722     1.0447  4.4721 7.744e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## -2 log L: 59.2114

Latihan

Pendugaan Regresi Linier Menggunakan MLE

Model regresi linier sederhana:
\(Y_{i} = B_{0} + B_{1}X_{i} + \varepsilon_{i}\) dimana \(i = 1,2,...n\)

1. Membangkitkan Bilangan Acak

# Membangkitkan bilangan acak X ~ N(mean=12, sd=4)
data.x <- rnorm(n = 50, mean = 12, sd = 4)
# Membangkitkan nilai Y dengan menentukan nilai B0 adalah 8 dan B1 adalah 3
B0 <- 8
B1 <- 3

# Menentukan error dari data dengan sd=1
error.sd <- 1
error <- rnorm(n = 50, mean = 0, sd = 1)
data.y <- B0 + B1 * data.x + error

# Plot 
plot(data.x, data.y)

Dari plot di ata terlihat bahwa data X dan Y yang dibangkitkan memiliki hubungan linier.


2. Pendugaan Regresi Linier dengan MLE Diasumsikan mengikuti sebaran Normal \(\mu = B_{0} + B_{1}X_{1}; \sigma=sd\).
Akan dicari nilai penduga \(B_{0}\) dan \(B_{1}\)

# MLE

lm.loss <- function(parameter){
  
  B0.parameter <- parameter[1]  # parameter B0/intercept
  B1.parameter <- parameter[2]  # parameter B1
  error.sigma <- parameter[3]   # standar deviasi
  
  # prosedur optimasi: dibuat nilai toleranis jika standar deviasi tidak valid 
  # misal negatif, maka dikembalikan kepada deviance yang sangat tinggi
  if(error.sigma < 0){ deviance <- 10000000}
  
  if(error.sigma > 0){
    likelihoods <- dnorm(data.y, 
                         mean = B0.parameter + B1.parameter*data.x, 
                         sd = error.sigma)
    
    log.likelihoods <- log(likelihoods)
    
    deviance <- -2 * sum(log.likelihoods)
  }
  return(deviance)
}


parameter.fits <- optim(par = c(1,1,10),
                        fn = lm.loss,
                        hessian = T)

parameter.fits
## $par
## [1] 8.187312 2.976223 1.135375
## 
## $value
## [1] 154.6134
## 
## $counts
## function gradient 
##      144       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL
## 
## $hessian
##              [,1]         [,2]         [,3]
## [1,]  77.57495533   953.843970   0.01977574
## [2,] 953.84397028 13472.658955   0.23635901
## [3,]   0.01977574     0.236359 155.25934269

Nilai $parameter pada output di atas menyatakan
\(B_{0} = 7.847328\) \(B_{1} = 3.020590\) \(sd = 1.075546\)


3. Membandingkan dengan FUngsi lm() Membandingkan hasil optimasi MLE poin 2 dengan fungsi regresi linier menggunakan fungsi lm().

regresi <- lm(data.y ~ data.x)
regresi
## 
## Call:
## lm(formula = data.y ~ data.x)
## 
## Coefficients:
## (Intercept)       data.x  
##       8.187        2.976

Diperoleh bahwa nilai penduga \(B_{0}\) dan \(B_{1}\) yang dihasilkan dengan optimasi MLE tidak jauh berbeda (hampir sama) dengan yang dihasilkan fungsi lm().

Referensi

Dito, G.A. Praktikum Komputasi Statistika. Diakses dari https://bookdown.org/gerryalfadito/praktikum-komputasi-statistika/

Raharjo, Mulianto. (2021). STA 561 Pemrograman Statistika: Pemrograman Fungsi & OOP dan Optimasi . Diakses dari https://newlms.ipb.ac.id/

Soleh. A.M. (2021). STA 561 Pemrograman Statistika: Pemrograman Fungsi R dan Optimasi secar Numerik. Diakses dari https://newlms.ipb.ac.id/