Praktikum Pemrograman Statistika

Optimize It Right: Pendekatan Numerik untuk Menemukan Yang Terbaik

Mengapa Kita Butuh Optimasi?

Di balik berbagai metode statistik —mulai dari regresi logistik hingga estimasi parameter berbasis likelihood—ada satu proses fundamental yang sering luput dari perhatian: optimasi numerik. Optimasi adalah upaya mencari titik terbaik (minimum atau maksimum) dari suatu fungsi. Sayangnya, tidak semua fungsi dapat diselesaikan secara analitik atau melalui turunan manual, terutama jika bentuknya kompleks atau tidak memiliki solusi eksak. Dalam situasi seperti ini, optimasi numerik menjadi solusi: pendekatan iteratif berbasis komputer untuk menemukan nilai optimal suatu fungsi.

Dalam statistika, optimasi muncul hampir di setiap langkah penting, misalnya:

• Menentukan koefisien regresi yang meminimalkan squared error,
  
• Memaksimumkan fungsi likelihood dalam Maximum Likelihood Estimation (MLE),
  
• Melatih model machine learning untuk menekan nilai fungsi loss sekecil mungkin.

Pertanyaannya, bagaimana kita menemukan titik optimal dari fungsi yang bahkan sulit digambarkan grafiknya? Di sinilah optimasi numerik memainkan perannya: menyediakan mekanisme komputasional untuk menjelajahi ruang parameter dan mendekati solusi terbaik ketika metode kalkulus tidak lagi memadai.

Melalui praktikum ini, kita tidak hanya akan memanfaatkan fungsi-fungsi optimasi dalam R, tetapi juga membangun pemahaman konseptual dan intuisi di baliknya. Dengan begitu, kita tidak sekadar menjalankan kode, melainkan mengerti apa yang sedang dicari dan bagaimana proses pencariannya berlangsung.

Representasi Bilangan dan Presisi Komputasi

Sebelum masuk ke algoritma optimasi, kita perlu menyadari satu hal mendasar: komputer tidak menyimpan bilangan seperti cara kita menuliskannya di kertas. Berbeda dengan notasi matematika yang bisa menuliskan bilangan hingga 100 atau 1000 digit di atas kertas, komputer menyimpan angka dalam bentuk floating-point berbasis biner, yang hanya mampu menyimpan sejumlah digit terbatas. Akibatnya, banyak angka yang hanya bisa direpresentasikan secara hampiran, bukan secara eksak. Bahkan angka sesederhana 0.1 atau 0.2 atau akar dari 2 tidak direpresentasikan secara eksak di dalam memori komputer, tetapi dalam bentuk pendekatan.

Sebagai ilustrasi:

0.1 + 0.2 == 0.3

[1] FALSE

Hasilnya mungkin mengejutkan: FALSE, padahal secara matematis 0.1 + 0.2 = 0.3. Mengapa demikian? Karena 0.1 + 0.2 secara internal tidak persis sama dengan 0.3, akibat pembulatan biner dalam penyimpanan.

Kita dapat menggunakan fungsi seperti print() dan sprintf() untuk mengeksplorasi ketelitian penyimpanan angka:

x <- sqrt(2)^2
print(x)

[1] 2

sprintf("%.20f", x)

[1] "2.00000000000000044409"

Dari sini terlihat bahwa meskipun secara konseptual kita mengharapkan hasil 2, namun pada level penyimpanan bit, masih terdapat perbedaan kecil yang bisa berpengaruh besar dalam proses komputasi berulang—seperti yang terjadi dalam algoritma optimasi.

Fenomena-fenomena di atas penting dalam optimasi numerik. Ketika algoritma berhenti saat perubahan nilai sangat kecil, bisa jadi yang dianggap “cukup kecil” oleh komputer sebenarnya belum cukup bagi model. Sebaliknya, bisa juga algoritma tidak pernah berhenti karena perubahan nilainya terus-menerus dianggap signifikan—padahal itu hanya noise dari presisi penyimpanan. Pada intinya:

• Jangan anggap komputer menghitung dengan sempurna.
• Jangan terlalu percaya pada angka di layar: ada batas presisi di baliknya.
• Saat algoritma optimasi gagal “konvergen”, bisa jadi bukan rumus atau modelnya yang keliru, tetapi batasan presisi komputasi yang menghalangi algoritma menemukan titik optimal.

Fungsi Kalkulus: Diferensial

Optimasi numerik sangat erat kaitannya dengan turunan dalam kalkulus. Melalui turunan, kita bisa mengetahui bagaimana suatu fungsi berubah ketika variabel inputnya berubah. Inilah yang menjadi dasar untuk mencari titik minimum atau maksimum suatu fungsi.

Turunan Satu Variabel

