Diferensial

Untuk mendapatkan turunan dari suatu fungsi pada pemrograman R dapat menggunakan fungsi D(expr, simbol atau deriv(~fungsi, simbol). Fungsi D dapat digunakan jika hasil turunan merupakan suatu fungsi. Fungsi deriv dapat digunakan jika akan memasukkan nilai dari hasil turunan pada suatu fungsi. Sebagai contoh, akan dicari turunan dari fungsi-fungsi berikut.

\[\mathrm{f(x)} = exp (x)^2\]

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

exp(x^2) * (2 * x)

\[\mathrm{f(x)} = x^2\], untuk nilai x = 2

xturunan <- deriv (~x^2, "x")
x<-2
eval(xturunan)

[1] 4
attr(,"gradient")
     x
[1,] 4

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.

Berikut akan diberikan contoh penggunakan fungsi integrate pada beberapa kasus.

\[\int{ x^2} dx\]

fs <- function (x)x^2
integrate(fs,0,1)

0.3333333 with absolute error < 3.7e-15

\[\int{ x^2} + 4x dx\]

f1 <- function(x) x^2 + 4*x
integrate(f1, lower=-10, upper=10)

666.6667 with absolute error < 7.6e-12

\[\int{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

Berikut akan diberikan contoh penggunakan fungsi yac_str pada beberapa kasus.

\[\int{ x^2} + 4x dx\]

library(Ryacas)
yac_str("Integrate(x) x^2 + 4*x")

[1] "x^3/3+2*x^2"

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

Golden Section Search

Metode golden section search digunakan untuk mencari nilai minimum suatu fungsi yang dibatasi dari dua buah nilai, yaitu sebuah selang a dan b. Secara umum, jika harus menyertakan dua buah nilai terkadang tidak diketahui apakah di selang tersebut terdapat nilai maksimum atau kah minimumnya. Untuk mengatasi permasalahan tersebut terdapat teknik newton raphson.

Algoritma untuk teknik ini adalah sebagai berikut:

• Mulai dengan selang [a,b] yang memuat minimum
  
• Perkecil selang [a’, b’] yang memuat minimum
  
• Berhenti sampai |b’ - a’| lebih kecil dari nilai tolerans

Pemilihan nilai a’ dan b’ adalah sebagai berikut:

• Nilai antara [a,b] memiliki sifat golden ratio
  
• Tentukan \(x_1\) dan \(x_2\)"
  \(x_1\) \(= b - (b-a)/golden ratio\)
  \(x_2\) \(= a + (b-a)/golden ratio\)
  
• Hitung f(\(x_1\)) dan f(\(x_2\))
• Jika f(\(x_1\)) > f(\(x_2\)) maka [a’, b’] = [\(x_1\), b]
  
• Jika f(\(x_1\)) < f(\(x_2\)) maka [a’, b’] = [a, \(x_2\)]

Membuat fungsi golden section search

golden <- function (f, a, b, tol=0.0000001){
  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 (f2>f1){
      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)
  }

\[\mathrm{f(x)} = |x-3.5| + (x-2)^2\]

# Membuat fungsi f(x)
f <- function(x) {abs(x-3.5)+(x-2)^2}

# Membuat plot
curve(f, 1,5)

# Menghitung nilai optimum
golden(f,1,2)

[1] 2

golden(f,1,5)

[1] 2.5

golden(f,3,5)

[1] 3

Jika dilihat secara plot grafik, maka terlihat nilai minimum ada di antara nilai x= 2 dan x = 3. Dengan menggunakan fungsi golden section search dengan interval 1 <x<5 terlihat bahwa nilai minimum berada di x = 2.5.

\[\mathrm{f(x)} = 2 sin x - \frac {x^2} {10}\]

# Membuat fungsi f(x)
f <- function(x) {(2 * sin(x)) - ((x^2) /10)}
fa <- function(x) -1 * f(x) # fungsi dikali dengan -1 karena fungsi yang dibuat meminimumkan

# Membuat plot
curve(f, 0,4)

# Menghitung nilai optimum
golden(fa,0,4)

[1] 1.427552

Jika dilihat dari plot tersebut, pada interval 0 < x < 4 nilai optimum terlihat di antara x = 1 dan x = 2, jika menggunakan fungsi golden section search ditemukan bahwa nilai x yang memaksimumkan fungsi adalah pada saat x = 1.47

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.

Membuat 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.
\[\mathrm{f(x)} = 4x^2-3x-7\]

fx <- expression (4*x^2-3*x-7)
newtonr(fx,3)

[1] 0.375

\[\mathrm{f(x)} = e^-x+x^4\]

fx <- expression (exp(-x)+x^4)
newtonr(fx)

[1] 0.5282519

\[\mathrm{f(x)} = x^2-x\]

fx <- expression (x^2-x)
newtonr(fx)

[1] 0.5

Fungsi Optimasi Bult-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, dan optim untuk lebih dari satu peubah.

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

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

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
optim (c(-1.2,1),fr) # argumen pertama adalah nilai inisial, karena menduga x vektor berukuran 2 maka dimasukkan nilai inisialnya

$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

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

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

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.8897314 1.9060933 3.0902045

lm(y~x1+x2)

Call:
lm(formula = y ~ x1 + x2)

Coefficients:
(Intercept)           x1           x2  
      0.893        1.906        3.090  

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. Bisa menggunakan logaritma karena merupakan fungsi monoton.

Sebagai contoh, jika \(x_1\), \(x_2\), …, \(x_n\) berasal dari peubah acak X~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: 5.042608 4.584582

cat ("Hasil Manual:",c(mean (x),sd(x))) # pembanding

Hasil Manual: 5.046133 4.833554

Dapat dilihat pendugaan \(\mu\) dan \(\sigma\) menggunakan MLE atau pun perhitungan manual memberikan hasil yang sama, yaitu \(\mu\) = 2.64 dan \(\sigma\) = 5.13.

Pendugaan untuk distribusi eksponensial, dengan menggunakan fungsi mle pada library stats4. Fungsi loglikelihood dari distribusi eksponensial adalah sebagai berikut.

library(stats4)
set.seed(1)
xexp = rexp(n=500, rate = 1) # pembangkitan bilangan acak
lexp = function(rate) {
  n = length(xexp)
  x = xexp
  flog <- (n*log(rate)-rate*sum(x))
  return(-flog)
}
estexp = mle(minuslogl = lexp, start=list(rate=1))
summary(estexp)

Maximum likelihood estimation

Call:
mle(minuslogl = lexp, start = list(rate = 1))

Coefficients:
     Estimate Std. Error
rate    1.047 0.04682324

-2 log L: 954.0697 

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. Misalkan ada sebuah data pengamatan antara respon dan peubah penjelas sebagai berikut:

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:

Cara pertama (penjabaran):

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