Resume Perkuliahan Komputasi Statistik Minggu 11

Fikri Dwi Alpian - 120450022 - RB

Optimasi Numeruk

Beberapa metode statistik menggunakan metode pendugaan nilai optimum dari suatu fungsi tujuan

Maximum Likelihood Method: Mencari nilai maksimum dari fungsi kemungkinan.
Least Square Method: Mencari nilai minimum dari fungsi jumlah kuadrat galat.

Mendapatkan nilai optimum dari suatu fungsi merupakan suatu teknik optimasi numerik. Beberapa metode yang sudah dikembangkan: - Golden section search - Newton-Raphson - Nelder-Mead - dan lainnya.

Golden Section Search

Mencari nilai minimum untuk fungsi variabel tunggal dari suatu selang yang diketahui.
Dimulai dengan selang [a,b] yang memuat nilai minimum pada suatu fungsi
Perkecil selang [a,b] menjadi selang [a’,b’] yang memuat minimum
Ini merupakan langkah iteratif, dengan pemberhentian iterasi dilakukan ketika |b’-a’| sudah lebih kecil dari batas toleransi.
Pemilihan selang yang baru dengan bantuan bilangan golden ratio 𝜙=(√5+1)/2
Hitung bilangan x1 dan x2:
𝑥1=𝑏−((𝑏−𝑎)/𝜙)
𝑥2=𝑎+((𝑏−𝑎)/𝜙)
Hitung 𝑓(𝑥1) dan 𝑓(𝑥2 )
Jika 𝑓(𝑥1 )>𝑓(𝑥2), maka [𝑎′,𝑏′ ]=[𝑥1,𝑏]
Jika 𝑓(𝑥1 )<𝑓(𝑥2 ), maka [𝑎′,𝑏′ ]=[𝑎,𝑥2]
Iterasi selanjutnya akan dilakukan jika |^′−𝑎′ |>ε
Jika |𝑏′−𝑎′ |<ε, maka nilai minimum berada pada titik (𝑎′+^′)/2

# Golden Section Search
golden<-function(a,b,f,eps){
  while(abs(b-a)>eps){
    rat<-(sqrt(5)+1)/2
    x1<-b-(b-a)/rat
    x2<-a+(b-a)/rat
    if(f(x1)>=f(x2)){
      a=x1
    } else {
      b=x2
    }
  }
  return ((a+b)/2)
}

Untuk menguji coba, masukkan fungsi f<-function(x) x^2
Kemudian jalankan fungsi golden(-4,5,f,10^-6)

f=function(x) x^2
golden(-4,5,f,10^-6)

## [1] 2.03118e-07

golden(3,6,f,10^-7)

## [1] 3

f=function(x) 2*x-x^2-1
golden(0,9,f,10^-6)

## [1] 9

golden<-function(a,b,f,eps){
  while(abs(b-a)>eps){
    rat<-2
    x1<-b-(b-a)/rat
    x2<-a+(b-a)/rat
    if(f(x1)>=f(x2)){
      a=x1
    } else {
      b=x2
    }
  }
  return ((a+b)/2)
}
f=function(x) 2*x-x^2-1
golden(0,9,f,10^-7)

## [1] 9

Untuk kendala memaksimumkan, cukup mengganti tanda dalam menentukan selang selanjutnya
Jika 𝑓(𝑥1 )<𝑓(𝑥2), maka [𝑎′,𝑏′ ]=[𝑥1,𝑏]
Jika 𝑓(𝑥1 )>𝑓(𝑥2 ), maka [𝑎′,𝑏′ ]=[𝑎,𝑥2]
Latihan: Carilah nilai maksimum dari fungsi 𝑓(𝑥)=2𝑥−𝑥^2−1 menggunakan GSS dengan selang awal [0,9]

golden_max<-function(a,b,f,eps){
  while(abs(b-a)>eps){
    rat<-(sqrt(5)+1)/2
    x1<-b-(b-a)/rat
    x2<-a+(b-a)/rat
    if(f(x1)<=f(x2)){
      a=x1
    } else {
      b=x2
    }
  }
  return ((a+b)/2)
}
f=function(x) 2*x-x^2-1
golden_max(0,9,f,10^-7)

## [1] 1

Jika suatu fungsi memiliki turunan pertama dan kedua, maka nilai minimum dapat dicari menggunakan metode Newton-Raphson
𝑥𝑛=𝑥(𝑛−1)−(𝑓′ (𝑥(𝑛−1) ))/(𝑓′′ (𝑥(𝑛−1) ) )
Metode Newton-Raphson lebih cepat dibandingkan Golden Section Search
Ini merupakan metode iteratif, iterasi berhenti jika 𝑓′ (𝑥(𝑛−1) ) dekat dengan 0 atau lebih kecil dari nilai toleransi

f=expression(x^2-x)
D(f,'x')

## 2 * x - 1

g=D(f,'x')
# g(8)  #karena tidak fungsi dasar eval maka turunan tidak dapat ditemukan

