Optimasi merupakan suatu proses untuk mecari kondisi yang optimum, dalam arti paling maksimal, atau paling minimal. Optimasi dapat dilakukan dengan meminimumkan atau memaksimumkan suatu fungsi

Fungsi Kalkulus

Diferensial

Untuk mendapatkan turunan dari suatu fungsi pada pemrograman R dapat menggunakan fungsi

  • D(expr, simbol) : jika hasil turunan merupakan suatu fungsi
  • deriv(~fungsi, simbol) : jika akan memasukkan nilai dari hasil turunan pada suatu fungsi

Contoh, akan dicari turunan dari fungsi berikut

\(\begin{aligned} f(x) = e^{(x^2)} \end{aligned}\)

xfs <- expression(exp(x^2))
D(xfs, "x")
## exp(x^2) * (2 * x)

\(\begin{aligned} f(x) = x^2 \end{aligned}\)

untuk nilai \(x=3\)

xturunan <- deriv(~x^2, "x")
x<-3
eval(xturunan)
## [1] 9
## attr(,"gradient")
##      x
## [1,] 6

Integral

Fungsi integral yang dapat digunakan pada R adalah integrate(fungsi, lower, upper). Selain itu, terdapat fungsi yac_str() pada package Ryacas untuk mencari integral tak tentu.

Contoh menggunakan fungsi integrate:

\(\begin{aligned} \int x^2 dx \end{aligned}\)

fs <- function (x) x^2
integrate(fs,0,1)
## 0.3333333 with absolute error < 3.7e-15

\(\begin{aligned} \int x^2 +4x dx \end{aligned}\)

f1 <- function(x) x^2 + 4*x
integrate(f1, lower=-10, upper=10)
## 666.6667 with absolute error < 7.6e-12

\(\begin{aligned} \int t^4 e^{-t} dt \end{aligned}\)

f2 <- function(t)t^4 * exp(-t)
integrate(f2, lower=0, upper=Inf)
## 24 with absolute error < 2.2e-05

Contoh dengan fungsi yac_str pada paket 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

\(\begin{aligned} \int x^2 + 4xdx \end{aligned}\)

yac_str("Integrate(x) x^2 + 4*x")
## [1] "x^3/3+2*x^2"

\(\begin{aligned} \int t^4 e^{-t} dt \end{aligned}\)

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)"

Optimasi Numerik

Mendapatkan nilai optimum dari suatu fungsi merupakan suatu teknik optimasi numerik. Optimasi dapat dilakukan pada satu variabel atau pun lebih dari satu variabel. Fungsi dengan satu variabel, misalkan \(y=f(x)\), kemudian akan dicari nilai \(x\)yang memaksimumkan/ meminimumkan nilai \(y\). Sehingga, \(x\) yang diperoleh merupakan nilai \(x\) yang paling optimum. Beberapa metode yang dapat dilakukan untuk optimasi satu variabel adalah metode golden section search, metode newton raphson. Fungsi dengan banyak variabel, misalkan \(y=f ( x_1, x_2, …, x_n)\), kemudian akan dicari nilai \(x_i\) yang meminimumkan fungsi \(y\).

Paket R untuk optimisasi

Problem type Package Routine
General purpose (1-dim) built-in optimize(...)
General purpose (n-dim) built-in optim(...)
Linear Programming lpSolve lp(...)
Quadratic Programming quadprog solve.QP(...)
Non-Linear Programming optimize optimize(...)
optimx optimx(...)
General Interface ROI ROI_solve(...)

Newton Raphson

Jika suatu fungsi memiliki turunan pertama dan kedua, maka nilai minimum dapat menggunakan metode Newton Raphson. Kelebihan metode ini adalah hanya memerlukan satu nilai untuk inisial. Kelemahannya adalah kita harus yakin f(x) memiliki turunan pertama dan turunan kedua. Jika di golden section tidak perlu ada turunan pertama dan turunan kedua.

Fungsi Newton Raphson

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

Hitung nilai minimum untuk fungsi- fungsi berikut.

\(\begin{aligned} f(x) = 4x^2 - 3x - 7 \end{aligned}\)

fx <- expression(4*x^2-3*x-7)
newtonr(fx,3)
## [1] 0.375

\(\begin{aligned} f(x) = e^{-x}+x^4 \end{aligned}\)

fx <- expression(exp(-x)+x^4)
newtonr(fx)
## [1] 0.5282519

\(\begin{aligned} f(x) = x^2 - x \end{aligned}\)

fx <- expression(x^2-x)
newtonr(fx)
## [1] 0.5

Fungsi Optimasi Built-in

Algoritma Nelder Mead adalah salah satu metode optimasi untuk fungsi yang memiliki lebih dari satu variabel. Di dalam R, fungsi optimasi dengan salah satu algoritma tersebut adalah

  • optimize atau optimise untuk menduga parameter/ mencari nilai minimum dari satu peubah
  • optim untuk lebih dari satu peubah