Untuk fungsi satu variabel \(f(x)\), turunan pertama \(f'(x)\) menunjukkan laju perubahan fungsi terhadap \(x\). Secara intuitif:

• Jika \(f'(x) > 0\), maka fungsi sedang meningkat.
• Jika \(f'(x) < 0\), maka fungsi sedang menurun.
• Jika \(f'(x) = 0\), kita berada di titik stasioner.

Titik stasioner inilah yang menjadi kandidat untuk minimum atau maksimum.

Misalnya, untuk fungsi kuadrat:

\[ f(x) = x^2 - 4x + 3 \]

Turunan pertamanya adalah:

\[ f'(x) = 2x - 4 \]

Dengan menyelesaikan persamaan \(f'(x) = 0\), kita peroleh titik stasioner \(x = 2\). Nilai \(f(2) = -1\) adalah titik minimum/maksimum dari fungsi tersebut.

Turunan Kedua

Selain turunan pertama, turunan kedua juga penting dalam optimasi.

• Jika ( f’’(x) > 0 ), maka titik stasioner adalah minimum lokal.
  
• Jika ( f’’(x) < 0 ), maka titik stasioner adalah maksimum lokal.
  
• Jika ( f’’(x) = 0 ), kita perlu analisis lebih lanjut (bisa jadi titik belok).

Dengan pemahaman konsep dasar diferensial ini, kita akan lebih mudah mengerti bagaimana algoritma optimasi bekerja.

Menghitung Turunan di R

Selain memahami konsep turunan secara matematis, kita juga bisa menghitung turunan langsung dalam R. R menyediakan fungsi bawaan untuk melakukan diferensiasi simbolik:

• D(expr, simbol) yang digunakan jika kita ingin menurunkan suatu fungsi secara simbolik.
  
• deriv(~fungsi, simbol) yang digunakan jika kita ingin menghitung nilai turunan pada fungsi tertentu.

Misalkan kita ingin mencari turunan dari:

\[ f(x) = x^2 - 4x + 3 \]

Kita bisa menghitung turunannya dengan D():

# Turunan simbolik dengan D()
D(expression(x^2 - 4*x + 3), "x")

2 * x - 4

Hasil dari perintah di atas adalah 2 * x - 4. Sehingga turunan fungsi dapat dituliskan sebagai:

\[ f'(x) = 2x - 4 \]

Untuk mencari titik stasioner, kita perlu menyelesaikan persamaan \(f'(x) = 0\):

\[ 2x - 4 = 0 \quad \Rightarrow \quad x = 2 \]

Selain D(), kita juga bisa menggunakan fungsi deriv() untuk mendapatkan turunan sekaligus mengevaluasi nilainya pada titik tertentu.

# Turunan pertama dengan deriv()
turunan1 <- deriv(~x^2 - 4*x + 3, "x")

# Evaluasi pada x = 2
eval(turunan1, list(x = 2))

[1] -1
attr(,"gradient")
     x
[1,] 0

Berdasarkan hal tersebut, kita mendapatkan nilai fungsi pada \(x=2\), yaitu \(f(2)=−1\) dan nilai turunan pertama pada \(x=2\), yaitu \(f'(2)=0\).

Untuk memperoleh turunan kedua, kita cukup menerapkan D() kembali pada hasil turunan pertama:

# Turunan kedua secara simbolik
D(D(expression(x^2 - 4*x + 3), "x"), "x")

[1] 2

Hasilnya adalah \(f''(x)=2\). Karena \(f''(2)=2>0\), maka titik \((2,−1)\) dipastikan merupakan titik minimum dari fungsi \(f(x)=x^2−4x+3\).

Fungsi Kalkulus: Integral

Selain diferensial, konsep penting lain dari kalkulus yang berperan dalam optimasi adalah integral. Jika diferensial menggambarkan laju perubahan suatu fungsi, maka integral menggambarkan akumulasi perubahan tersebut.

Secara sederhana, integral dapat dipahami sebagai luas daerah di bawah kurva suatu fungsi. Untuk fungsi \(f(x)\), integral tak tentu dituliskan sebagai:

\[ \int f(x)\,dx \]

Sedangkan integral tentu, yang dihitung pada selang tertentu \([a,b]\), dituliskan sebagai:

\[ \int_a^b f(x)\,dx \]

Integral tentu ini memiliki interpretasi geometris sebagai luas daerah di bawah kurva \(f(x)\) dari \(x=a\) hingga \(x=b\).

Integral Tentu dengan integrate()

R menyediakan fungsi integrate(fungsi, lower, upper) untuk menghitung integral tentu dari suatu fungsi pada interval tertentu.

Contoh:

\[ \int_0^1 x^2 \, dx \]

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

0.3333333 with absolute error < 3.7e-15

Contoh lain:

\[ \int_{-10}^{10} \left( x^2 + 4x \right) \, dx \]

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

666.6667 with absolute error < 7.6e-12

Untuk kasus integral tak hingga, misalnya:

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

Integral Tak Tentu dengan Ryacas

Untuk menghitung integral secara simbolik, kita bisa menggunakan fungsi yac_str() dari library Ryacas.

library(Ryacas)

Warning: package 'Ryacas' was built under R version 4.4.3

Attaching package: 'Ryacas'

The following object is masked from 'package:stats':

    integrate

The following objects are masked from 'package:base':

    %*%, det, diag, diag<-, lower.tri, upper.tri

# Integral tak tentu dari x^2 + 4x
yac_str("Integrate(x) x^2 + 4*x")

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

Contoh lain:

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

Dengan integrate() kita bisa menghitung integral tentu secara numerik, sedangkan Ryacas membantu kita mencari bentuk simboliknya. Keduanya melengkapi pemahaman kita tentang integral dalam konteks komputasi.

Mari kita lihat contoh fungsi berikut:

\[ f(x) = x^2 + 4x \]

Kita cari integral tak tentu terlebih dahulu:

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

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

Hasilnya adalah

\[ \int (x^2 + 4x)\, dx = \frac{x^3}{3} + 2x^2 + C \]

Kemudian, jika batas integralnya dari \(-10\) hingga \(10\), maka:

\[ \int_{-10}^{10} (x^2 + 4x)\, dx = \left[ \frac{x^3}{3} + 2x^2 \right]_{-10}^{10} \]

Kita selesaikan dengan menghitung secara manual

\[ \left( \frac{10^3}{3} + 2 \cdot 10^2 \right) - \left( \frac{(-10)^3}{3} + 2 \cdot (-10)^2 \right) \]

\[ = \left( \frac{1000}{3} + 200 \right) - \left( -\frac{1000}{3} + 200 \right) = \frac{2000}{3} \]

atau sekitar \(666.67\).

Sekarang kita cek hasilnya dengan integrate():

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

666.6667 with absolute error < 7.6e-12

Dengan demikian, baik perhitungan simbolik (Ryacas) maupun numerik (integrate()) menghasilkan hasil yang konsisten.

Optimasi Numerik

Banyak metode statistik bergantung pada proses mencari nilai optimum dari suatu fungsi tujuan. Misalnya:

• Metode kemungkinan maksimum (Maximum Likelihood Estimation, MLE): mencari nilai yang memaksimumkan fungsi likelihood,
  
• Metode kuadrat terkecil (Least Squares): mencari nilai yang meminimalkan jumlah kuadrat galat.

Proses mendapatkan nilai optimum ini dikenal dengan istilah optimasi numerik.

Untuk fungsi satu variabel, misalkan \(y = f(x)\), tujuan optimasi adalah mencari nilai \(x\) yang menghasilkan \(y\) maksimum atau minimum. Nilai \(x\) tersebut kemudian dianggap sebagai solusi optimal dari permasalahan.

Berbagai metode optimasi telah dikembangkan, mulai dari teknik sederhana hingga algoritma kompleks. Dua metode klasik yang akan dibahas lebih dulu adalah:

• Golden Section Search, yang mencari titik minimum fungsi satu variabel dengan memperkecil interval pencarian secara bertahap,
  
• Newton–Raphson, yang memanfaatkan informasi turunan pertama dan kedua untuk mempercepat pencarian titik optimum.

Kedua metode ini memberikan fondasi awal untuk memahami bagaimana komputer “menjelajahi” ruang parameter dalam rangka menemukan solusi optimal.

Golden Section Search

Golden Section Search adalah metode optimasi numerik sederhana untuk mencari titik minimum dari suatu fungsi satu variabel pada interval tertentu. Metode ini tidak memerlukan turunan, hanya membutuhkan nilai fungsi pada titik-titik tertentu.

Ide dasarnya adalah mempersempit interval pencarian secara bertahap dengan memanfaatkan rasio emas (golden ratio):

\[ \varphi = \frac{\sqrt{5}-1}{2} \approx 0.618 \]

Langkah algoritma secara umum:

1. Tentukan fungsi yang akan diminimalkan
   Misalkan \(f(x)\) didefinisikan pada interval \([a,b]\).

2. Gunakan sifat rasio emas
   Rasio emas didefinisikan sebagai:

   \[ \varphi = \frac{\sqrt{5} + 1}{2} \approx 1.618 \]

3. Hitung dua titik dalam interval
   \[ x_1 = b - \frac{b-a}{\varphi}, \quad x_2 = a + \frac{b-a}{\varphi} \]

4. Evaluasi fungsi pada kedua titik
   Hitung \(f(x_1)\) dan \(f(x_2)\).

5. Perbarui interval pencarian

   • Jika \(f(x_1) > f(x_2)\), maka interval baru adalah \([x_1, b]\).
     
   • Jika \(f(x_1) < f(x_2)\), maka interval baru adalah \([a, x_2]\).

6. Ulangi langkah 3–5 sampai selisih dari panjang interval \(|b-a|\) lebih kecil dari toleransi yang ditentukan.

7. Ambil solusi optimum. Titik minimum hampiran diperoleh sebagai rata-rata dari ujung interval terakhir:

   \[ x^\ast = \frac{a+b}{2} \]

Berdasarkan algoritma Golden Section Search tersebut, kita dapat membuat fungsi di R untuk mencari titik minimum suatu fungsi satu variabel. Fungsi ini akan menerima masukan berupa fungsi \(f(x)\), interval awal \([a,b]\), dan toleransi penghentian iterasi.

golden_search <- function(f, a, b, tol = 1e-6) {
  ratio <- (sqrt(5) + 1) / 2  # golden ratio
  
  # Inisialisasi dua titik awal di dalam interval [a, b]
  x1 <- b - (b - a) / ratio
  x2 <- a + (b - a) / ratio
  
  # Ulangi sampai panjang interval lebih kecil dari toleransi
  while (abs(b - a) > tol) {
    # Bandingkan nilai fungsi di kedua titik
    if (f(x1) > f(x2)) {
      # Minimum ada di kanan, geser interval ke [x1, b]
      a <- x1
      x1 <- x2
      x2 <- a + (b - a) / ratio
    } else {
      # Minimum ada di kiri, geser interval ke [a, x2]
      b <- x2
      x2 <- x1
      x1 <- b - (b - a) / ratio
    }
  }
  
  # Titik minimum hampiran adalah titik tengah interval terakhir
  return((a + b) / 2)
}

Kita gunakan fungsi berikut untuk ilustrasi:

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

# Definisi fungsi
f <- function(x) abs(x - 3.5) + (x - 2)^2

Mari kita lihat bentuk kurva fungsi \(f(x)\) pada rentang \([0,6]\).

curve(f, from = 0, to = 6, n = 1000,
      main = "Kurva f(x) = |x - 3.5| + (x - 2)^2",
      xlab = "x", ylab = "f(x)")

Dari grafik terlihat bahwa fungsi memiliki minimum global di sekitar \(x=2.5\)

Sekarang kita coba gunakan golden_search() dengan interval berbeda-beda.

# Interval luas yang mencakup minimum
golden_search(f, 0, 6)

[1] 2.5

Algoritma menemukan titik minimum global di \(x^\ast=2.5\)

# Interval yang masih mencakup minimum
golden_search(f, 2, 5)

[1] 2.5

Hasilnya sama, karena interval masih mencakup titik minimum global.

# Interval yang tidak mencakup minimum (salah inisiasi)
golden_search(f, 0, 2)

[1] 2

Karena interval tidak mencakup titik minimum global (2.5), algoritma berhenti di batas kanan interval, yaitu \(x=2\).

Percobaan pada beberapa interval di atas menunjukkan bahwa penentuan interval awal sangat memengaruhi hasil Golden Section Search.

• Jika interval mencakup titik optimum global (misalnya \([0,6]\) atau \([2,5]\)), algoritma akan menemukan titik minimum global di sekitar \(x^\ast = 2.5\).
• Jika interval tidak mencakup titik optimum global (misalnya \([0,2]\)), algoritma tetap berhenti, tetapi hasilnya bukan solusi global. Pada kasus ini, algoritma berhenti di batas kanan \(x=2\), dengan nilai \(f(2) = 1.5\).

Oleh karena itu, pemilihan interval awal \([a,b]\) menjadi tahap krusial dalam algoritma berbasis pencarian interval. Cara praktis untuk menentukan interval awal yang tepat adalah melakukan visualisasi terlebih dahulu, lalu menentukan interval secara intuitif berdasarkan bentuk kurva yang terlihat. Dengan pendekatan ini, kita dapat memastikan bahwa Golden Section Search bekerja sesuai tujuan, yaitu menemukan titik minimum global dari fungsi yang diteliti.

Berikut contoh visualisasi fungsi dengan dua interval berbeda:

• Interval benar: \([0,6]\), yang mencakup titik minimum global.
  
• Interval salah: \([0,2]\), yang tidak mencakup titik minimum global.

# Plot fungsi
curve(f, from = 0, to = 6, n = 1000,
      main = "Pengaruh Pemilihan Interval Awal",
      xlab = "x", ylab = "f(x)")

# Tandai minimum global di sekitar x=2.5
abline(v = 2.5, col = "red", lty = 2)
text(2.5, 3, "Minimum global", pos = 4, col = "red")

# Tandai interval benar [0,6]
rect(0, 0, 6, max(f(seq(0,6,0.01))), border = "darkgreen", lwd = 2)

# Tandai interval salah [0,2]
rect(0, 0, 2, max(f(seq(0,2,0.01))), border = "blue", lwd = 2)

legend("topright",
       legend = c("Minimum global", "Interval benar [0,6]", "Interval salah [0,2]"),
       col = c("red", "darkgreen", "blue"),
       lty = c(2,1,1), lwd = c(1,2,2), bty = "n")

Sejauh ini kita sudah bisa menemukan nilai \(x\) yang meminimumkan fungsi. Namun, dalam praktik seringkali kita juga ingin tahu berapa nilai minimum fungsi itu sendiri pada titik optimum tersebut. Mari kita coba sekaligus menghitung keduanya.

# Gunakan golden_search untuk mencari minimum global di [0,6]
xmin <- golden_search(f, 0, 6)

# Hitung nilai fungsi pada titik tersebut
fmin <- f(xmin)

xmin   # Titik x yang meminimumkan fungsi

[1] 2.5

fmin   # Nilai minimum fungsi

[1] 1.25

Dari perhitungan, kita peroleh titik minimum di \(x = 2.5\) dengan nilai minimum dari fungsi \(f(x)\) adalah \(1.25\).

Metode Golden Section Search pada dasarnya dirancang untuk mencari nilai minimum dari suatu fungsi. Lalu, bagaimana jika kita ingin mencari nilai maksimum?. Caranya sangat sederhana, kita cukup meminimalkan fungsi negatif \(-f(x)\). Dengan begitu, titik minimum dari \(-f(x)\) akan menjadi titik maksimum dari \(f(x)\).

Misalkan kita ingin mencari nilai maksimum dari fungsi berikut:

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

# Definisi fungsi
f <- function(x) 2*sin(x) - (x^2)/10

Mari kita lihat bentuk kurva fungsi \(f(x)\) pada rentang \([0,5]\).

curve(f, from = 0, to = 5, n = 1000,
      main = "Kurva f(x) = 2*sin(x) - x^2/10",
      xlab = "x", ylab = "f(x)")

Dari grafik terlihat bahwa fungsi mencapai puncak (maksimum) di sekitar \(x≈1.4\) atau antara \(1\) sampai \(2\). Oleh karena itu, kita bisa memilih interval \([0,2]\) yang mencakup titik maksimum tersebut. Kita juga dapat memilih interval yang lebih lebar jika kita ragu apakah suatu interval memuat titik optimum dari suatu fungsi.

Karena kita mencari maksimum, kita definisikan fungsi negatifnya:

f_neg <- function(x) -1 * f(x)

# Cari titik maksimum pada interval [1, 2]
xmax <- golden_search(f_neg, 1, 2)

# Hitung nilai maksimumnya
fmax <- f(xmax)

xmax   # Titik maksimum hampiran

[1] 1.427552

fmax   # Nilai maksimum fungsi

[1] 1.775726

Dengan modifikasi sederhana ini, Golden Section Search dapat digunakan bukan hanya untuk mencari minimum, tetapi juga untuk mencari maksimum dengan meminimalkan fungsi negatifnya.

Newton-Raphson Method

Metode Newton-Raphson adalah salah satu algoritma numerik yang digunakan untuk mencari titik optimum dari suatu fungsi.
Prinsip dasarnya adalah mencari titik \(x\) yang membuat turunan pertama fungsi bernilai nol (\(f'(x) = 0\)), dengan bantuan turunan kedua untuk memperbaiki tebakan.

Newton-Raphson menggunakan pendekatan deret Taylor untuk melakukan iterasi:

\[ x_{k+1} = x_k - \frac{f'(x_k)}{f''(x_k)} \]

• \(f'(x)\) : turunan pertama,
  
• \(f''(x)\) : turunan kedua,
  
• \(x_k\) : tebakan awal,
  
• Iterasi berlanjut sampai selisih antar iterasi cukup kecil.

Kelebihan metode Newton-Raphson adalah hanya memerlukan satu nilai inisialisasi (\(x_0\)) untuk memulai iterasi. Selain itu, metode ini memiliki laju konvergensi yang sangat cepat apabila tebakan awal cukup dekat dengan titik optimum. Namun, kelemahannya adalah metode ini membutuhkan turunan pertama dan turunan kedua dari fungsi, sehingga tidak bisa digunakan jika fungsi sulit untuk diturunkan. Newton-Raphson juga sensitif terhadap tebakan awal; jika tebakan awal kurang tepat, metode dapat gagal konvergen atau menuju ke titik yang salah. Dibandingkan dengan Golden Section Search, metode Newton-Raphson umumnya lebih cepat, tetapi kurang fleksibel karena adanya syarat keberadaan turunan pertama dan kedua.

Algoritma dari metode Newton-Raphson adalah

1. Tentukan fungsi \(f(x)\) dan hitung turunannya \(f'(x)\) dan \(f''(x)\).
   
2. Pilih tebakan awal \(x_0\).
   
3. Hitung nilai baru menggunakan rumus:
   \[ x_{k+1} = x_k - \frac{f'(x_k)}{f''(x_k)} \]
   
4. Ulangi sampai konvergen ((|x_{k+1} - x_k| < )).
5. Ambil (x_{k+1}) sebagai titik minimum hampiran

Berdasarkan algoritma tersebut, kita dapat membuat fungsi di R untuk mencari titik minimum suatu fungsi satu variabel. Fungsi ini akan menerima masukan berupa turunan pertama ((f’(x))), turunan kedua ((f’’(x))), tebakan awal ((x_0)), serta toleransi konvergensi.

newton_raphson <- function(fx, x0 = 1, tol = 1e-6) {
  # Hitung turunan pertama dan kedua dari fungsi simbolik fx
  f_prime <- deriv(fx, "x")                # turunan pertama f'(x)
  f_double <- deriv(D(fx, "x"), "x")       # turunan kedua f''(x)
  
  # Inisialisasi error
  error <- 1000
  
  # Iterasi Newton-Raphson
  while (error > tol) {
    x <- x0
    
    # Evaluasi turunan pertama dan kedua pada x
    f1 <- attr(eval(f_prime, list(x = x)), "gradient")[1]
    f2 <- attr(eval(f_double, list(x = x)), "gradient")[1]
    
    # Perbarui titik menggunakan formula Newton-Raphson
    x_new <- x - f1 / f2
    
    # Update error (selisih antara titik baru dan lama)
    error <- abs(x_new - x)
    
    # Update nilai awal untuk iterasi berikutnya
    x0 <- x_new
  }
  
  # Hasil titik optimum hampiran
  return(x0)
}

Dalam implementasi di atas, perulangan Newton-Raphson akan terus berjalan sampai nilai error \(|x_{k+1} - x_k| < \varepsilon\) terpenuhi. Namun, jika fungsi atau tebakan awal yang dipilih menyebabkan algoritma tidak konvergen, maka perulangan tidak akan pernah berhenti (infinite loop).

Kondisi ini dapat terjadi, misalnya karena:
- Turunan kedua \(f''(x)\) mendekati nol sehingga pembaruan \[ x_{k+1} = x_k - \frac{f'(x_k)}{f''(x_k)} \] menjadi tidak stabil,
- Tebakan awal \(x_0\) terlalu jauh dari titik optimum sehingga iterasi menyimpang,
- Fungsi tidak memiliki bentuk yang sesuai untuk dicari optimum dengan Newton-Raphson.

Untuk mengatasi hal tersebut, dalam praktik biasanya ditambahkan batas jumlah iterasi maksimum. Jika jumlah iterasi telah mencapai batas tersebut tetapi belum konvergen, maka algoritma dihentikan secara paksa dan dianggap gagal menemukan solusi. Dengan cara ini, kita dapat menghindari masalah infinite loop dan menjaga agar program tidak berjalan tanpa henti.

Misalkan kita ingin mencari nilai minimum dari fungsi:

\[ f(x) = 4x^2 - 3x - 7 \]

Pertama, kita definisikan fungsi dalam R:

fx1 <- expression(4*x^2 - 3*x - 7)

Selanjutnya, kita gunakan metode Newton-Raphson untuk mencari titik minimum:

xmin1 <- newton_raphson(fx1, x0 = 1)   # tebakan awal x0 = 1
fmin1 <- eval(fx1, list(x = xmin1))    # nilai fungsi pada titik minimum

xmin1   # titik minimum hampiran

[1] 0.375

fmin1   # nilai minimum fungsi

[1] -7.5625

Dengan demikian, kita memperoleh titik minimum di \(x=0.375\) dengan nilai fungsi minimum \(f(x)\) adalah \(-7.5625\)

Sebagai tambahan, mari kita terapkan metode Newton-Raphson pada dua fungsi lain.

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

fx2 <- expression(exp(-x) + x^4)

xmin2 <- newton_raphson(fx2, x0 = 0)   # tebakan awal
fmin2 <- eval(fx2, list(x = xmin2))

xmin2   # titik minimum hampiran

[1] 0.5282519

fmin2   # nilai minimum fungsi

[1] 0.6675038

\[ f(x) = x^2 - x \]

fx3 <- expression(x^2 - x)

xmin3 <- newton_raphson(fx3, x0 = 0.5)   # tebakan awal
fmin3 <- eval(fx3, list(x = xmin3))

xmin3   # titik minimum hampiran

[1] 0.5

fmin3   # nilai minimum fungsi

[1] -0.25

Fungsi Optimasi Built-In

Selain menggunakan algoritma khusus seperti Golden Section Search atau Newton-Raphson, R juga menyediakan fungsi bawaan (built-in) untuk melakukan optimasi secara lebih praktis. Salah satu algoritma yang banyak digunakan adalah Nelder-Mead, yaitu metode optimasi yang cocok untuk fungsi dengan lebih dari satu variabel. Dalam R, fungsi-fungsi optimasi yang bisa kita gunakan antara lain:

• optimize() atau optimise() yang digunakan untuk mencari nilai minimum atau maksimum dari fungsi satu variabel.
  
• optim() yang digunakan untuk kasus lebih dari satu variabel, dengan algoritma default Nelder-Mead.

Fungsi optimize() atau optimise()

Untuk kasus satu peubah, R menyediakan fungsi optimize() (atau optimise()) untuk mencari nilai minimum maupun maksimum suatu fungsi.

Sebagai contoh, misalkan kita memiliki fungsi:

\[ f(x) = \left(x - \tfrac{1}{3}\right)^2 \]

Kita definisikan fungsi dalam R:

f <- function(x) (x - 1/3)^2

Selanjutnya, kita gunakan optimize() untuk mencari titik minimum pada interval tertentu, misalnya [−2,2]:

result <- optimize(f, interval = c(-2, 2))
result

$minimum
[1] 0.3333333

$objective
[1] 3.081488e-33

Hasil optimize() akan memberikan dua informasi utama, yaitu minimum (nilai \(x\) yang meminimumkan fungsi) dan objective (nilai minimum dari fungsi pada titik tersebut). Mari kita gambarkan fungsi tersebut untuk memperjelas:

curve(f, from = -2, to = 2, n = 1000,
      main = expression(f(x) == (x - 1/3)^2),
      xlab = "x", ylab = "f(x)")
abline(v = result$minimum, col = "red", lty = 2)     # garis vertikal titik minimum
points(result$minimum, result$objective, col = "blue", pch = 19)  # titik minimum

Sekarang, kita akan mencari nilai minimum dan maksimum lokal maupun global dari fungsi berikut

\[ f(x) = \sin(x) + \sin(2x) + \cos(3x) \]

Pertama, kita definisikan fungsi di R:

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

Lalu kita gambarkan kurva pada interval \([0, 2\pi]\) untuk melihat pola naik-turun fungsi:

curve(f, from = 0, to = 2*pi, n = 2000,
      main = "f(x) = sin(x) + sin(2x) + cos(3x)",
      xlab = "x", ylab = "f(x)")

Fungsi \(f(x) = \sin(x) + \sin(2x) + \cos(3x)\) merupakan kombinasi dari beberapa gelombang sinus dan kosinus dengan frekuensi yang berbeda. Hasilnya, grafik fungsi memperlihatkan pola naik-turun atau berosilasi sepanjang interval \([0, 2\pi]\). Pola osilasi ini menyebabkan munculnya beberapa titik puncak (kandidat maksimum lokal) dan beberapa titik lembah (kandidat minimum lokal). Di antara titik-titik lokal tersebut terdapat satu puncak tertinggi yang menjadi maksimum global serta satu lembah terdalam yang menjadi minimum global. Oleh karena itu, dari visualisasi kurva dapat kita simpulkan bahwa fungsi ini memiliki lebih dari satu titik ekstrim, baik minimum maupun maksimum, yang akan dianalisis lebih lanjut menggunakan metode optimasi bawaan R.

Dari grafik terlihat adanya lembah di sekitar \(x \approx 1.15\). Untuk menemukan titik ini secara numerik, kita dapat menggunakan fungsi optimize() dengan membatasi pencarian pada interval \([0.8, 1.4]\).

min_lokal <- optimize(f, interval = c(0.8, 1.4))
min_lokal$minimum      # minimum lokal

[1] 1.152872

min_lokal$objective    # f(x) minimum lokal

[1] 0.7056316

Pada grafik \(f(x)=\sin x+\sin(2x)+\cos(3x)\) terlihat sangat jelas naik-turun berkali-kali, sehingga terdapat beberapa lembah (minimum lokal) dan puncak (maksimum lokal). Jika kita langsung memanggil optimize() pada rentang penuh \([0,2\pi]\)

min_02pi <- optimize(f, interval = c(0, 2*pi))  # asumsi unimodal (tidak selalu valid di sini)
min_02pi$minimum      # x* yang ditemukan

[1] 3.033129

min_02pi$objective    # f(x*) yang ditemukan

[1] -1.054505

Hasil yang diperoleh sering kali adalah minimum lokal, bukan jaminan menghasilkan minimum global. Alasannya, fungsi optimize() bekerja paling baik pada interval yang hanya naik dan turun satu kali. Ketika interval terlalu lebar dan fungsi memiliki banyak osilasi, algoritma dapat “mendarat” pada lembah terdekat yang ditemuinya, bukan pada lembah terdalam. Dari grafik terlihat bahwa lembah terdalam (kandidat minimum global) tampak berada pada kuadran kanan, kira-kira di \(x \in [4, 2\pi]\). Karena itu, kita memangkas interval agar pencarian dilakukan pada segmen yang memuat lembah terdalam tersebut.

min_global <- optimize(f, interval = c(4, 2*pi))
min_global$minimum      # perkiraan x* minimum global

[1] 5.273383

min_global$objective    # nilai minimum global f(x*)

[1] -2.741405

Dengan pemangkasan interval berbasis inspeksi grafik, optimize() bekerja sesuai asumsi dan lebih andal menemukan minimum global dibanding menjalankan pencarian di rentang lebar.

Ketika kita ingin mencari maksimum dari suatu fungsi dengan optimize(), kita cukup menambahkan parameter maximum = TRUE. Dengan begitu, algoritma yang semula mencari titik minimum akan otomatis beralih mencari titik maksimum.

Pertama, kita coba mencari maksimum pada interval penuh \([0, 2\pi]\)

max_02pi <- optimize(f, interval = c(0, 2*pi), maximum = TRUE)
max_02pi$maximum      # x* yang ditemukan

[1] 4.0598

max_02pi$objective    # f(x*) yang ditemukan

[1] 1.096473

Hasil yang diperoleh sering kali merupakan maksimum lokal, alasannya sama seperti pada kasus minimum. Dari plot terlihat bahwa puncak tertinggi berada di sekitar \(x \in [0, 1.5]\). Karena itu, kita batasi interval pencarian pada rentang tersebut

max_global <- optimize(f, interval = c(0, 1.5), maximum = TRUE)
max_global$maximum      # perkiraan x* maksimum global

[1] 0.3323289

max_global$objective    # nilai maksimum global f(x*)

[1] 1.485871

Dengan demikian, pemilihan interval yang sesuai berdasarkan grafik sangat penting agar optimize() dapat menemukan maksimum global, bukan hanya maksimum lokal.

Visualisasi berikut menggambarkan posisi titik minimum dan maksimum yang ditemukan algoritma:

# Plot fungsi
curve(f, from = 0, to = 2*pi, n = 2000,
      main = "Minimum dan Maksimum f(x) = sin(x) + sin(2x) + cos(3x)",
      xlab = "x", ylab = "f(x)")

# Tandai minimum lokal
points(min_lokal$minimum, min_lokal$objective, pch = 19, col = "blue", cex = 1.2)

# Tandai minimum global
points(min_global$minimum, min_global$objective, pch = 19, col = "darkblue", cex = 1.2)

# Tandai maksimum lokal
points(max_02pi$maximum, max_02pi$objective, pch = 17, col = "red", cex = 1.2)

# Tandai maksimum global
points(max_global$maximum, max_global$objective, pch = 17, col = "darkred", cex = 1.2)

# Tambahkan legenda
legend("topright",
       legend = c("Minimum lokal", "Minimum global",
                  "Maksimum lokal", "Maksimum global"),
       col = c("blue", "darkblue", "red", "darkred"),
       pch = c(19, 19, 17, 17), bty = "n")

Dari grafik terlihat jelas perbedaan antara titik ekstrim lokal dan global. Pemangkasan interval secara intuitif berdasarkan grafik awal kita terbukti sangat membantu agar optimize() menemukan titik optimum global, bukan sekadar optimum lokal.

Fungsi optim()

Untuk kasus dengan lebih dari satu variabel, kita bisa menggunakan fungsi optim() di R. Fungsi ini sangat fleksibel karena mendukung berbagai algoritma optimasi. Secara default, optim() menggunakan metode Nelder–Mead, yang bekerja tanpa memerlukan turunan fungsi. Metode ini cocok untuk fungsi non-linear atau fungsi yang sulit diturunkan.

Bentuk umum dari fungsi ini adalah:

optim(par, fn, method = "Nelder-Mead")

dengan:

• par : tebakan awal berupa vektor parameter,
• fn : fungsi yang ingin diminimalkan,
• method : algoritma optimasi yang digunakan (“Nelder-Mead” adalah default).

Sebagai contoh, kita ingin mencari minimum dari fungsi dua variabel:

\[ f(x,y) = (x-2)^2 + (y+1)^2 \]

Fungsi ini memiliki titik minimum global di \((x,y) = (2,-1)\).

# Definisi fungsi dua variabel
fxy <- function(par) {
  x <- par[1]
  y <- par[2]
  (x - 2)^2 + (y + 1)^2
}

# Gunakan optim dengan tebakan awal (0,0)
result_optim <- optim(par = c(0, 0), fn = fxy)
result_optim$par       # titik minimum hampiran (x*, y*)

[1]  1.9998933 -0.9999019

result_optim$value     # nilai minimum f(x*, y*)

[1] 2.101683e-08

Dengan optim(), kita dapat melakukan optimasi untuk fungsi yang lebih kompleks dan berdimensi lebih tinggi dibandingkan optimize(). Dalam contoh ini, hasil optimasi mendekati titik \((2,−1)\) dengan nilai minimum mendekati nol.

Contoh lainnya, misalkan kita ingin mencari nilai \(x_1\) dan \(x_2\) yang meminimalkan fungsi berikut

\[ f(x_1, x_2) = 100(x_2 - x_1^2)^2 + (1 - x_1)^2 \]

Fungsi ini memiliki lembah sempit yang memanjang, sehingga cukup menantang bagi banyak metode numerik. Titik minimum global dari fungsi ini berada pada \((x_1, x_2) = (1,1)\) dengan nilai minimum \(f(x_1,x_2)=0\).

Kita dapat mendefinisikan fungsi ini dalam R dan menggunakan optim() untuk mencari minimum:

# Definisi fungsi
fxx <- function(par) {
  x1 <- par[1]
  x2 <- par[2]
  100 * (x2 - x1^2)^2 + (1 - x1)^2
}

# Gunakan optim dengan tebakan awal (0,0)
resultfxx <- optim(par = c(0, 0), fn = fxx)
resultfxx$par      # titik minimum hampiran (x1*, x2*)

[1] 0.9999564 0.9999085

resultfxx$value    # nilai minimum f(x1*, x2*)

[1] 3.729052e-09

Hasil yang diperoleh mendekati titik \((1,1)\) dengan nilai minimum mendekati nol, sesuai dengan solusi analitik. Contoh ini memperlihatkan bagaimana optim() dapat digunakan untuk fungsi dua variabel yang lebih kompleks dibandingkan kasus kuadrat sederhana.

resultfxx

$par
[1] 0.9999564 0.9999085

$value
[1] 3.729052e-09

$counts
function gradient 
     169       NA 

$convergence
[1] 0

$message
NULL

Dalam penggunaannya, fungsi optim() juga memungkinkan kita memilih berbagai metode optimasi melalui argumen method, di antaranya adalah “Nelder-Mead”, “BFGS”, “CG”, “L-BFGS-B”, “SANN” dan “Brent”.

Ordinary Least Squares

Selain menggunakan fungsi bawaan lm(), estimasi koefisien regresi OLS juga bisa dilakukan dengan pendekatan optimasi numerik.
Prinsipnya, kita meminimalkan jumlah kuadrat galat (Sum of Squared Errors, SSE):

\[ SSE = \sum_{i=1}^n \left(y_i - \hat{y}_i\right)^2 \]

Misalkan kita ingin menduga parameter pada model:

\[ y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \varepsilon \]

Kita buat data simulasi:

set.seed(123)
x1 <- runif(10, 1, 10)
x2 <- runif(10, 1, 10)
galat <- rnorm(10, 0, 0.5)
y <- 1 + 2*x1 + 3*x2 + galat   # model data simulasi

Definisikan fungsi SSE

SSE <- function(par, y, x) {
  X <- cbind(1, x)                     # tambahkan intercept
  yhat <- X %*% as.matrix(par)         # prediksi
  error2 <- sum((y - yhat)^2)          # jumlah kuadrat galat
  return(error2)
}

hasilOLS <- optim(par = c(1, 1, 1), fn = SSE, y = y, x = cbind(x1, x2))
hasilOLS$par   # estimasi koefisien (β0, β1, β2)

[1] 1.856706 1.865194 3.014691

Sebagai pembanding, hasil lm() adalah

lm(y ~ x1 + x2)

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

Coefficients:
(Intercept)           x1           x2  
      1.856        1.866        3.014  

Kedua metode memberikan hasil yang hampir identik, menunjukkan bahwa optim() dapat digunakan untuk mereplikasi estimasi OLS dengan meminimalkan SSE secara numerik.

Maximum Likelihood Estimator

Maximum Likelihood Estimation atau Penduga Peluang Maksimum adalah metode yang sangat penting dalam statistika untuk menduga parameter suatu distribusi probabilitas, berdasarkan data yang diamati. Intuisi dasar dari MLE sangat sederhana namun kuat: MLE mencari parameter yang membuat data yang kita amati paling mungkin terjadi.

Misalkan kita memiliki data hasil pengamatan sebanyak \(n\) nilai:

\[ x_1, x_2, \dots, x_n \]

dan kita mengasumsikan bahwa data tersebut berasal dari distribusi probabilitas \(f(x \mid \theta)\), tetapi parameter \(\theta\) belum diketahui. Tugas kita adalah menduga nilai \(\theta\) berdasarkan data tersebut.

Fungsi likelihood didefinisikan sebagai:

\[ L(\theta \mid x_1, x_2, \dots, x_n) = \prod_{i=1}^{n} f(x_i \mid \theta) \]

Nilai \(\hat{\theta}\) yang memaksimumkan fungsi likelihood disebut penduga MLE:

\[ \hat{\theta}_{\text{MLE}} = \arg\max_\theta L(\theta \mid x_1, \dots, x_n) \]

Karena bentuk produk sering menyulitkan perhitungan, kita menggunakan bentuk log-likelihood:

\[ \ell(\theta) = \log L(\theta) = \sum_{i=1}^{n} \log f(x_i \mid \theta) \]

Memaksimumkan \(\ell(\theta)\) setara dengan memaksimumkan \(L(\theta)\).

Kita dapat mencari MLE dengan cara

1. Secara analitik: jika turunan log-likelihood dapat dihitung dan diselesaikan secara aljabar,

   \[ \frac{d\ell(\theta)}{d\theta} = 0 \]

   dan solusi \(\theta\) pada titik stasioner menjadi kandidat MLE.

2. Secara numerik: jika fungsi log-likelihood terlalu rumit untuk diselesaikan secara manual, kita gunakan metode optimasi numerik untuk meminimalkan:

   \[ -\ell(\theta) \]

Sebagai contoh, misalkan \(x_1, x_2, \dots, x_n\) merupakan sampel acak dari peubah acak \(X \sim \mathcal{N}(\mu, \sigma)\). Kita ingin menduga parameter \(\mu\) dan \(\sigma\) dari data tersebut menggunakan pendekatan MLE.

PDF dari distribusi Normal adalah:

\[ f(x \mid \mu, \sigma) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left(-\frac{(x - \mu)^2}{2\sigma^2}\right) \]

Fungsi log-likelihood untuk data \(x_1, \dots, x_n\) menjadi:

\[ \ell(\mu, \sigma) = \sum_{i=1}^n \log f(x_i \mid \mu, \sigma) = -\frac{n}{2} \log(2\pi) - n \log \sigma - \frac{1}{2\sigma^2} \sum_{i=1}^n (x_i - \mu)^2 \]

Untuk mencari nilai \(\mu\) dan \(\sigma\) yang memaksimalkan \(\ell(\mu, \sigma)\), kita gunakan pendekatan numerik. Karena fungsi optim() di R mencari minimum, maka kita minimalkan bentuk negatif dari log-likelihood:

# Fungsi log-likelihood negatif
negloglik <- function(para, data) {
  mu <- para[1]
  sigma <- para[2]
  -sum(dnorm(data, mean = mu, sd = sigma, log = TRUE))
}

Simulasi data dari distribusi Normal \(\mathcal{N}(2, 25)\):

set.seed(123)
x <- rnorm(100, mean = 2, sd = sqrt(25))

hasil <- optim(par = c(1, 1), fn = negloglik, data = x)
hasil$par  # Estimasi (μ, σ) dengan MLE

[1] 2.453112 4.540807

Sebagai pembanding, estimasi dengan metode fungsi mean dan standar deviasi dari data

c(mean(x), sd(x))

[1] 2.452030 4.564079

Misalkan kita ingin melihat bagaimana bentuk log-likelihood berubah terhadap \(\mu\) jika \(\sigma\) dianggap tetap:

mu_grid <- seq(0, 4, length.out = 100)
loglik_values <- sapply(mu_grid, function(mu) -negloglik(c(mu, sd(x)), x))

plot(mu_grid, -loglik_values, type = "l",
     main = "Log-Likelihood terhadap μ (σ tetap)",
     xlab = expression(mu), ylab = "Log-Likelihood")
abline(v = hasil$par[1], col = "red", lty = 2)

Dalam regresi linear klasik, kita memodelkan hubungan antara variabel respon \(y\) dan satu atau lebih variabel prediktor melalui persamaan linear:

\[ y_i = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \cdots + \beta_p x_{ip} + \varepsilon_i \]

untuk \(i = 1, 2, \dots, n\), dengan asumsi bahwa galat \(\varepsilon_i\) bersifat independen dan identik terdistribusi Normal, yaitu \(\varepsilon_i \sim \mathcal{N}(0, \sigma^2)\). Berdasarkan asumsi ini, nilai \(y_i\) juga mengikuti distribusi Normal:

\[ y_i \sim \mathcal{N}(\mu_i, \sigma^2) \quad \text{dengan} \quad \mu_i = \boldsymbol{x}_i^\top \boldsymbol{\beta} \]

di mana \(\boldsymbol{x}_i\) adalah vektor prediktor untuk observasi ke-\(i\), termasuk elemen 1 sebagai intercept. Dengan demikian, kita dapat membentuk fungsi likelihood berdasarkan distribusi Normal dari seluruh observasi:

\[ L(\boldsymbol{\beta}, \sigma^2 \mid \boldsymbol{y}) = \prod_{i=1}^n \frac{1}{\sqrt{2\pi \sigma^2}} \exp\left(-\frac{(y_i - \boldsymbol{x}_i^\top \boldsymbol{\beta})^2}{2\sigma^2}\right) \]

Untuk mempermudah perhitungan, fungsi likelihood ini biasanya diubah menjadi bentuk log-likelihood:

\[ \ell(\boldsymbol{\beta}, \sigma^2) = \log L(\boldsymbol{\beta}, \sigma^2 \mid \boldsymbol{y}) = - \frac{n}{2} \log(2\pi) - \frac{n}{2} \log(\sigma^2) - \frac{1}{2\sigma^2} \sum_{i=1}^n (y_i - \boldsymbol{x}_i^\top \boldsymbol{\beta})^2 \]

Tugas utama dalam estimasi maksimum likelihood adalah mencari nilai parameter \(\boldsymbol{\beta}\) dan \(\sigma^2\) yang memaksimumkan fungsi log-likelihood tersebut. Karena bentuk log-likelihood bergantung pada jumlah kuadrat galat, maka memaksimalkannya ekuivalen dengan meminimalkan:

\[ S(\boldsymbol{\beta}) = \sum_{i=1}^n (y_i - \boldsymbol{x}_i^\top \boldsymbol{\beta})^2 \]

yang merupakan fungsi kuadrat sisaan total (Sum of Squared Errors/SSE). Proses ini identik dengan metode kuadrat terkecil (Ordinary Least Squares/OLS). Maka, estimasi MLE untuk parameter regresi diperoleh dengan menyelesaikan sistem persamaan normal berikut:

\[ \hat{\boldsymbol{\beta}} = (X^\top X)^{-1} X^\top \boldsymbol{y} \]

dan penduga untuk variansi galat diberikan oleh:

\[ \hat{\sigma}^2 = \frac{1}{n} \sum_{i=1}^n (y_i - \boldsymbol{x}_i^\top \hat{\boldsymbol{\beta}})^2 \]

Perlu dicatat bahwa pembagi \(n\) pada \(\hat{\sigma}^2\) adalah khas pendekatan likelihood, sedangkan dalam OLS biasanya digunakan pembagi \(n - p\) untuk memperoleh penduga yang tak bias terhadap \(\sigma^2\). Dengan demikian, pada regresi linear dengan asumsi normalitas galat, metode Maximum Likelihood menghasilkan penduga \(\boldsymbol{\beta}\) yang sama dengan OLS, namun memberikan perbedaan pada pendugaan varians.

Berikut ditunjukkan dua cara implementasi pendugaan parameter regresi linear menggunakan pendekatan Maximum Likelihood dalam R.

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

# Matriks desain (X) dengan intercept
X <- cbind(1, x)
n <- nrow(X)

# Pendekatan Pertama: MLE manual via log-likelihood fungsi normal
regloglik_1 <- function(beta, y, X) {
  # beta[1] = intercept, beta[2] = slope, beta[3] = varians galat (σ^2)
  e <- y - X %*% beta[1:2]
  var_lin <- beta[3]
  
  loglik <- -0.5 * n * log(2 * pi) -
            0.5 * n * log(var_lin) -
            (t(e) %*% e) / (2 * var_lin)
  
  return(-loglik)  # fungsi optim akan meminimalkan
}

# Estimasi parameter MLE
mle_1 <- optim(par = c(1, 1, 1), fn = regloglik_1, y = y, X = X)
mle_1$par  # [Intercept, Slope, σ²]

[1] -1.2669724  2.0286218  0.4698408

# Pendekatan Kedua: Langsung gunakan dnorm() dengan log-likelihood
regloglik_2 <- function(par, y, x) {
  mu <- par[1] + par[2] * x
  sigma <- par[3]
  
  loglik <- sum(dnorm(y, mean = mu, sd = sigma, log = TRUE))
  return(-loglik)  # Negatif karena ingin memaksimalkan
}

# Estimasi parameter MLE
mle_2 <- optim(par = c(1, 1, 1), fn = regloglik_2, y = y, x = x)

# Hasil MLE (Pendekatan 2)
mle_2$par  # tampilkan parameter (β0, β1, σ)

[1] -1.2664693  2.0285134  0.6853899

Kedua pendekatan tersebut menghasilkan estimasi parameter regresi linear menggunakan prinsip Maximum Likelihood, namun dibangun dengan cara yang berbeda. Pendekatan pertama menyusun log-likelihood secara manual berdasarkan bentuk eksplisit fungsi densitas Normal, di mana parameter ketiga diasumsikan sebagai varians galat (\(\sigma^2\)). Sementara itu, pendekatan kedua memanfaatkan langsung fungsi dnorm() dengan parameter simpangan baku (\(\sigma\)), sehingga formulanya lebih ringkas dan lebih dekat dengan pemanfaatan fungsi probabilitas bawaan dalam R. Meskipun struktur fungsi berbeda, kedua pendekatan tersebut akan memberikan hasil estimasi parameter yang sangat mirip, karena didasarkan pada asumsi distribusi galat yang sama. Pilihan antara keduanya bergantung pada preferensi formulasi: eksplisit melalui bentuk matematis penuh, atau implisit menggunakan fungsi probabilitas.

Muhammad Yusran
Program Studi Magister Statistika dan Sains Data
Sekolah Sains Data, Matematika, dan Informatika
IPB University

Email: muhammadyusran@apps.ipb.ac.id
LinkedIn: Muhammad Yusran
Instagram: mhammad.yusran_