Tugas Praktikum Ketigas STA561 Pemrograman Statistika
Klik disini untuk ke halaman rpubs.
Fungsi dan OOP (Praktikum Pertemuan Keenam)
Fungsi
R telah menyiapkan banyak fungsi yang dapat dimanfaatkan untuk komputasi dan grafik. Fungsi yang tidak ada dalam R karena berkembangnya statistika dapat diciptakan sendiri.
Membuat sekumpulan mekanisme dasar yang dijalankan secara simultan secara umum digunakan Syntax:
function( arglist ) { expr return(value) }
#Return value
angka_acak1=function(n,pw)
{ x=runif(n)
y=runif(n)
z=(x+y)^pw
return(z)
}
angka_acak1(10,2)## [1] 0.37910386 0.65399185 0.61955430 0.09504381 2.08552152 1.69692971
## [7] 0.33464538 1.40952429 0.19876157 0.29361407#Return value
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.06925621 0.21146990 0.79076226 0.75413825 0.30374451 0.71856188
## [7] 0.39196244 0.37674184 0.18063171 0.75784270
##
## $y
## [1] 0.6834656 0.9850363 0.9601427 0.1307271 0.8079957 0.3476402 0.8804274
## [8] 0.2528720 0.6430956 0.6706402
##
## $z
## [1] 0.5665901 1.4316270 3.0656682 0.7829868 1.2359662 1.1367868 1.6189758
## [8] 0.3964136 0.6785267 2.0405635#Return value
angka_acak3=function(n=10,pw=2)
{ x=runif(n)
y=runif(n)
z=(x+y)^pw
return(z)
}
angka_acak3()## [1] 0.9528312 1.3862525 0.8747391 1.2777245 2.0903869 0.6351781 0.5438862
## [8] 0.5139846 0.9744523 0.2760594#Return value
angka_acak4=function()
{ x=runif(n)
y=runif(n)
z=(x+y)^pw
return(z)
}
n <- 5; pw <- 3
angka_acak4()## [1] 0.67211248 1.04457223 1.76399170 0.03836678 1.33339392Latihan 1
Buatlah fungsi untuk mencari median dari suatu vektor!
#mencari median
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.5Latihan 2
Buatlah fungsi untuk mencari modus dari suatu vektor!
#Mencari modus
modus <- function(vect)
{ v <- unique(vect)
f <- NULL
for(i in v)
{ byk <- sum(vect==i)
f <- c(f,byk)
}
fmax <- max(f)
vf <- cbind(v,f)
mode <- vf[f==fmax,]
return(mode)
}
modus(x1)## v f
## [1,] 3 2
## [2,] 7 2Latihan 3
Buatlah fungsi untuk menduga parameter pada regresi berganda!
#Membuat penduga parameter regresi berganda
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)
}
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.01368033model<-lm(Pendapatan~Biaya.Iklan+Jumlah.Warung)
model$coefficients## (Intercept) Biaya.Iklan Jumlah.Warung
## -0.21381852 0.89843390 0.01745279model$fitted.values## 1 2 3 4 5 6 7
## 3.094910 3.277176 2.830539 3.184754 3.988185 2.738116 2.286320model$residuals## 1 2 3 4 5 6
## 0.40508982 -0.07717642 0.16946108 -0.28475357 0.01181483 -0.23811608
## 7
## 0.01368033Latihan 4
Buat fungsi bernama three.M yang digunakan untuk menghitung mean, median, dan modus dari suatu vektor modus (tanpa menggunakan fungsi mean, quantile, ataupun fungsi “instan” lain yang sudah tersedia sebelumnya di R)!
Hitung mean, median, modus dari X dibawah ini dengan menggunakan fungsi three.M tersebut!
set.seed(123) X <- rbinom(100,10,0.5)
#Three.M
three.M <- function(vect) {
n <- length(vect) # banyak data
# menghitung rataan
jumlah <- sum(vect)
rataan <- jumlah/n
# menhitung median
vects <- sort(vect) # urutkan
if(n%%2 == 1) {m <- vects[(n+1)/2]}
else {m <- (vects[n/2]+vects[(n/2)+1])/2}
# menghitung modus
v <- unique(vect)
f <- NULL
for(i in v)
{ byk <- sum(vect==i)
f <- c(f,byk)
}
fmax <- max(f)
vf <- cbind(v,f)
mode <- vf[f==fmax,]
# output
my_list <- list("mean" = rataan, "median" = m, "modus" = mode)
return(my_list)
}
# hitung mean median modus distribusi binomial
set.seed(123)
x1 <- rbinom(100,10,0.5)
three.M(x1)## $mean
## [1] 4.99
##
## $median
## [1] 5
##
## $modus
## v f
## [1,] 6 24
## [2,] 5 24Object Oriented Programming (OOP)
Object Oriented Programming:
A class defines the behavior of objects by describing their attributes and their relationship to other classes.
A method is functions that behave differently depending on the class of their input.
Classes are usually organized in a hierarchy: if a method does not exist for a child, then the parent’s method is used instead; the child inherits behavior from the parent.
Pengembangan awal objek di R menggunakan Class System S3 yang tidak terlalu ketat.
Pendefinisian yang ketat secara formal, R menggunakan Class System S4
Misal:
Objeknya: mobil
kelasnya: kendaraan roda empat.
Metodenya: bisa rem, mengisi bensin, mengatur kecepatan.
Attribut: bahan bakar apa, dimensi kendaraan, komponen apa saja yang ada.
Konsep OOP dalam R:
• Encapsulation
• Inheritance
numeric(), matrix(), list()
• Polymorphism
plot.acf(…), plot.ecdf(…), plot.data.frame(…), plot.lm(…)
Object-S3
A1 <- c(1:10)
class(A1)## [1] "integer"A2 <- matrix(A1,2,5)
class(A2)## [1] "matrix" "array"A3 <- 1:12
A4 <- letters[1:12]
B1 <- data.frame(A3,A4)
class(B1)## [1] "data.frame"B1$A4## [1] "a" "b" "c" "d" "e" "f" "g" "h" "i" "j" "k" "l"A5 <- 10+A3+rnorm(12)
B2 <- lm(A5~A3) #membuat model linear
class(B2)## [1] "lm"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 codesummary(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-08names(B2)## [1] "coefficients" "residuals" "effects" "rank"
## [5] "fitted.values" "assign" "qr" "df.residual"
## [9] "xlevels" "call" "terms" "model"B2$coefficients## (Intercept) A3
## 10.3946172 0.9661381Mengubah menjadi Class:
class(obj) <- “class.name”
Mobil1 <- list(Nama="Toyota", Panjang=3.5,
Lebar=2, Kecepatan=180)
class(Mobil1)## [1] "list"class(Mobil1) <- "mobil"class(Mobil1) <- "mobil"
Mobil2 <- list(Nama="Suzuki", Panjang=1,
Lebar=1.8, Kecepatan=150)
class(Mobil2) <- "mobil"
class(Mobil2)## [1] "mobil"Mengubah menjadi Class (Menggunakan Fungsi Konstruktor):
Lebih direkomendasikan menggunakan fungsi konstruktor
Fungsi konstruktor menambahkan screening sebelum menambahkan class
#contoh input data yang sesuai kriteria
Mobil <- function(Nama,Panjang,Lebar,Kecepatan){
if(Panjang<2 || Lebar<1.5 || Kecepatan<80)
stop("atribut tidak sesuai")
Mobil <- list(Nama=Nama, Panjang=Panjang,
Lebar=Lebar, Kecepatan=Kecepatan)
class(Mobil) <- "mobil"
Mobil
}
Mobil3 <- Mobil("Daihatsu", 2.1, 1.9, 120)
Mobil3## $Nama
## [1] "Daihatsu"
##
## $Panjang
## [1] 2.1
##
## $Lebar
## [1] 1.9
##
## $Kecepatan
## [1] 120
##
## attr(,"class")
## [1] "mobil"#contoh input data yang tidak sesuai kriteria
#Mobil <- function(Nama,Panjang,Lebar,Kecepatan){
#if(Panjang<2 || Lebar<1.5 || Kecepatan<80)
#stop("atribut tidak sesuai")
#Mobil <- list(Nama=Nama, Panjang=Panjang,
#Lebar=Lebar, Kecepatan=Kecepatan)
#class(Mobil) <- "mobil"
#Mobil
#}
#Mobil4 <- Mobil("Proton", 2, 1.8, 70)
#Mobil4Fungsi Aksesor
Cara langsung (tidak direkomendasikan):
Mobil2$Nama
Mobil3$Panjang
Mobil2$Nama## [1] "Suzuki"Mobil3$Panjang## [1] 2.1Dengan fungsi aksesor: nama <- function(objek) objek$Nama
kecepatan <- function(objek) objek$Kecepatan
nama(Mobil1)
kecepatan(Mobil3)
nama <- function(objek) objek$Nama
kecepatan <- function(objek) objek$Kecepatan
nama(Mobil1)## [1] "Toyota"kecepatan(Mobil3)## [1] 120Fungsi Generik Print Method
function(…){…}
print.mobil <- function(objek) {
print(cat("Nama : ", nama(objek), "\n",
"Kecepatan : ", kecepatan(objek),
sep="")
)
}
Mobil1## Nama : Toyota
## Kecepatan : 180NULLMenciptakan Fungsi Generik
Method hanya dapat didefinisikan untuk fungsi yang generik. Membuat nama method baru dengan menciptakan fungsi generik
fungsibaru <- function (objek) UseMethod("fungsibaru")
Object-S4
• Mengatasi masalah dalam sistem objek S3 dengan sistem objek lebih formal
• Mendefinisikan class baru :
setClass(“class.name”, representation(x=“type”), prototype(x=“…”))
• Mendefinisikan objek baru dari kelas tertentu :
new(class.name,…)
#Membuat Object S4
setClass("car",
representation(Nama="character",
Panjang="numeric",
Lebar="numeric",
Kecepatan="numeric"))
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#Membuat Object S4 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)
class(Car2)## [1] "car"
## attr(,"package")
## [1] ".GlobalEnv"class(Mobil1)## [1] "mobil"Car2## An object of class "car"
## Slot "Nama":
## [1] "Suzuki"
##
## Slot "Panjang":
## [1] 2.4
##
## Slot "Lebar":
## [1] 1.8
##
## Slot "Kecepatan":
## [1] 150Akses terhadap Slot
Cara langsung (tidak direkomendasikan)
# Akses terhadap slot dengan Cara langsung
Car1@Nama## [1] "Toyota"Car2@Kecepatan## [1] 150- Dengan fungsi aksesor:
# Akses terhadap slot dengan fungsi aksesor
nama1 <- function(objek) objek@Nama
kecepatan1 <- function(objek) objek@Kecepatan
nama1(Car1)## [1] "Toyota"kecepatan1(Car2)## [1] 150Class Method
setMethod(“method”,“class.name”, function(…){…})
# Class Method
setMethod(show, "car", function(object) {
print(cat("Nama : ", nama1(object), "\n",
"Kecepatan : ", kecepatan1(object),
sep="")
)}
)
Car2## Nama : Suzuki
## Kecepatan : 150NULLMenciptakan Fungsi Generik S4:
setGeneric(“fungsibaru”, function(objek) standardGeneric(“fungsibaru”))
setGeneric("fungsibaru",
function(objek)
standardGeneric("fungsibaru"))## [1] "fungsibaru"Optimasi (Praktikum Pertemuan Ketujuh)
#Operasi yang kurang tepat bilangan riil dalam basis 2
(.3 - .1)== .2## [1] FALSE.2 - (.3 - .1) #jika pangkat dari sudah mencapai -6 maka data sudah dianggap 0 ## [1] 2.775558e-17#Mengatasi Eror basis 2
isTRUE (all.equal(.2 , .3 - .1) )## [1] TRUE#Mengatasi Eror basis 2
all.equal(.2 , .3) # not a logical value## [1] "Mean relative difference: 0.5"#Mengatasi Eror basis 2
isTRUE (all.equal(.2 , .3) ) # always a logical value## [1] FALSEDiferensial
#Diferensial
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,] 4Integral
#Integral
fs<-function(x) x ^2
integrate (fs ,0 ,1)## 0.3333333 with absolute error < 3.7e-15Latihan 1
#Latihan 1
#Integral
#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.triyac_str("Integrate(x) x^2 + 4*x")## [1] "x^3/3+2*x^2"f1 <-function(x) x^2 + 4*x
integrate(f1,lower = -10,upper = 10)## 666.6667 with absolute error < 7.6e-12Latihan 2
#Latihan 2
#Integral
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)"f2 <-function(t) t^(4) * exp(-t)
integrate(f2,lower = 0,upper = Inf)## 24 with absolute error < 2.2e-05gamma(5)## [1] 24Golden Section Search
#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)
}
f<-function(x) {
+ abs(x -3.5) + (x -2) ^2
}
curve(f, from = 1, to = 5)
golden (f ,1 ,2)## [1] 2golden (f ,1 ,5)## [1] 2.5golden (f ,3 ,5)## [1] 3Newton Raphson
#Newton Raphson
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 (x - x0 )
x <- x0 - f1/f2
x0 <- x
}
return( x )
}
fx<-expression(4*x ^2 - 3*x - 7)
newtonr (fx ,3)## [1] 0.375fx<-expression(exp(-x) + x^4)
newtonr ( fx )## [1] 0.5282519fx<-expression(x^2 - x)
newtonr ( fx )## [1] 0.5Optimasi Satu Dimensi (Optimize)
Digunakan untuk mendapatkan nilai minimum dari suatu fungsi dengan satu peubah.
Misal akan dicari nilai minimum dari fungsi: f(x) = (x - 1/3)^2
#Optimize/Optimise
f<-function(x , a) (x - a)^2
xmin<- optimize(f ,c(0 , 1) ,
tol = 0.0001 , a = 1/3)
xmin## $minimum
## [1] 0.3333333
##
## $objective
## [1] 0fx <- function(x)(x-(1/3))^2
curve(fx,from=-1.5, to=2) 
Latihan 3
#Latihan 3
#Optimasi Satu Dimensi
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.054505optimize(f3, interval = c(4, 2*pi)) #minimum global## $minimum
## [1] 5.273383
##
## $objective
## [1] -2.741405optimize(f3, interval = c(0, 2*pi), maximum = T) #maksimum lokal ## $maximum
## [1] 4.0598
##
## $objective
## [1] 1.096473#default maximum = F untuk mencari titik minimum
optimize(f3, interval = c(0, 1.5), maximum = T) #maksimum global## $maximum
## [1] 0.3323289
##
## $objective
## [1] 1.485871Optimasi Fungsi Polinomial (Optim)
Untuk mencari nilai minimum dari fungsi lebih dari satu peubah.
Misal mencari nilai x1 dan x2 yang membuat f(x1,x2) = 100(x2- x1^2 ) ^2 + (1-x1) ^2 minimum.
#Optim
fr <- function(x) {
x1 <- x[1]
x2 <- x[2]
+ 100*(x2 - x1^2)^2
+ (1 - x1)^2
}
optim(c( -1.2 ,1) , fr )## $par
## [1] 1.0002837 -0.1642825
##
## $value
## [1] 8.050091e-08
##
## $counts
## function gradient
## 47 NA
##
## $convergence
## [1] 0
##
## $message
## NULLLatihan 4
#Latihan 4
#Optimasi Fungsi Polinomial
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)## Warning in optim(par = c(-0.5), fn = f4): one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly## $par
## [1] 0.728418
##
## $value
## [1] -1.832126
##
## $counts
## function gradient
## 36 NA
##
## $convergence
## [1] 0
##
## $message
## NULL#Metode Kuadrat Terkecil dengan Optim
f<-function( para , y , x){
X<- cbind(1 , x)
yhat<- X %*%as.matrix( para )
sisa2<- sum((y - yhat ) ^2)
return( sisa2 )
}
x1<-runif(10 ,1 ,10)
x2<-runif(10 ,1 ,10)
galat<-rnorm(10 ,0 ,0.5)
y<-1 + 2*x1 + 3*x2 + galat
hasil <- optim(c(1 ,1 ,1) ,f ,y=y ,x=cbind(x1,x2))
hasil$par## [1] 1.881706 1.914983 2.927032lm(y~x1 + x2 )##
## Call:
## lm(formula = y ~ x1 + x2)
##
## Coefficients:
## (Intercept) x1 x2
## 1.875 1.915 2.928Regresi Linier
Latihan 5
#Latihan 5
#Regresi Linier
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.028449hasil2$coefficients## (Intercept) x
## -1.266667 2.028571hasil1$value## [1] 2.819048sum(hasil2$residuals^2)## [1] 2.819048(.3 - .1)## [1] 0.2MLE: Maximum Likelihood Estimator
Metode ini merupakan metode yang paling sering digunakan untuk menduga parameter sebaran
Teladan: Jika x1,x2,… xn berasal dari peubah acak X~N (miyu,sigma), tentukan penduga bagi miyu dan sigma menggunakan MLE.
#MLE
negloglik<-function( para , xd ){
nilai<--1*sum(dnorm(xd ,mean= para [1] ,sd= para [2] ,log=
TRUE ))
+ return( nilai )
}
x<-rnorm(10 ,2 ,5)
hasil<-optim(c(1 ,1) , negloglik , xd =x)
hasil$par## [1] 4.141031 2.730255c(mean(x) ,sd(x)) # pembanding## [1] 4.140905 2.879107Tugas
Tentukan pendugaan parameter dalam regresi linier dengan menggunakan kemungkinan maksimum (MLE)!
#Regresi Linier dengan MLE
regloglik<-function( para, y , x ){
nilai<--1*sum(dnorm(y ,mean= para [1] + x*para [2] ,sd= para [3] ,log=
TRUE ))
return( nilai )
}
x<-c(2,6,3,8,13)
y<-c(8,17,15,11,19)
hasil<-optim(c(1 ,1,1) , regloglik , y=y, x=x)
hasil$par## [1] 9.8544869 0.6476157 3.0848436lm(y~x) # pembanding##
## Call:
## lm(formula = y ~ x)
##
## Coefficients:
## (Intercept) x
## 9.8549 0.6477Diperoleh nilai penduga B0 dan B1 yang sama antara metode langsung menggunakan lm() dan metode MLE.
Eksplorasi Mandiri (Pendugaan Regresi Linier dengan Menggunakan MLE)
Menentukan pendugaan regresi linier dengan menggunakan kemungkinan maksimum (MLE).
Model Regresi Linier Sederhana: Yi = B0 + B1Xi + Ei ; i = 1,2,…,n
Membangkitkan Bilangan Acak
#Membangkitkan Bilangan Acak X dari distribusi normal
data.x <- rnorm(n = 100, mean = 10, sd = 2)#Membangkitkan nilai Y dengan menentukan nilai B0 adalah 8 dan B1 adalah 3
B0 <- 8
B1 <- 3
Y <- B0 + data.x*B1
#Lalu menentukan noise atau error dari data dengan sd = 1
noise.err.sd. <- 1
noise_error <- rnorm(100,0,1)
data.y <- Y + noise_error
plot(data.x, data.y)
Terlihat bahwa data X dan Y yang dibangkitkan memiliki hubungan linier.
Pendugaan Regresi Linier dengan MLE
Diasumsikan berdistribusi Normal (miyu = B0 + B1X1, sigma=standar deviasi). Akan dicari nilai penduga beta nol dan beta satu. Selain itu, akan dicari juga nilai standar deviasi.
#MLE
lm.loss <- function(par) {
B0.par <- par[1] # parameter B0 atau intercept
B1.par <- par[2] # parameter B1
err.sigma <- par[3] #standard deviasi
if(err.sigma < 0) {deviance <- 10000000}
# prosedur optimalisasi: dibuat nilai toleransi jika standar deviasi tidak valid misal negatif, maka akan dikembalikan kepada deviance yang sangat tinggi.
if(err.sigma > 0) {
likelihoods <- dnorm(data.y, mean = B0.par + data.x *B1.par, sd = err.sigma) #Spesifikasi Model
log.likelihoods <- log(likelihoods) # Mengkalkulasi nilai likelihood untuk setiap data
deviance <- -2 * sum(log.likelihoods) # Menyusun deviance (-2*sum of log-likelihoods)
}
return(deviance)
}
parameter.fits <- optim(par = c(1, 1, 10),
fn = lm.loss, hessian = T
)
parameter.fits## $par
## [1] 8.0168814 3.0034131 0.9391614
##
## $value
## [1] 271.2381
##
## $counts
## function gradient
## 188 NA
##
## $convergence
## [1] 0
##
## $message
## NULL
##
## $hessian
## [,1] [,2] [,3]
## [1,] 226.75114386 2.243185e+03 -0.01707136
## [2,] 2243.18529325 2.304490e+04 0.04962114
## [3,] -0.01707136 4.962114e-02 453.53369504Nilai $par pada output di atas secara berturut-turut menunjukkan nilai dari penduga B0, B1, dan standar deviasi.
Membandingkan dengan fungsi lm()
Membandingkan hasil optimasi MLE sebelumnya dengan fungsi regresi linier biasa dengan menggunakan lm()
built.in <- lm(data.y ~ data.x)
built.in##
## Call:
## lm(formula = data.y ~ data.x)
##
## Coefficients:
## (Intercept) data.x
## 8.016 3.004Diperoleh nilai penduga B0 dan B1 yang sama dengan metode optimasi MLE yang telah dilakukan sebelumnya.
TERIMA KASIH
Tugas Praktikum 3 STA561 Pemrograman Statistika, Reni Amelia, Mahasiswa Pascasarjana Statistika dan Sains Data, IPB University, reniamelia@apps.ipb.ac.id↩︎