Fungsi optimize/optimise

\(\begin{aligned} f(x) = \left( x - \frac{1}{3} \right)^2 \end{aligned}\)

f <- function(x) ((x-(1/3))^2) # membuat fungsi tujuan
curve(f)

xmin <- optimize(f, c(0,1), tol=0.0001) # tolerance optional
xmin
## $minimum
## [1] 0.3333333
## 
## $objective
## [1] 0

Carilah titik maksimum dan minimum dari fungsi berikut:

\(\begin{aligned} f(x) = \mathrm{sin}(x) + \mathrm{sin}(2x) + \mathrm{cos}(3x) \end{aligned}\)

f3 <- function(x) sin(x) + sin (2*x) + cos(3*x)
curve(f3, from =0, to = 2* pi)

optimize(f3, interval = c(0, 2*pi)) #minimum lokal
## $minimum
## [1] 3.033129
## 
## $objective
## [1] -1.054505
optimize(f3, interval = c(4, 2*pi)) #minimum global
## $minimum
## [1] 5.273383
## 
## $objective
## [1] -2.741405
optimize(f3, interval = c(0, 2*pi), maximum = T) #maksimum lokal
## $maximum
## [1] 4.0598
## 
## $objective
## [1] 1.096473
optimize(f3, interval = c(0, 1.5), maximum = T) #maksimum global
## $maximum
## [1] 0.3323289
## 
## $objective
## [1] 1.485871

Fungsi optim

Digunakan untuk mencari nilai minimum dari fungsi yang lebih dari satu peubah. Contoh mencari nilai \(x_1\) dan \(x_2\), yang mebuat \(f(x_1,x_2) = 100(x_2 - x_1^2)^2 + (1-x_1)^2\)

fr <- function (x){ # tetap dituliskan dalam sebuah vektor, akan diduga x
  x1<- x[1]
  x2 <- x[2]
  100 * (x2-x1^2)^2 + (1-x1)^2 # ini adalah nilai fungsi objetivenya
} 

# argumen pertama adalah nilai inisial, karena menduga x vektor berukuran 2 maka dimasukkan nilai inisialnya
optim (c(-1.2,1),fr) 
## $par
## [1] 1.000260 1.000506
## 
## $value
## [1] 8.825241e-08
## 
## $counts
## function gradient 
##      195       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL
optim(par=2, fn=f3)
## $par
## [1] 3.033154
## 
## $value
## [1] -1.054505
## 
## $counts
## function gradient 
##       32       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL
optim(par=4, fn=f3)
## $par
## [1] 5.273389
## 
## $value
## [1] -2.741405
## 
## $counts
## function gradient 
##       32       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

Jika akan mencari nilai maksimum, maka function dikali -1.

f3a <- function(x) -1 * f3(x)
optim(par=4, fn=f3a)
## $par
## [1] 4.059814
## 
## $value
## [1] -1.096473
## 
## $counts
## function gradient 
##       28       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL
optim(par=5.5, fn=f3a)
## $par
## [1] 6.615576
## 
## $value
## [1] -1.485871
## 
## $counts
## function gradient 
##       28       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL
optim(par=1, fn=f3a)
## $par
## [1] 0.3323242
## 
## $value
## [1] -1.485871
## 
## $counts
## function gradient 
##       32       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

Fungsi polinomial

\(\begin{aligned} f(x) = 4x^4 - 2x^3 - 3x \end{aligned}\)

f4 <- function(x) 4*x^4-2*x^3-3*x
curve(f4, from = -1, to = 1.5)

optim(par = c(-0.5), fn = f4)
## $par
## [1] 0.728418
## 
## $value
## [1] -1.832126
## 
## $counts
## function gradient 
##       36       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

Metode Kuadrat Terkecil dengan Optim

Metode kuadrat terkecil pada regresi adalah mencari parameter \(\beta\) sehingga jumlah kuadrat sisaan adalah minimum. Jadi meminimumkan dari fungsi kuadrat sisaan. Sehingga objective function-nya adalah kuadrat sisaan.

f <- function (para,y,x){
  X <- cbind (1,x)
  yhat <- X %*% as.matrix(para)
  sisa2 <- sum((y-yhat)^2) # kuadrat sisaan
  return(sisa2) 
}
# simulasi
x1 <- runif (10,1,10)
x2 <- runif (10,1,10)
galat <- rnorm (10,0,0.5)
y <- 1+2*x1+3*x2+galat # untuk menduga nilai beta
hasil <- optim (c(1,1,1),f,y=y,x= cbind (x1,x2))
hasil$par
## [1] 0.6582865 2.1109194 2.9339999
lm(y~x1+x2)
## 
## Call:
## lm(formula = y ~ x1 + x2)
## 
## Coefficients:
## (Intercept)           x1           x2  
##      0.6634       2.1106       2.9336

Contoh:

Lakukan pendugaan paramater regresi dengan meminimumkan jumlah kuadrat galat (residual sum of square) dari data berikut! Kemudian bandingkan hasilnya dnegan output fungsi lm!

x y
1 1
2 3
3 5
4 6
5 8
6 12 
data5=data.frame(x=c(1,2,3,4,5,6),
                 y=c(1,3,5,6,8,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)
hasil2 <- lm(y~x, data = data5)

plot(data5)
abline(hasil1$par,col=4)

hasil1$par
## [1] -1.266302  2.028449
hasil2$coefficients
## (Intercept)           x 
##   -1.266667    2.028571
hasil1$value
## [1] 2.819048
sum(hasil2$residuals^2)
## [1] 2.819048
# Menggunakan fungsi LM
model1 <- lm(y~x, data=data5)
summary(model1)
## 
## Call:
## lm(formula = y ~ x, data = data5)
## 
## Residuals:
##       1       2       3       4       5       6 
##  0.2381  0.2095  0.1810 -0.8476 -0.8762  1.0952 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -1.2667     0.7815  -1.621 0.180388    
## x             2.0286     0.2007  10.109 0.000539 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8395 on 4 degrees of freedom
## Multiple R-squared:  0.9623, Adjusted R-squared:  0.9529 
## F-statistic: 102.2 on 1 and 4 DF,  p-value: 0.000539

MLE (Maximum Likelihood Estimator)

Metode ini adalah metode yang paling sering digunakan untuk menduga parameter sebaran. Mencari nilai maksimum dari suatu fungsi tujuan yang berupa fungsi likelihood-nya. Untuk mendapatkan likelihood merupakan perkalian dari fungsi sebaran, lebih mudah menduga teta dengan mentransformasi fungsi likelihood menjadi negatif log likelihood

Sebagai contoh, jika \(x_1, x_2, ..., x_n\) berasal dari peubah acak \(X \sim N(\mu, \sigma)\), tentukan penduga \(\mu\) dan \(\sigma\) menggunakan MLE. Jika log likehood saja meminimumkan, maka menjadi negatif log likelihood suapaya dimaksimumkan.

negloglik <- function (para,xd){
  nilai <- -1* sum (dnorm (xd, mean =para[1], sd=para[2], log= TRUE)) # penjumlahan fungsi likelihood distribusi normal
  return (nilai)
}

x <- rnorm (10,2,5)

hasil <- optim (c(1,1),negloglik,xd=x)
cat ("Hasil MLE:", hasil $par)
## Hasil MLE: -0.7738515 3.884483
c(mean(x), sd(x)) # pembanding
## [1] -0.7733304  4.0954121

Pendugaan untuk regresi linear dengan menggunakan metode kemungkinan maksimum.

Sebagai contoh, akan diduga parameter dari regresi linear sederhana dengan menggunakan MLE. Contoh soal sebelumnya akan digunakan kembali.

x <- matrix(c(1,2,3,4,5,6))
y <- matrix(c(1,3,5,6,8,12))

X <- cbind(1,x) # membuat maxtrix [1 x1]

Fungsi log likelihood dari Regresi Linear Sederhana adalah sebagai berikut:

\(\begin{aligned} \mathrm{log}L(\beta_0, \beta_1, \sigma^2) = -\frac{n}{2}\mathrm{log}(2\pi) - \frac{n}{2}\mathrm{log}(\sigma^2) - \frac{1}{2\sigma^2}\sum_{i=1}^n (y_i - \beta_0 - \beta_1x_i)^2 \end{aligned}\)

Cara pertama:

regloglik <- function (beta, y, X) {
  n <- nrow(X) # jumlah baris pada matrix x
  k <- ncol(X) # jumlah kolom pada matrix x
  beta0 <- beta[1]
  beta1 <- beta[2]
  var.lin <- beta[3]
  e <- y - X %*% beta[1:2] # menghitung sisaan
  logl <- -0.5*n*log(2*pi) - 0.5*n*log(var.lin) - ((t(e)%*%e)/(2*var.lin))
  return(-logl)
}
par.reg <- optim(c(1,1,1), regloglik, y=y, X=X)
par.reg$par
## [1] -1.2669724  2.0286218  0.4698408

Cara kedua:

regloglik <- function (para,y, x){
  # fungsi terdiri dari 3 parameter dalam argumen para
  # beta0 = para [1]
  # beta1 = para [2]
  # sigma = para [3]
  # nilai tengah dari regresi linear adalah beta0 + b1x1
  # standar deviasi adalah sigma
  nilai <- -1* sum (dnorm (y, mean =para[1]+ x*para[2], sd=para[3], log= TRUE))
  return (nilai)
}  

x=c(1,2,3,4,5,6)
y=c(1,3,5,6,8,12)

hasil <- optim (c(1,1,1),regloglik,y =y, x=x)
cat ("Hasil MLE:", hasil $par)
## Hasil MLE: -1.266469 2.028513 0.6853899

  1. Mahasiswa Pascasarjana IPB, ↩︎