f=expression(2*x^2+1)
fx=deriv(f,'x')
x=6
attr(eval(fx),'gradient')[1]

## [1] 24

Untuk mendapatkan derivative suatu fungsi pada nilai tertentu, gunakan fungsi dasar deriv, kemudian menggunakan fungsi dasar eval.

f=expression(4*x^2+2*x)
fx=deriv(f,'x')
x=7
attr(eval(fx),'gradient')[1]

## [1] 58

attr(eval(fx),'gradient')[1]

## [1] 58

newrap=function(f,x0,eps){
  f1x=deriv(f,'x')
  f2x=deriv(D(f,'x'),'x')
  e=eps+1
  while(e>eps){
    x=x0
    f1=attr(eval(f1x),'gradient')[1]
    f2=attr(eval(f2x),'gradient')[1]
    e=abs(f1)
    x0=x0-f1/f2
  }
  return(x0)
}

f=expression(x^2)
newrap(f,3,10^-5)

## [1] 0

Cobalah mencari nilai minimum fungsi 𝑓(𝑥)=4𝑥^2−3𝑥−7

f=expression(4*x^2-3*x-7)
newrap(f,4,10^-5)

## [1] 0.375

Latihan: Carilah nilai minimum dari fungsi berikut menggunakan metode Newton-Raphson - 𝑓(𝑥)=𝑒^{(−𝑥)+𝑥}4

f1=expression(exp(1)^(-x)+x^(4))
newrap(f1,10,10^-5)

## [1] 0.5282519

𝑔(𝑥)=𝑥^2−𝑥

f2=expression(x^2-x)
newrap(f2,7,10^-5)

## [1] 0.5

Fungsi Optimasi Built-in

Optimasi Golden Section dan Newton-Raphson sebelumnya digunakan untuk mendapatkan nilai optimum dari satu peubah.
Algoritma Nelder-Mead adalah salah satu metode optimasi untuk fungsi yang memiliki lebih dari satu peubah
R memiliki fungsi dasar dengan salah satu algoritmanya adalah Nelder-Mead:
optimize atau optimize (satu variabel)
optim (lebih dari satu variabel)
Digunakan untuk mendapat nilai minimum dari suatu fungsi dengan satu peubah.
Misal akan dicari nilai minimum fungsi: 𝑓(𝑥)=𝑥^2−4𝑥+4

f<-function(x) x^2-4*x+4
xmin<-optimize(f,c(6,7),tol=10^-6)
xmin

## $minimum
## [1] 6
## 
## $objective
## [1] 16

Latihan: Nilai minimum fungsi: 𝑓(𝑥)=𝑥^3−2𝑥2+2𝑥−1 pada [-1,4]

f<-function(x) x^3-2*x^2+2*x-1
xmin<-optimize(f,c(-1,4),tol=10^-6)
xmin

## $minimum
## [1] -0.9999994
## 
## $objective
## [1] -5.999994

Untuk mencari nilai minimum dari fungsi lebih dari satu peubah.
Misalkan mencari nilai 𝑥_1 dan 𝑥_2 untuk meminimumkan
𝑓(𝑥1,𝑥2 )=100(𝑥2−𝑥1^2 )^2+(1−𝑥1 )^2

f=function(x){
  100*(x[1]-x[2]^2)^2+(1-x[1])^2
}
f(c(-1,2))

## [1] 2504

optim(c(-1.2,1),f)

## $par
## [1]  0.9993547 -0.9996045
## 
## $value
## [1] 2.537102e-06
## 
## $counts
## function gradient 
##       65       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

f<-function(x1,x2){
  100*(x1-x2^2)^2+(1-x1)^2
}
f(-1,2)

## [1] 2504

Contoh aplikasi optim: Regresi
Meminimumkan residu jumlah kuadrat dari data
Misalkan fungsi yang akan menjadi hasil regresi adalah y=par1 + par2 * x
Buat data sebagai berikut: dat=data.frame(x=c(1,2,3,4,5,6),y=c(1,3,5,6,8,12))

dat=data.frame(x=c(1,2,3,4,5,6), y=c(1,3,5,6,8,12))
LR=function(data,par){
  with(data,sum((par[1]+par[2]*x-y)^2))
}
result=optim(par=c(0,1),fn=LR,data=dat)
result

## $par
## [1] -1.266846  2.028620
## 
## $value
## [1] 2.819048
## 
## $counts
## function gradient 
##       89       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

Latihan: Dengan data yang sama, berikan regresi kuadratiknya

dat=data.frame(x=c(1,2,3,4,5,6), y=c(1,3,5,6,8,12))
LR=function(data,par){
  with(data,sum((par[1]+par[2]*x+par[3]*x^2-y)^2))
}
result=optim(par=c(0,1,0),fn=LR,data=dat)
result

## $par
## [1] 0.4002723 0.7782111 0.1786361
## 
## $value
## [1] 1.628572
## 
## $counts
## function gradient 
##      100       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL