Pemrograman Statistika (Visualisasi Data & Pembangkitan Bilangan Acak)

Pokok Bahasan Praktikum

Pokok bahasa yang dilatihkan dalam praktikum STA561 pertemuan 4 dan 5 adalah sebagai berikut

Praktikum 4: Visulisasi Data
- Grafik Dasar
  - Plot
  - Warna
  - Menambahkan label dan title
  - Customize Warna
  - Menambahkan amatan
  - Menambahkan garis
  - Menambahkan tanda panah
  - Menambahkan tulisan
- Grafik dengan packages ggplot2
  - Scatter plot
  - Histogram
  - Mengelola legend
  - Aesthetic Mapping
Praktikum 5: Pembangkitan Bilangan Acak
- Pembangkitan bilangan acak
- Mencari nilai peluang kumulatif peubah acak
- Mencari kuantil dari peluang peubah acak
- Mencari nilai density peubah acak
- Invers Transform Method
- Acceptance-Rejection Method
- Direct Transformation

Praktikum 4: Visualisasi Data

Selain digunakan untuk analisis statistika, R juga menyediakan berbagai fungsi untuk membuat grafik. Packages grafik standar adalahgraphics, namun terdapat beberapa package graphics lain diantaranya lattice,grid, dan ggplot2

Grafik Dasar

Fungsi umum untuk membuat grafik di R adalah plot(x,y, ...) yang selanjutkan diikuti oleh argumen-argumen.

Plot

Tipe-tipe yang biasanya digunakan untuk mem-plot:

“p” : points
“l” : lines
“b” : both
“c” : the lines part alone of “b”
“o” : overplotted
“h” : histogram like
“s” : stair steps
“S” : other steps
“n” : no plotting

Warna

Beberapa cara yang dapat digunakan untuk mengatur pewarnaan grafik, diantaranya:

Numeric 1-8 (recycle)
Character “colorname”(colors())
Fungsi rainbow(…), rgb(…)
Package “colorspace”, “colourpicker”, "RColorBrewer“

Ilustrasi 1

Plot type=“p”, menampilkan point/titik pada grafik

# plot dasar
x <- 1:40
y <- rnorm(40,5,1)
plot(x,y,type="p")

2. Plot type=“o”, menampilkan titik dan garis pada grafik

plot(x,y,type="o")

3. Plot type=“n” menampilkan kanvas kosong

plot(x,y,type="n")

Menambahkan label dan title

Secara default, pch(point character) terdapat 25 karakter

plot(x,y,type="p",xlab="Sumbu x",ylab="Sumbu y",
main="Bilangan Acak Normal",col=2,pch=16)

Customize warna

Syntax di atas menampilkan point dengan sequence warna merah hingga ungu sebanyak 40 point; dengan cex(ukuran point) sebesar 2; dengan range sumbu x 0-50 dan sumbu y 2.5-7.5

plot(x,y,type="p",xlab="Sumbu x",ylab="Sumbu y",
     main="Bilangan Acak Normal",col=rainbow(40),
     pch=16,cex=2,xlim=c(0,50),ylim=c(2.5,7.5))

Ilustrasi 2

Menambahkan Amatan

Plot di bawah ini menambahkan amatan baru pada plot sebelumnya dengan menggunakan fungsi points()

plot(x,y,type="p",xlab="Sumbu x",ylab="Sumbu y",
     main="Bilangan Acak Normal",col=rainbow(40),
     pch=16,cex=2,xlim=c(0,50),ylim=c(2.5,7.5))
x1 <- 41:50
y1 <- rnorm(10,5,1)
points(x1,y1,cex=2)

Menambahkan Garis

Penambahan garis pada plot dilakukan dengan menggunakan fungsi lines().

plot(x,y,type="p",xlab="Sumbu x",ylab="Sumbu y",
     main="Bilangan Acak Normal",col=rainbow(40),
     pch=16,cex=2,xlim=c(0,50),ylim=c(2.5,7.5))

# titik
x1 <- 41:50
y1 <- rnorm(10,5,1)
points(x1,y1,cex=2)

# garis
x2 <- rep(40.5,20)
y2 <- seq(min(c(y,y1)),max(c(y,y1)),length=20)
lines(x2,y2,col=4,lwd=2,lty=2)
abline(h=mean(y),col="red",lwd=2.5)
abline(a=2,b=1/10,col="maroon3",lwd=2,lty=2)

Menambahkan Tanda Panah

Penambahan tanda panah pada plot dilakukan dengan menggunakan fungsi arrows().

plot(x,y,type="p",xlab="Sumbu x",ylab="Sumbu y",
main="Bilangan Acak Normal",col=rainbow(40),
pch=16,cex=2,xlim=c(0,50),ylim=c(2.5,7.5))

# point
x1 <- 41:50
y1 <- rnorm(10,5,1)
points(x1,y1,cex=2)

# garis
x2 <- rep(40.5,20)
y2 <- seq(min(c(y,y1)),max(c(y,y1)),length=20)
lines(x2,y2,col=4,lwd=2,lty=2)
abline(h=mean(y),col="red",lwd=2.5)
abline(a=2,b=1/10,col="maroon3",lwd=2,lty=2)

# tanda panah
arrows(x0=30,y0=3.5,x1=40,y1=mean(y)-.1,lwd=2)

Menambahkan Tulisan

Penambahan tanda panah pada plot dilakukan dengan menggunakan fungsi text().

plot(x,y,type="p",xlab="Sumbu x",ylab="Sumbu y",
main="Bilangan Acak Normal",col=rainbow(40),
pch=16,cex=2,xlim=c(0,50),ylim=c(2.5,7.5))

# point
x1 <- 41:50
y1 <- rnorm(10,5,1)
points(x1,y1,cex=2)

# garis
x2 <- rep(40.5,20)
y2 <- seq(min(c(y,y1)),max(c(y,y1)),length=20)
lines(x2,y2,col=4,lwd=2,lty=2)
abline(h=mean(y),col="red",lwd=2.5)
abline(a=2,b=1/10,col="maroon3",lwd=2,lty=2)

# tanda panah
arrows(x0=30,y0=3.5,x1=40,y1=mean(y)-.1,lwd=2)

# tulisan
text(x=29,y=3.3,labels="Titik potong",cex=0.7)
text(x=3,y=7.3,labels="Data awal",cex=0.7)
text(x=46,y=7.3,labels="Data baru",cex=0.7)

Latihan 1

Membuat plot dari fungsi sin dengan nilai x: p < x < 2pi

plot(sin,-pi, 2*pi)

Membuat barplot dari fungsi sebaran poisson dengan lambda = 5

plot(table(rpois(100,5)),type="h",col="red",
lwd=1,main="rpois(100,lambda=5)")

Latihan 2

Menampilkan plot dengan membangkitkan bilangan acak sebaran normal sebnayak 25; dengan pch(point character) berbentuk w

a1 <- 1:25
a2 <- rnorm(25,4,2)
plot(a1,a2,pch="w",main="W")

Latihan 3

Dari plot latihan 2, ditambahkan text, label, dan modifikasi `col(warna=rainbow), cex(ukuran point=0.8)

a1 <- 1:25
a2 <- rnorm(25,4,2)
plot(a1,a2,type="n",main="W")
text(a1,a2,labels=paste("w",1:25,sep=""),col=rainbow(25),cex=0.8)

Latihan 4

Membuat program untuk menampilkan grafik dengan 100 amatanX~x^2(4) [data sebaran chi-square, derajat bebas=4].Untuk membangkitkan data yang berdistribusi chi-square dapat menggunakan fungsi rchisq()

x <- rchisq(100,df=4)
hist(x,freq=FALSE,ylim=c(0,0.2))
curve(dchisq(x,df=4),col=2,lty=2,lwd=2,add=TRUE)

Untuk menampilkan plot pada window menjadi beberapa bagian, kita dapat menggunakan fungsi 1. par(mfrow)

par(mfrow=c(2,2))

plot(1:40,y,type="p",xlab="Sumbu x",ylab="Sumbu y",
main="1. Bilangan Acak Normal",col=2,pch=16)

# plot 2
plot(sin,-pi, 2*pi, main="2. sin(x)")

plot(table(rpois(100,5)),type="h",col="red",
lwd=1,main="3. rpois(100,lambda=5)")

plot(a1,a2,type="n",main="4. W")
text(a1,a2,labels=paste("w",1:25,sep=""),
col=rainbow(25),cex=0.8)

atau 2. par(mfcol)

par(mfcol=c(2,2))
plot(1:40,y,type="p",xlab="Sumbu x",ylab="Sumbu y",
main="1. Bilangan Acak Normal",col=2,pch=16)

plot(sin,-pi, 2*pi, main="2. sin(x)")

plot(table(rpois(100,5)),type="h",col="red",lwd=1,
main="3. rpois(100,lambda=5)")

plot(a1,a2,type="n",main="4. W")
text(a1,a2,labels=paste("w",1:25,sep=""),
col=rainbow(25),cex=0.8)

Latihan 5

Membuat 4 garfik dalam 1 window dengan data 100 bilangan acak yang berdistribusi Normal(5,1)

windows()
yb <- rnorm(100,5,1)
split.screen(c(2,2))

## [1] 1 2 3 4

screen(3)

## [1] FALSE

boxplot(yb)
title("Boxplot Bilangan Acak Normal",cex.main=0.7)

screen(4)

## [1] FALSE

xb <- 1:100
plot(xb,yb,type="l",lwd=2,col="blue")
title("Line Plot Bilangan Acak Normal",cex.main=0.7)

screen(2)

## [1] FALSE

hist(yb,freq=FALSE,main=NULL,ylim=c(0,0.5))
x <- yb
curve(dnorm(x,5,1),col="red",lty=2,lwd=2,add=TRUE)
title("Histogram Bilangan Acak Normal",cex.main=0.7)

screen(1)

## [1] FALSE

plot(xb,yb,pch=16,col=rainbow(100))
title("Scatter Plot Bilangan Acak Normal",cex.main=0.7)

Grafik dengan packages ggplot2

ggplot2 merupakan implementasi dari Grammar of Graphics yang dilakukan oleh Leland Wilkinson dan merupakan buatan Hadley Wickham. ggplot2 merupakan sistem ketiga dari grafik R (sebelumnya: base dan lattice) yang tersedia di CRAN. ggplot2 memiliki aragam kelebihan untuk visualisasi data yang elegan dan kompleks.

Ilustrasi

Menenggunakan dataset Cars93. Berikut ini sekilas info mengenai dataset tersebut.

library(MASS)
library(ggplot2)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following object is masked from 'package:MASS':
## 
##     select

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

data("Cars93")
glimpse(Cars93)

## Rows: 93
## Columns: 27
## $ Manufacturer       <fct> Acura, Acura, Audi, Audi, BMW, Buick, Buick, Buick,~
## $ Model              <fct> Integra, Legend, 90, 100, 535i, Century, LeSabre, R~
## $ Type               <fct> Small, Midsize, Compact, Midsize, Midsize, Midsize,~
## $ Min.Price          <dbl> 12.9, 29.2, 25.9, 30.8, 23.7, 14.2, 19.9, 22.6, 26.~
## $ Price              <dbl> 15.9, 33.9, 29.1, 37.7, 30.0, 15.7, 20.8, 23.7, 26.~
## $ Max.Price          <dbl> 18.8, 38.7, 32.3, 44.6, 36.2, 17.3, 21.7, 24.9, 26.~
## $ MPG.city           <int> 25, 18, 20, 19, 22, 22, 19, 16, 19, 16, 16, 25, 25,~
## $ MPG.highway        <int> 31, 25, 26, 26, 30, 31, 28, 25, 27, 25, 25, 36, 34,~
## $ AirBags            <fct> None, Driver & Passenger, Driver only, Driver & Pas~
## $ DriveTrain         <fct> Front, Front, Front, Front, Rear, Front, Front, Rea~
## $ Cylinders          <fct> 4, 6, 6, 6, 4, 4, 6, 6, 6, 8, 8, 4, 4, 6, 4, 6, 6, ~
## $ EngineSize         <dbl> 1.8, 3.2, 2.8, 2.8, 3.5, 2.2, 3.8, 5.7, 3.8, 4.9, 4~
## $ Horsepower         <int> 140, 200, 172, 172, 208, 110, 170, 180, 170, 200, 2~
## $ RPM                <int> 6300, 5500, 5500, 5500, 5700, 5200, 4800, 4000, 480~
## $ Rev.per.mile       <int> 2890, 2335, 2280, 2535, 2545, 2565, 1570, 1320, 169~
## $ Man.trans.avail    <fct> Yes, Yes, Yes, Yes, Yes, No, No, No, No, No, No, Ye~
## $ Fuel.tank.capacity <dbl> 13.2, 18.0, 16.9, 21.1, 21.1, 16.4, 18.0, 23.0, 18.~
## $ Passengers         <int> 5, 5, 5, 6, 4, 6, 6, 6, 5, 6, 5, 5, 5, 4, 6, 7, 8, ~
## $ Length             <int> 177, 195, 180, 193, 186, 189, 200, 216, 198, 206, 2~
## $ Wheelbase          <int> 102, 115, 102, 106, 109, 105, 111, 116, 108, 114, 1~
## $ Width              <int> 68, 71, 67, 70, 69, 69, 74, 78, 73, 73, 74, 66, 68,~
## $ Turn.circle        <int> 37, 38, 37, 37, 39, 41, 42, 45, 41, 43, 44, 38, 39,~
## $ Rear.seat.room     <dbl> 26.5, 30.0, 28.0, 31.0, 27.0, 28.0, 30.5, 30.5, 26.~
## $ Luggage.room       <int> 11, 15, 14, 17, 13, 16, 17, 21, 14, 18, 14, 13, 14,~
## $ Weight             <int> 2705, 3560, 3375, 3405, 3640, 2880, 3470, 4105, 349~
## $ Origin             <fct> non-USA, non-USA, non-USA, non-USA, non-USA, USA, U~
## $ Make               <fct> Acura Integra, Acura Legend, Audi 90, Audi 100, BMW~

Scatter plot variabel Horsepower terhadap MPG.city

# COntoh 1
plot(MPG.city ~ Horsepower, data=Cars93)

# Contoh 2
ggplot(Cars93, aes(x=Horsepower, y=MPG.city)) +
geom_point()

Menampilkan histogram dari variabel MPG.city

# Contoh 1
hist(Cars93$MPG.city)

# Contoh 2
ggplot(Cars93, aes(x=MPG.city)) +
geom_histogram(binwidth=5)

Menambahkan legend, custom jenis plot yang berbeda berdasarkan 2 katoegori Origin

# Contoh 1
par(mar=c(4,4,.1,.1), cex.axis = 1.5, cex.lab = 1.7)

plot(MPG.city ~ Horsepower, data=Cars93,
     col=Cars93$Origin, pch=c(1,2))

legend("topright", c("USA", "non-USA"),
       title="Origin", pch=c(1,2),
       col=c("black", "red"))

# Contoh 2
ggplot(Cars93, aes(x=Horsepower, 
                   y=MPG.city,
                   color=Origin)) + geom_point(size = 3)

# Contoh 3
ggplot(Cars93, aes(x=Horsepower, y=MPG.city, color=Origin)) +
  geom_point(size = 2) +
  facet_wrap(~Cylinders)

Aesthetic Mapping

Grafik semakin interaktif dan semakin mendukung prinsip visualisasi data (ACCENT: Apprehention, Clarity, Consistency, Efficiency, Necesity, Truthfulness) dengan adanya fungsi aes() untuk melakukan penyesuaian bentuk, warna, ukuran, jenis garis, dll.

1. Aesthetic mapping terhadap objek geometri menggunakan fungsi `aes()`

ggplot(Cars93, aes(x = Horsepower, y = MPG.city)) +
  geom_point(aes(color = Cylinders, shape = Cylinders))

2. Aesthetic mapping terhadap objek geometri dengan fungsi `aes()`; dilanjutkan `geom() +`

Setiap jenis geom hanya menerima bagian aesthetic (misal untuk mengatur bentuk di aes() & tidak termasuk ketika dipetakan ke dalam geom_bar

3. Visualisasi dengan lebih dari satu `geom`

Contoh 1

g1 <- ggplot(Cars93, aes(x=Horsepower, y=MPG.city))
g2 <- g1 + geom_point(aes(color=Weight)) + geom_smooth()
g2

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

Contoh 2

g1 + geom_text(aes(label=substr(Manufacturer,1,3), size=EngineSize)) +
  geom_smooth()

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

Contoh 3 nilai di luar aes() – assigment

g1 + geom_point(color="red", size=3) + geom_smooth()

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

Contoh 4: variabel di dalam aes() – aesthetic mapping

g1 + geom_point(aes(color=Weight, shape=Origin)) +
geom_smooth()

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

4. Setiap blok ggplot2 dapat didefinisikan dengan parameter.

Default values for geom:

g <- ggplot(Cars93, aes(x = Type, y = MPG.city))
g + geom_boxplot() + geom_point()

Default values for statistics:

g <- ggplot(Cars93, aes(x = Type, y = MPG.city))
g + stat_boxplot() + geom_point()

ggplot(Cars93, aes(x = MPG.city)) + geom_histogram(binwidth=3)

5. Contoh, Flexiblity

g <- ggplot(Cars93, aes(x = MPG.city, y = Horsepower))
g + geom_point()

g + geom_line()

6. Contoh, Grouping

ggplot(Cars93, aes(MPG.city, Horsepower, shape = Cylinders)) + 
  geom_point() +
  stat_smooth(method = "lm")

## `geom_smooth()` using formula 'y ~ x'

## Warning in qt((1 - level)/2, df): NaNs produced

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf

7. Scales

Scales: Aesthetic for variabels

Scales dapat mendefinisikan beberapa indikator seperti: warna, ukuran, bentuk, jenis garis. Scales didefinisikan dan dimodifikasi dengan menggunakan fungsi scale_aesthetic_type

Contoh 1

g <- ggplot(Cars93, aes(x = Type, y = MPG.city))
g <- g + geom_point(aes(color = Weight))
g

Contoh 2

g <- g + scale_x_discrete("Cylinders",
                          labels =  c("three","four","five",
                                      "six","eight","rotary"))
g

Contoh 3

m <- max(Cars93$MPG.city)
g <- g + scale_y_continuous(breaks=c(10,20,30,m),
                            labels=c("ten (10)", "twenty (20)", "thirty (30)",
                                     paste("max (",m,")"))); g

Grafik standar di setiap grup dalam data:

Pengelompokan satu variabel : facet_wrap()
Pengelompokan dua variabel : facet_grid()

Praktikum 5: Pembakitan Bilangan Acak

Pendahuluan

Pembangkitan bilangan acak merupakan alat yang diperlukan dalam komputasi statistik dan umumnya untuk simulasi. Bilangan acak yang dibangkitkan merupakan pseudorandom (acak semu) yang memennuhi sebaran statistik tertentu. Metode bilangan acak tergantung dari pembangkitan bilangan acak uniform (seragam).

Fungsi umumnya dimulai dengan huruf r diikuti dengan nama sebaran, contoh:

Data sebaran pseudo normal: `rnorm``
Data sebaran pseudo seragam: `runif``

Fungsi-fungsi peluang suatu sebaran:

Fungsi density/mass (pdf/pmf) dimulai dengan huruf d diikuti dengan nama sebarannya, contoh: dnorm, dunif
Fungsi peluang sebaran kumulatif (cdf) dimulai huruf p diikuti dengan nama sebarannya, contoh: pnorm, punif
Fungsi quantile/invers cdf dimulai huruf q diikuti dengan nama sebarannya, contoh: qnorm, qunif

Teknik pembangkitan bilangan acak:

Invers Transform Method
Acceptance-Rejection Method
Direct Transformation

Penerapan Pembangkitan Bilangan Acak

Contoh 1: Pembangkitan bilangan acak

Pembangkitan bilangan acak normal sebanyak 10 sampel dengan mean=0 dan sd=1

x <- rnorm(10)          #x~N(0,1)
x

##  [1] -0.002043579  0.101848769 -0.290971528  0.389890530 -0.103561788
##  [6]  0.582718936  0.472764512  0.016180301 -0.885039524  0.262181940

Pembangkitan bilangan acak normal sebanyak 10 sampel dengan mean=3 dan sd=2

x1 <- rnorm(10,3,2)     #x1~N(3,sd=2)
x1

##  [1]  2.5928209  0.2954697  4.7410194  4.9746858  2.3380652  5.4524955
##  [7]  0.9093570  4.6397474  3.6966808 -1.2452946

Pembangkitan bilangan acak binomial sebanyak 10 sample

x2 <- rbinom(10,1,0.4)  #x2~bernoulli(0.4)
x2

##  [1] 0 1 0 0 0 0 1 0 0 0

Contoh 2: Mencari nilai peluang kumulatif peubah acak

Menghitung nilai peluang kumulatif dari distribusi normal P(Z<1.645)=0.95

p1 <- pnorm(1.645) 
p1

## [1] 0.9500151

Menghitung nilai peluang kumulatif dari distribusi normal P(Z<1.96)=0.975

p2 <- pnorm(1.96) 
p2

## [1] 0.9750021

p3 <- pnorm(-1.96)
p3

## [1] 0.0249979

Menghitung nilai peluang kumulatif dari distribusi f

p4 <- pf(15,df1=10,df2=15)
p4

## [1] 0.9999955

Contoh 3: Mencari nilai kuantil dari peluang peubah acak

q1 <- qnorm(0.975)
q1

## [1] 1.959964

#X~N(2,1), P(X<x)=0.95
q2 <- qnorm(0.95,2,1) 
q2

## [1] 3.644854

Contoh 4: Mencari nilai density peubah acak

# Plot density normal
a <- seq(-4,4,length=1000)
da <- dnorm(a)
plot(a,da)

Teknik Invers Transform Method

Ide dasar:
- 0 ≤ 𝐹(𝑥) ≤ 1 sehingga 𝑈=𝐹(𝑥)~𝑈𝑛𝑖𝑓𝑜𝑟𝑚(0,1)
- Jika 𝑋=𝐹^-1(𝑢) maka 𝑋~𝑓(𝑥)
Algoritma:
- Menentukan fungsi sebaran kumulatif 𝐹(𝑥) dari sebaran yang diinginkan
- Menentukan 𝐹⁻¹(𝑥)
- Membangkitkan 𝑈~𝑈𝑛𝑖𝑓𝑜𝑟𝑚(0,1)
- Transformasi 𝑋=𝐹⁻¹(𝑢)

Contoh

Membangkitkan bilangan acak Eksponensial(λ)

# Sebaran Eksponensial(lambda=3)
eks <- function(n,lambda){
   U <- runif(n)
   X <- -log(1-U)/lambda
   return(X)
}

#inverse transform method
yy1 <- eks(1000,3)   

#fungsi bawaan R
yy2 <- rexp(1000,rate=3) 

par(mfrow=c(1,2))
hist(yy1,main="Exp dari Inverse Transform")
hist(yy2,main="Exp dari fungsi rexp")

Teknik Acceptance-Rejection Method

Algoritma untuk membangkitkan 𝑋~𝑓(𝑥); 𝑥₀ < 𝑥 < 𝑥₁ :
- Membangkitkan 𝑥~𝑈(𝑥₀, 𝑥₁)
- Menentukan 𝑐 sehingga 𝑓(𝑥) ≤ 𝑐 untuk 𝑥₀ < 𝑥 < 𝑥₁
- Membangkitkan 𝑦₁~𝑈(0, 𝑐)
- Tentukan 𝑦₂ = 𝑓(𝑥)
- Jika 𝑦₁ ≤ 𝑦₂, terima 𝑥

Contoh

Bangkitkan bilangan acak yang memiliki fungsi kumulatif probability 𝑓_Y(𝑦)=3𝑦²; 0 < 𝑦 < 1 menggunakan acceptance-rejection method!

ar <- function(n,x0,x1,f) {
   xx <- seq(x0,x1,length=10000)
   c <- max(f(xx))
   terima <- 0; iterasi <- 0
   hasil <- numeric(n)
   while(terima<n) {
      x <- runif(1,x0,x1)
      y1<- runif(1,0,c)
      y2<- f(x)
      if(y1<=y2) {
         terima <- terima+1
         hasil[terima]<-x }
      iterasi <- iterasi+1 }
   list(hasil=hasil,iterasi=iterasi)
}

set.seed(10)
f <- function(x) {3*x^2}
yyy <- ar(100,0,1,f)
yyy

## $hasil
##   [1] 0.8647212 0.7751099 0.8382877 0.7707715 0.5355970 0.8613824 0.2036477
##   [8] 0.7979930 0.7438394 0.3443435 0.9837322 0.6935082 0.6331153 0.8880315
##  [15] 0.7690405 0.6483695 0.8795432 0.9360689 0.7233519 0.7620444 0.9868082
##  [22] 0.8760261 0.7240640 0.8140516 0.5588949 0.8900940 0.7456896 0.8480646
##  [29] 0.8703302 0.8223331 0.8508123 0.7709219 0.8953595 0.5803863 0.5982260
##  [36] 0.9235285 0.7367755 0.6898170 0.8301572 0.9293209 0.9095163 0.5347576
##  [43] 0.3478601 0.8759762 0.7286815 0.8749293 0.6988356 0.8312562 0.5572238
##  [50] 0.6647687 0.7400502 0.9806898 0.3800746 0.7553169 0.5184889 0.8879149
##  [57] 0.9177773 0.8084086 0.8537441 0.4232184 0.7604306 0.3405763 0.3886568
##  [64] 0.4774175 0.5387605 0.9485434 0.7124685 0.9081691 0.9457656 0.7716899
##  [71] 0.6946655 0.5368832 0.8481593 0.8242752 0.5123742 0.3152032 0.9924487
##  [78] 0.9327120 0.9892809 0.6283590 0.5254605 0.8810815 0.5291748 0.5765517
##  [85] 0.7231807 0.8761180 0.3995670 0.8986123 0.9335217 0.7859216 0.7784128
##  [92] 0.6955333 0.9060413 0.9916424 0.4729846 0.9770567 0.9386110 0.9959093
##  [99] 0.8543663 0.8309547
## 
## $iterasi
## [1] 322

Grafik

par(mfrow=c(1,1))
hist(yyy$hasil, 
main="f(x)=3(x^2)", prob=T)
x <- seq(0, 1, 0.01)
lines(x, f(x), lwd=2, col=4)

Teknik Direct Transformation

Ide: memanfaatkan beberapa fungsi transformasi dari berbagai sebaran yang ada

Contoh

Jika X~𝑈(0,1) menggunakan 𝑌 = (𝑏−𝑎)𝑋+𝑎maka 𝑌~𝑈(𝑎,𝑏)

Membangkitkan bilangan acak 𝑌~𝑈(3,5) sehingga 𝑌= (5-3)𝑋+3 -> 𝑌= 2𝑋+3

x <- runif(1000)
y <- 2*x+3
y

##    [1] 4.595651 4.168212 3.865580 3.368448 3.685894 3.168929 4.232466 3.471525
##    [9] 3.816635 4.432023 3.654097 3.025882 4.017800 3.645034 4.035468 4.220792
##   [17] 4.369285 4.777494 4.474976 4.379829 3.045420 3.458973 3.489782 4.082726
##   [25] 3.896056 4.563540 3.730657 3.726081 3.958200 3.236057 4.076578 3.582958
##   [33] 4.210074 3.726026 4.420964 3.973593 4.783335 3.698059 4.224458 3.630964
##   [41] 3.646383 3.750261 4.467030 3.371470 3.469069 4.693028 4.894643 4.280374
##   [49] 3.537450 4.790113 4.299423 4.519895 4.965529 3.051509 3.303683 3.906797
##   [57] 3.375193 4.141702 4.796323 3.377092 3.108475 3.419335 4.738254 4.080452
##   [65] 4.969106 4.856840 4.331528 3.915997 4.770318 3.073207 4.697783 4.413223
##   [73] 4.486807 4.246549 4.955789 4.387593 4.726321 3.386577 4.532412 4.993235
##   [81] 3.161616 3.375557 3.197817 4.397423 3.332837 3.218048 3.082443 4.578633
##   [89] 4.594420 4.758882 3.447881 4.580076 4.968844 3.084786 3.672106 3.884082
##   [97] 3.995492 3.071297 4.311110 3.813428 3.714606 3.236832 3.667819 3.942326
##  [105] 4.193309 4.102471 4.973566 4.951496 4.246684 3.942961 3.572900 4.492391
##  [113] 3.137360 3.982451 4.322077 3.751640 3.597455 4.001769 3.725986 4.395659
##  [121] 4.534990 4.707527 3.577784 4.937220 3.147125 3.566941 3.292801 3.264311
##  [129] 3.494709 3.345439 3.539915 3.026146 4.941346 4.853478 3.884387 3.050453
##  [137] 3.522547 3.310321 3.680943 4.031377 4.268245 4.580220 3.257111 3.599641
##  [145] 4.052678 4.027126 4.030312 4.403521 4.642983 4.886053 4.967236 3.922478
##  [153] 3.548622 3.377231 4.786665 4.642731 3.936669 4.018760 4.762696 3.363010
##  [161] 4.971195 4.327953 4.315386 3.614140 3.460296 3.062097 3.050408 3.744086
##  [169] 3.051171 4.354773 3.346706 3.725468 4.769000 4.526478 3.532470 3.834602
##  [177] 3.894296 4.798221 3.986260 3.234679 3.645782 3.086577 4.859232 3.050443
##  [185] 4.969942 3.187620 3.002687 4.369433 3.076397 4.112787 3.715231 3.806230
##  [193] 4.291993 3.315377 3.217174 4.120446 4.365137 3.564720 4.782446 3.341252
##  [201] 3.553095 4.455447 4.548107 3.554820 4.459492 4.290764 3.819135 3.092542
##  [209] 3.480910 3.737291 4.781122 3.456868 4.871115 4.384660 4.821405 3.713408
##  [217] 3.901560 3.308500 4.649646 4.521387 4.015555 4.857621 4.593943 4.782217
##  [225] 3.760751 3.087453 3.952728 4.439325 4.686256 4.284194 3.168267 4.869526
##  [233] 3.315281 4.150718 3.981427 4.437502 4.878076 4.042365 3.447635 4.644413
##  [241] 4.685966 4.115470 4.216081 3.375351 3.935933 4.295949 3.366233 3.346195
##  [249] 4.363890 3.628338 4.835012 3.815494 4.027718 4.390839 3.196155 4.508844
##  [257] 4.755563 3.952660 3.127663 3.486304 3.011791 4.950004 3.479238 4.126829
##  [265] 3.773781 3.515035 3.663339 3.868293 3.731539 4.564783 3.968643 4.505326
##  [273] 4.624903 3.351114 3.718886 3.149059 3.466030 3.876820 4.920182 3.604485
##  [281] 3.931436 3.851931 4.731431 4.749438 3.984699 3.443737 3.135317 4.509235
##  [289] 4.993079 3.632973 3.457061 4.919519 4.713212 4.893936 3.974407 3.498106
##  [297] 3.589225 3.288184 3.054149 4.625688 4.686465 3.231876 3.780065 4.845464
##  [305] 3.211458 3.840869 4.788149 4.815876 3.847575 4.974608 4.638219 4.271901
##  [313] 3.030907 4.954875 4.802814 3.227271 3.668571 3.904605 3.798721 4.781529
##  [321] 4.815544 3.656805 4.598357 4.052858 3.142921 3.067473 3.537728 4.267681
##  [329] 4.084502 4.700150 3.065181 4.197046 3.106003 4.154466 4.180102 4.051311
##  [337] 4.643063 3.554346 3.740491 4.782274 4.801987 4.384187 4.045663 3.493019
##  [345] 3.027529 3.214891 3.601037 3.660628 4.720207 4.644854 3.807896 4.375866
##  [353] 3.363114 3.900362 3.385230 4.943052 4.615413 4.420358 3.496498 4.019566
##  [361] 4.137517 3.663136 3.873324 3.399394 4.572513 4.569445 4.909801 4.239187
##  [369] 4.785013 3.465037 4.023484 4.060217 4.902524 4.148773 4.742039 4.878590
##  [377] 3.160793 4.536739 4.124574 3.323276 4.614337 3.178618 3.973487 4.737529
##  [385] 4.659949 4.963217 4.741565 3.446988 4.306195 3.833228 4.293684 3.338345
##  [393] 3.801240 3.163310 4.385805 4.209765 3.183583 4.442544 3.905830 4.058087
##  [401] 3.526571 3.962961 4.502033 3.884139 3.765793 4.085434 4.469945 3.784377
##  [409] 4.906737 3.246651 3.016395 3.021085 3.177066 4.960241 4.770640 3.963674
##  [417] 4.804280 4.375579 3.952479 4.962703 3.493305 4.968655 3.146102 3.376815
##  [425] 4.804378 3.767460 4.760473 3.138694 4.612437 4.526004 3.021634 3.744841
##  [433] 3.163980 3.631048 3.165785 3.683686 3.283250 3.737113 3.346487 4.695052
##  [441] 3.557584 3.530328 4.699201 4.113038 4.353739 4.803387 3.232695 4.243600
##  [449] 4.571575 4.263495 3.076344 3.388704 3.556719 3.713973 3.199713 3.681986
##  [457] 3.754938 3.085458 3.484692 4.643762 3.008197 4.676406 4.590641 3.282191
##  [465] 4.241101 3.690541 3.953717 4.524971 4.645850 3.423704 3.826890 4.164745
##  [473] 4.679355 4.809141 3.054438 3.643403 3.292472 3.574535 4.626748 3.457478
##  [481] 4.999602 4.428614 4.541616 4.046089 3.835642 4.207878 4.142799 3.668253
##  [489] 4.846237 3.173952 3.657645 3.948126 4.286523 3.747317 4.513459 4.484815
##  [497] 4.649410 4.672014 4.864272 4.989662 4.753757 3.485692 4.179418 4.599828
##  [505] 3.195160 4.490126 3.949943 4.482142 4.679428 3.352180 3.314449 3.395827
##  [513] 3.452662 3.902255 4.918844 3.761364 3.234991 3.361696 3.676700 3.658317
##  [521] 3.116457 3.073477 3.779114 4.891029 3.558328 3.567479 4.490951 3.030119
##  [529] 3.523791 4.140343 3.589872 4.926635 4.658660 4.447436 3.915916 3.447092
##  [537] 3.708353 3.577917 3.860111 3.520046 4.217823 3.469950 4.186825 4.009659
##  [545] 4.960761 4.663596 3.197030 3.339365 4.432129 4.634356 4.925802 4.272933
##  [553] 4.168626 4.954075 4.623677 4.177371 3.367073 3.954232 3.956338 4.033501
##  [561] 3.174859 4.450297 4.248931 3.222602 4.837201 4.414659 3.225810 3.461793
##  [569] 3.838243 3.014711 3.120554 3.293754 4.395576 4.572241 3.675333 3.467193
##  [577] 3.087367 3.367471 4.986978 4.200105 4.381540 3.404015 3.500395 3.453344
##  [585] 3.375089 3.230314 3.083384 4.929557 4.454521 3.762533 3.383604 4.935798
##  [593] 4.531758 3.546295 3.975483 3.449228 3.906684 3.617831 4.078245 3.133104
##  [601] 3.498497 3.095699 3.960723 3.293161 3.574989 3.776852 4.983567 4.034965
##  [609] 3.928389 3.297234 3.411277 4.795526 3.774142 3.597192 4.104445 3.426783
##  [617] 3.288198 4.798747 4.536796 4.074385 4.253072 3.400720 3.727104 4.490435
##  [625] 4.800692 3.538431 3.310350 3.387320 3.610006 4.521122 3.867236 4.123544
##  [633] 3.762143 3.619385 4.957283 4.589340 4.632362 3.734138 4.107727 3.190931
##  [641] 3.840473 3.063554 3.419761 4.101204 4.480505 3.335551 3.736374 4.556347
##  [649] 3.083326 4.414239 4.206617 3.770268 3.359301 3.676023 3.569652 3.155290
##  [657] 3.386589 3.570812 4.590909 4.524650 3.423634 3.188275 4.670851 4.008845
##  [665] 4.993068 3.570497 3.892304 4.282780 3.915558 3.453340 3.178928 3.771504
##  [673] 4.342322 3.476966 4.663507 3.074694 3.809217 3.975703 4.066113 4.942717
##  [681] 3.717558 4.438593 4.261515 4.868353 4.209460 4.378453 4.997412 4.186681
##  [689] 3.730297 3.844101 4.712017 4.604720 3.465583 4.094532 4.045095 4.956168
##  [697] 4.360094 3.730014 3.709412 3.571868 4.166459 3.285475 4.802813 3.410562
##  [705] 3.103771 4.930671 3.433312 3.629858 4.390588 3.646413 4.776634 3.443542
##  [713] 4.239717 3.058778 4.166561 4.556808 4.360835 3.151751 3.056866 4.371048
##  [721] 4.784433 4.894728 4.304202 4.876747 3.975111 3.056864 3.105490 4.607493
##  [729] 4.348202 3.825592 3.489664 4.047145 4.840111 4.381992 3.088879 4.344586
##  [737] 4.152450 4.382661 3.482083 3.152124 3.418771 3.672935 3.523344 4.098084
##  [745] 4.917674 3.125032 3.023792 3.352062 3.106859 4.448480 4.486721 4.953225
##  [753] 4.633411 3.515097 4.549516 3.502145 3.820571 3.471089 3.486448 4.379599
##  [761] 4.387598 3.217931 4.462260 3.665351 3.095538 3.702529 3.736679 4.759911
##  [769] 4.387971 3.180551 3.089623 4.191134 3.836067 3.263177 3.690190 3.049972
##  [777] 3.234589 4.124177 4.349138 4.343698 4.164451 4.373908 3.226232 3.741214
##  [785] 4.557165 4.857548 3.549154 3.417801 3.118530 3.808807 3.319413 4.696115
##  [793] 4.725716 4.473617 4.640932 4.620981 3.439493 4.008431 3.634191 3.146057
##  [801] 3.258180 3.486712 3.668043 3.569197 3.401142 3.357883 3.547849 4.657281
##  [809] 3.438995 3.970205 3.372129 4.077384 3.971496 4.567775 3.760992 4.005694
##  [817] 3.002594 3.857844 3.044588 3.937043 3.963014 4.170063 3.719193 3.152707
##  [825] 3.845376 3.627579 3.129907 3.438821 3.372375 4.916585 3.388299 4.098514
##  [833] 4.419020 4.473513 3.850693 3.426815 3.745070 3.338124 3.591310 3.336311
##  [841] 3.710282 4.104672 4.530202 4.209216 4.772064 4.567791 4.763648 3.686085
##  [849] 4.497901 4.330342 4.760206 3.460885 4.149536 3.826504 4.382265 4.910413
##  [857] 3.051461 4.665443 4.326053 3.222537 4.759395 3.248831 3.148112 4.744803
##  [865] 3.713805 3.605298 3.543323 3.536265 4.250934 4.321638 4.728602 3.304895
##  [873] 3.684716 4.978813 4.984080 4.507777 4.753503 4.254134 4.264930 3.120719
##  [881] 4.351738 3.016930 3.753487 3.608348 4.839784 3.526077 3.681608 3.255639
##  [889] 3.929795 3.787575 4.254725 3.102797 4.764022 3.313598 3.798330 3.102817
##  [897] 4.623546 4.145303 3.625152 4.554414 3.466996 4.238963 3.234733 3.263493
##  [905] 3.590617 4.564968 4.095003 3.331838 4.812556 4.309262 4.851002 3.554039
##  [913] 4.196659 3.473372 4.585330 4.985140 3.462578 4.505041 4.203711 4.297979
##  [921] 3.393737 4.094990 4.170738 4.369566 3.478194 4.223525 4.522832 3.555257
##  [929] 3.654607 4.375362 3.735138 3.017994 4.395007 3.065458 4.078626 3.762759
##  [937] 4.542933 3.843434 3.780590 4.038624 4.410163 3.913584 3.929607 4.931056
##  [945] 3.455814 4.747601 3.118741 4.301843 4.721172 3.875542 3.019229 3.604486
##  [953] 4.163448 3.752417 4.555575 3.242685 3.033581 3.864228 4.397758 3.879327
##  [961] 3.962973 3.970135 3.250362 3.526578 4.470872 3.082856 3.467550 3.220450
##  [969] 3.699188 3.228758 3.501766 3.592546 4.632085 4.836780 4.624836 4.725689
##  [977] 4.303669 4.760358 3.353661 4.767130 4.544235 4.009297 4.379276 4.318141
##  [985] 3.982503 3.481575 4.443483 4.440020 4.550645 3.715663 3.621288 4.462299
##  [993] 3.463366 3.380312 3.999186 4.768694 4.579657 4.360228 4.732703 3.965204

Pembangkitan Bilangan Acak untuk Model

Analisis Regresi

Contoh 1 1. Membangkitkan data untuk variabel penjelas X1 dan X2 serta variabel respon Y

set.seed(123)

X1 <- runif(25,10,25)
X2 <- runif(25,90,200)
Y <- 10 + 2.3*X1 + 0.7*X2 + rnorm(25,0,9)

latihan <- data.frame(X1, X2, Y)
latihan

##          X1        X2        Y
## 1  14.31366 167.93835 145.2980
## 2  21.82458 149.84726 172.6297
## 3  16.13465 155.35562 157.2390
## 4  23.24526 121.80757 138.4862
## 5  24.10701 106.18250 151.0582
## 6  10.68335 195.93267 175.5627
## 7  17.92158 189.25289 181.0410
## 8  23.38629 165.97758 188.0289
## 9  18.27153 177.50142 184.1787
## 10 16.84922  92.70751 121.0427
## 11 24.35250 142.55756 171.9988
## 12 16.80001 173.43055 175.0267
## 13 20.16356 113.80487 135.4824
## 14 18.58950 124.99991 137.5021
## 15 11.54387 115.47884 113.9618
## 16 23.49737 105.70800 131.7872
## 17 13.69132 135.60010 134.5388
## 18 10.63089 135.50968 117.9193
## 19 14.91881 130.57300 155.2350
## 20 24.31755 106.76892 151.5403
## 21 23.34309 105.26867 127.2692
## 22 20.39205 115.63375 134.2194
## 23 19.60760 141.25587 149.7767
## 24 24.91405 119.25699 157.8019
## 25 19.83559 184.36105 183.9243

2. Estimasi Model regresi Linier

Model 1

model1 <- lm(Y ~ X1, data=latihan)
summary(model1)

## 
## Call:
## lm(formula = Y ~ X1, data = latihan)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -29.991 -17.738  -2.632  17.896  33.171 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  131.846     19.741   6.679 8.18e-07 ***
## X1             1.049      1.015   1.033    0.312    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 22.42 on 23 degrees of freedom
## Multiple R-squared:  0.04433,    Adjusted R-squared:  0.002775 
## F-statistic: 1.067 on 1 and 23 DF,  p-value: 0.3124

Model 2

model2 <- lm(Y ~ X2, data=latihan)
summary(model2)

## 
## Call:
## lm(formula = Y ~ X2, data = latihan)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -31.804  -4.850  -1.606  10.011  20.569 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 70.83913   13.86879   5.108 3.57e-05 ***
## X2           0.58213    0.09767   5.960 4.46e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.38 on 23 degrees of freedom
## Multiple R-squared:  0.607,  Adjusted R-squared:  0.5899 
## F-statistic: 35.52 on 1 and 23 DF,  p-value: 4.464e-06

Model 3

model3 <- lm(Y ~ X1 + X2, data=latihan)
summary(model3)

## 
## Call:
## lm(formula = Y ~ X1 + X2, data = latihan)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.8291  -4.5994  -0.4576   5.0602  20.5307 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.62881   13.40229  -0.047    0.963    
## X1           2.72951    0.40701   6.706 9.67e-07 ***
## X2           0.72459    0.06105  11.869 4.91e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.427 on 22 degrees of freedom
## Multiple R-squared:  0.8709, Adjusted R-squared:  0.8592 
## F-statistic: 74.21 on 2 and 22 DF,  p-value: 1.66e-10

Model 4

model4 <- lm(Y ~ X1 * X2, data=latihan)
summary(model4)

## 
## Call:
## lm(formula = Y ~ X1 * X2, data = latihan)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.4425  -4.5808  -0.6702   4.6431  20.2409 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 11.625896  37.941991   0.306   0.7623  
## X1           2.063573   1.967526   1.049   0.3062  
## X2           0.635870   0.263687   2.411   0.0251 *
## X1:X2        0.004905   0.014165   0.346   0.7326  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.6 on 21 degrees of freedom
## Multiple R-squared:  0.8716, Adjusted R-squared:  0.8533 
## F-statistic: 47.53 on 3 and 21 DF,  p-value: 1.548e-09

Model 5

model5 <- lm(Y ~ X1 : X2, data=latihan)
summary(model5)

## 
## Call:
## lm(formula = Y ~ X1:X2, data = latihan)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -20.929  -7.009  -1.430   2.856  37.374 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 80.645221  10.481428   7.694 8.32e-08 ***
## X1:X2        0.027491   0.003929   6.997 3.94e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 12.97 on 23 degrees of freedom
## Multiple R-squared:  0.6804, Adjusted R-squared:  0.6665 
## F-statistic: 48.96 on 1 and 23 DF,  p-value: 3.942e-07

Ringkasan Model

# Nilai Rquared
Rsquared <- matrix(c(summary(model1)[[8]], summary(model1)[[9]],
                     summary(model2)[[8]], summary(model2)[[9]],
                     summary(model3)[[8]], summary(model3)[[9]],
                     summary(model4)[[8]], summary(model4)[[9]],
                     summary(model5)[[8]], summary(model5)[[9]]),5, byrow=TRUE)
colnames(Rsquared) <-c("Rsquared","Adjusted Rquared"); Rsquared*100

##       Rsquared Adjusted Rquared
## [1,]  4.432592        0.2774876
## [2,] 60.699195       58.9904640
## [3,] 87.090357       85.9167529
## [4,] 87.163648       85.3298839
## [5,] 68.036265       66.6465370

4. Model terbaik

Terpilih model terbaik adalah model3

# Koefisien dari estimasi model 3
koef3 <- coef(model3)

# Selang kepercayaan
sk3 <- confint(model3)

cbind(koef3, sk3)

##                  koef3       2.5 %     97.5 %
## (Intercept) -0.6288095 -28.4234482 27.1658292
## X1           2.7295126   1.8854322  3.5735929
## X2           0.7245911   0.5979779  0.8512044

5. Uji Anova

anova(model3)

## Analysis of Variance Table
## 
## Response: Y
##           Df  Sum Sq Mean Sq  F value    Pr(>F)    
## X1         1   536.4   536.4   7.5538   0.01173 *  
## X2         1 10002.4 10002.4 140.8614 4.911e-11 ***
## Residuals 22  1562.2    71.0                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

6. Diagnostik Model

par(mfrow=c(2,2)); plot(model3)

Latihan

Contoh 1 Bangkitkan peubah Y, X1, X2, dan X3 berdasarkan model regresi liniear berganda berikut ini:

Y = 10 + 3X₁ + 5X₂ + 7X₃ + ϵ dengan mengasumsikan bahwa ϵ ~ N(0,1). Banyaknya amatan yang dibangkitkan adalah 1000 amatan.

Penyelesaian

Memnagkitkan residual berdasarkan distribusi normal ϵ ~ N(0,1)

n <- 1000
set.seed(12)
epsilon <- rnorm(n, mean=0, sd=1)
head(epsilon)

## [1] -1.4805676  1.5771695 -0.9567445 -0.9200052 -1.9976421 -0.2722960

Membangkitkan peubah penjelas X_i yang saling bebas

set.seed(10)
X <- replicate(3, rnorm(n, mean=0, sd=1.5))
head(X)

##             [,1]       [,2]        [,3]
## [1,]  0.02811926  1.5750205 -0.46179747
## [2,] -0.27637881  0.4291389  1.13712837
## [3,] -2.05699582  0.3608472 -0.86079512
## [4,] -0.89875157  1.2490578 -1.40811669
## [5,]  0.44181769 -0.3344748 -0.04154899
## [6,]  0.58469145  0.4325163 -1.59937298

Membangkitkan peubah respon Y dengan menggunakan model regresi

betas <- matrix(c(10,3,5,7), nrow=4, ncol=1)

# nilai peubah penjelas X_i
Xgab <- cbind(1,X)   
head(Xgab)

##      [,1]        [,2]       [,3]        [,4]
## [1,]    1  0.02811926  1.5750205 -0.46179747
## [2,]    1 -0.27637881  0.4291389  1.13712837
## [3,]    1 -2.05699582  0.3608472 -0.86079512
## [4,]    1 -0.89875157  1.2490578 -1.40811669
## [5,]    1  0.44181769 -0.3344748 -0.04154899
## [6,]    1  0.58469145  0.4325163 -1.59937298

Ket: Xgab <- cbind(1,X digunakan untuk mengabungkan peubah penjelas yang terdapat pada matrix X dengan vektor yang isinya 1. Vektor tersebut digunakan sebagai pengali konstanta model (dalam hal ini bernilai 10).

# nilai peubah respon
Y <- Xgab %*% betas + epsilon
head(Y)

##           [,1]
## [1,] 13.246310
## [2,] 20.853626
## [3,] -1.349062
## [4,]  2.772212
## [5,]  7.364594
## [6,]  2.448749

# menggabungkan X_i dan Y menjadi sbeuah dataframe (dataregresi)
dataregresi <- data.frame(Y,X)
colnames(dataregresi) <- c("Y", "X1", "X2", "X3")
head(dataregresi)

##           Y          X1         X2          X3
## 1 13.246310  0.02811926  1.5750205 -0.46179747
## 2 20.853626 -0.27637881  0.4291389  1.13712837
## 3 -1.349062 -2.05699582  0.3608472 -0.86079512
## 4  2.772212 -0.89875157  1.2490578 -1.40811669
## 5  7.364594  0.44181769 -0.3344748 -0.04154899
## 6  2.448749  0.58469145  0.4325163 -1.59937298

Membuktikan bahwa pembangkitan yang dilakukan benar, maka kita akan memodelkan data di atas dengan fungsi lm (linier modelling)

modelregresi <- lm(Y ~ X1 + X2 + X3, data=dataregresi)
summary(modelregresi)

## 
## Call:
## lm(formula = Y ~ X1 + X2 + X3, data = dataregresi)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.0364 -0.6052 -0.0207  0.6087  3.1774 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  9.97484    0.03033   328.8   <2e-16 ***
## X1           2.97225    0.02042   145.6   <2e-16 ***
## X2           4.97282    0.01939   256.5   <2e-16 ***
## X3           7.00493    0.01989   352.2   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9589 on 996 degrees of freedom
## Multiple R-squared:  0.9955, Adjusted R-squared:  0.9954 
## F-statistic: 7.275e+04 on 3 and 996 DF,  p-value: < 2.2e-16

Berdasarkan output regresi diatas, nilai koefisien regresi yang dihasilkan mendekati nilai koefisien regresi yang kita bangkitkan. Hal ini membuktikan bahwa data sudah sesuai dengan model regresi yang kita inginkan.

Nilai koefisien regresi yang dibangkitkan Y = 10 + 3X₁ + 5X₂ + 7X₃ + ϵ

Nilai koefisien regresi model Y = 9.97 + 2.98X1 + 4.98X2 + 7.01X3

Referensi

Raharjo, Mulianto. (16 September 2021). STA 561 Pemrograman Statistika: Pembangkitan Bilangan Acak. Diakses dari https://newlms.ipb.ac.id/

Soleh. A.M. (2021). STA 561 Pemrograman Statistika. Pembangkitan Bilangan Acak. Diakses dari https://newlms.ipb.ac.id/

Dito, G.A. Praktikum Komputasi Statistika. Diakses dari https://bookdown.org/gerryalfadito/praktikum-komputasi-statistika/

Pemrograman Statistika (Visualisasi Data & Pembangkitan Bilangan Acak)

Pokok Bahasan Praktikum

Praktikum 4: Visualisasi Data

Grafik Dasar

Plot

Warna

Ilustrasi 1

Ilustrasi 2

Latihan 1

Latihan 2

Latihan 3

Latihan 4

Latihan 5

Grafik dengan packages ggplot2

Ilustrasi

Aesthetic Mapping

1. Aesthetic mapping terhadap objek geometri menggunakan fungsi aes()

2. Aesthetic mapping terhadap objek geometri dengan fungsi aes(); dilanjutkan geom() +

3. Visualisasi dengan lebih dari satu geom

4. Setiap blok ggplot2 dapat didefinisikan dengan parameter.

5. Contoh, Flexiblity

6. Contoh, Grouping

7. Scales

Praktikum 5: Pembakitan Bilangan Acak

Pendahuluan

Penerapan Pembangkitan Bilangan Acak

Teknik Invers Transform Method

Contoh

Teknik Acceptance-Rejection Method

Contoh

Teknik Direct Transformation

Contoh

Pembangkitan Bilangan Acak untuk Model

Analisis Regresi

Latihan

Referensi

1. Aesthetic mapping terhadap objek geometri menggunakan fungsi `aes()`

2. Aesthetic mapping terhadap objek geometri dengan fungsi `aes()`; dilanjutkan `geom() +`

3. Visualisasi dengan lebih dari satu `geom`