Laporan 3 Praktikum Pemrograman Statistik - Tugas Satria June Adwendi

Praktikum Pertemuan 6 - Fungsi & OOP

Fungsi

Membuat sekumpulan mekanisme dasar yang dijalankan secara simultan

function( arglist) { expr return(value) }

Return

angka_acak1=function(n,pw)
{x=runif(n)
y=runif(n)
z=(x+y)^pw
return(z)
}
angka_acak1(10,2)

##  [1] 0.1840032 0.5650315 2.1324597 1.0922413 1.0634348 0.3607675 2.5355532
##  [8] 2.6799679 1.2458474 0.3034951

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.34771411 0.28272275 0.98891855 0.97539718 0.01094555 0.97042143
##  [7] 0.04266487 0.53694026 0.66027109 0.01821078
## 
## $y
##  [1] 0.2770898 0.1794267 0.3863038 0.6795290 0.7899223 0.6080055 0.2755837
##  [8] 0.9441639 0.7249228 0.8457920
## 
## $z
##  [1] 0.3903799 0.2135821 1.8912364 2.7387806 0.6413893 2.4914317 0.1012821
##  [8] 2.1936695 1.9187621 0.7465008

angka_acak3=function(n=10,pw=2)
{x=runif(n)
y=runif(n)
z=(x+y)^pw
return(z)
}
angka_acak3()

##  [1] 0.5076296 3.1583038 1.1074516 0.2111496 0.2117834 0.1122646 0.1005030
##  [8] 0.4399891 1.6763004 2.0489449

angka_acak4=function()
{x=runif(n)
y=runif(n)
z=(x+y)^pw
return(z)
}
n <-5; pw <-3
angka_acak4()

## [1] 0.03909105 0.85385876 0.11642387 0.05214659 1.93985231

Latihan 1

Buatlah fungsi untuk mencari median dari suatu vektor

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.5

Latihan 2

Buatlah fungsi untuk mencari modus dari suatu vektor

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 2

Latihan 3

Buatlah fungsi untuk menduga parameter pada regresi berganda!

p.est<-function(A){
if(!is.matrix(A))
stop("inputmustbeon 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.01368033

model<-lm(Pendapatan~Biaya.Iklan+Jumlah.Warung)
model$coefficients

##   (Intercept)   Biaya.Iklan Jumlah.Warung 
##   -0.21381852    0.89843390    0.01745279

model$fitted.values

##        1        2        3        4        5        6        7 
## 3.094910 3.277176 2.830539 3.184754 3.988185 2.738116 2.286320

model$residuals

##           1           2           3           4           5           6 
##  0.40508982 -0.07717642  0.16946108 -0.28475357  0.01181483 -0.23811608 
##           7 
##  0.01368033

Latihan 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 <- function(vect) {
n <- length(vect)

# menghitung rataan 
jumlah <- sum(vect)
rataan <- jumlah/n

# menghitung 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 24

Object Oriented Programming (OOP)

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.

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) #membuatmodel 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 code

summary(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-08

names(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.9661381

Mengubah menjadi Class

class(obj) <-“class.name”

Mobil1 <-list(Nama="Toyota", Panjang=3.5,
Lebar=2, Kecepatan=180)
class(Mobil1)

## [1] "list"

class(Mobil1) <-"mobil"
Mobil2 <-list(Nama="Suzuki", Panjang=1,
Lebar=1.8, Kecepatan=150)
class(Mobil2) <-"mobil"

Mengubah menjadi Class (Menggunakan Fungsi Konstruktor)

•Lebih direkomendasikan menggunakan fungsi konstruktor •Fungsi konstruktor menambahkan screening sebelum menambahkan class

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)
#Mobil4 <-Mobil("Proton", 2, 1.8, 70) , jika dimasukan ini akan erorr karena atribut tidak sesuai

Funsi Aksesor

•Cara langsung (tidak direkomendasikan)

Mobil2$Nama

## [1] "Suzuki"

Mobil3$Panjang

## [1] 2.1

•Dengan fungsi aksesor

nama <-function(objek) objek$Nama
kecepatan <-function(objek) objek$Kecepatan

nama(Mobil1)

## [1] "Toyota"

kecepatan(Mobil3)

## [1] 120

Fungsi Generik Print Method

function(…){…}

print.mobil <-function(objek) {
print(cat("Nama : ",nama(objek),"\n",
"Kecepatan : ",kecepatan(objek),
sep="")
)
}
Mobil1

## Nama : Toyota
## Kecepatan : 180NULL

Menciptakan Fungsi Generik

•Method hanya dapat didefinisikan untuk fungsi yang generik

•Membuat nama method baru dengan menciptakan fungsi generik >fungsibaru <-function (objek) UseMethod(“fungsibaru”)

fungsibaru <-function (objek) UseMethod("fungsibaru")
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=“…”))

•Medefinisikan objek baru dari kelas tertentu: >new(class.name,…)

setClass("car",
        representation(Nama="character",
        Panjang="numeric",
        Lebar="numeric",
        Kecepatan="numeric"))

Car1 <-new("car", Nama="Toyota",
            Panjang=3.5, Lebar=2,
            Kecepatan=180)

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"

Akses terhadap slot

•Cara langsung (tidak direkomendasikan)

Car1@Nama

## [1] "Toyota"

Car2@Kecepatan

## [1] 150

•Dengan fungsi aksesor

nama1 <-function(objek) objek@Nama
kecepatan1 <-function(objek) objek@Kecepatan
nama1(Car1)

## [1] "Toyota"

kecepatan1(Car2)

## [1] 150

Class Method

setMethod(“method”,“class.name”,function(…){…})

setMethod(show, "car", function(object) {
print(cat("Nama: ", nama1(object), "\n",
"Kecepatan: ", kecepatan1(object),
sep="")
)}
)
Car2

## Nama: Suzuki
## Kecepatan: 150NULL

MenciptakanFungsiGenerikS4

setGeneric("fungsibaru",
function(objek)
standardGeneric("fungsibaru"))

## [1] "fungsibaru"

Praktikum Pertemuan 7 - OPTIMASI

Latihan 1 (Integral)

library(Ryacas)

## Warning: package 'Ryacas' was built under R version 4.0.4

## 
## 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

yac_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-12

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-05

gamma(5)

## [1] 24

Latihan 3 (Optimasi Satu Dimensi)

Carilah titik maksimum dan minimum dari fungsi berikut >𝑓(𝑥)=sin(𝑥)+ sin(2𝑥)+ cos(3𝑥)!

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

Latihan 4 (Optimasi Fungsi Polinomial)

Carilah titik minimum dari fungsi >𝑓(𝑥)= 4𝑥^4 − 2𝑥^3 − 3𝑥!

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

Latihan 5 (Regresi Linear)

Lakukan pendugaan parameter regresi dengan meminimumkan jumlah kuadrat galat (residual sum of square) dari data berikut!

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

Kemudian bandingkan hasilnya dengan output dari fungsi lm!

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

Demikian, Terima Kasih

Satria June Adwendi, IPB University, sjadwendi@apps.ipb.ac.id ↩︎