Tipe-tipe plot

Akan dibuat plot antara data x bilangan berurutan 1 sampai 40 dan y yaitu angka acak dengan rata-rata 5 dan standar deviasi 1

x <- 1:40
y <- rnorm(40,5,1)

Tipe Titik

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

Tipe Garis

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

Tipe both (Titik dan Garis Terpisah)

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

Tipe c (Titik kosong dihubungkan dengan garis)

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

Tipe Overplot (Titik dan Garis)

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

Tipe Tangga : Datar terlebih dahulu

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

Tipe S (S huruf kapital) : Tipe tangga lainnya

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

Tipe Histogram

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

Tipe Noplotting : Membuat kanvas kosong

Skala kanvas mengikuti codingan titik

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

Pilihan Garis pada Grafik

Bentuk garis ini dapat diatur dengan menggunakan argumen lty.

Berikut adalah daftar dari isi argumen tersebut beserta keterangannya:

Pilihan Warna Grafik

Menggunakan angka integer 1-8

Warna ini akan berulang jika lebih dari 8, misal angka 9 akan kembali ke warna hitam, 10 akan berwarna merah, dan seterusnya.

Cara akses dengan dengan syntax col=angka, misal : col=1

Menggunakan Fungsi colors()

Ada 657 warna. Dapat diakses menggunakan “colorname” sesuai dengan nama warna yang tersedia pada fungsi colors()

Menggunakan fungsi rainbow()

Memberikan gradasi warna dari merah sampai ungu tergantung pada banyaknya nilai yang kita masukkan.

rainbow(10)

##  [1] "#FF0000" "#FF9900" "#CCFF00" "#33FF00" "#00FF66" "#00FFFF" "#0066FF"
##  [8] "#3300FF" "#CC00FF" "#FF0099"

Akan menghasilkan 10 gradasi warna dari merah sampai ungu.

Menggunakan fungsi rgb()

Ada pilihan 0-255 warna. Caranya dengan memasukkan angka untuk masing-masing warna red, green , dan blue yang akan menghasilkan code hexacode untuk mewarnai grafik.

rgb(red=0.30, green=0.50, blue=0.250)

## [1] "#4D8040"

Contoh Penggunaan

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

cex digunakan untuk mengatur ukuran point. defaultnya adalah 1. semakin besar angka akan semakin besar titiknya xlim dan ylim untuk mengatur skala yang ada di grafik. nilai dari x adalah 0-100, misalnya mau dibuat x nya sampai 150

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

Menambahkan Objek pada Grafik

Menambah Amatan

Akan ditambahkan amatan dalam plot yang sudah ada berupa titik (dengan fungsi points() )sebanyak 10 titik antara x1 yaitu angka berurutan 41 sampai 50 dan angka acak random normal dengan mean=5 dan sd=1

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

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

Data baru adalah data titik dengan warna putih yang ada di sebelah kanan.

Menambahkan garis

digunakan fungsi lines. Penambahan garis berguna untuk membedakan atau memberi batasan pada data, misalnya: data training dan data testing

lwd : line width lty : integer

Akan ditambahkan garis x2 yaitu angka x 40.5 sebanyak 20 titik. dan y2 adalah angka berurut dari y1 namun hanya diambil 20 titik pertama.

x2 <- rep(40.5,20)
y2 <- seq(min(c(y,y1)), max(c(y,y1)), length=20)
x2

##  [1] 40.5 40.5 40.5 40.5 40.5 40.5 40.5 40.5 40.5 40.5 40.5 40.5 40.5 40.5 40.5
## [16] 40.5 40.5 40.5 40.5 40.5

y2

##  [1] 2.256281 2.545201 2.834121 3.123041 3.411962 3.700882 3.989802 4.278722
##  [9] 4.567642 4.856562 5.145482 5.434402 5.723322 6.012242 6.301162 6.590082
## [17] 6.879002 7.167922 7.456842 7.745762

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

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

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

menambahkan garis dengan fungsi abline. Akan ditambahkan garis kedua yaitu garis horisontal rata-rata dari y dengan warna merah. Opsi abline untu membuat garis horisontal adalah dengan fungsi h (garis horisontal x=sekian)

y

##  [1] 6.336363 6.304346 5.803345 6.843525 6.406804 7.001023 5.753731 5.568186
##  [9] 4.954513 3.909676 5.007701 5.656115 5.821170 7.745762 3.989398 4.747912
## [17] 4.317158 5.076852 5.510257 5.337317 4.676833 5.571877 6.222197 4.689398
## [25] 4.774655 6.443215 5.642894 5.541299 3.223438 3.730646 5.798377 4.438596
## [33] 6.021566 4.320261 5.579155 3.981210 4.482418 4.589406 4.699340 2.256281

h=mean(y)
h

## [1] 5.219355

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

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

#Menambah garis 1
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)

#Menambah garis 2
abline(h=mean(y), col="red", lwd=2.5)

Akan ditambahkan garis ke-3 yaitu suatu garis dengan persamaan a+bx dimana a=2 dan b=1/10 dengan nilai x yang berada pada rentang kanvas dan diberi warna kuning

#PLOT AWAL
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))

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

#Menambah garis 1
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)

#Menambah garis 2
abline(h=mean(y), col="red", lwd=2.5)

#Menambah garis 3
abline(a=2,b=1/10, col="yellow", lwd=2, lty=2)

Menambahkan Tanda Panah

arrows perlu ditentukan koordinat awal dan akhir

#PLOT AWAL
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))

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

#Menambah garis 1
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)

#Menambah garis 2
abline(h=mean(y), col="red", lwd=2.5)

#Menambah garis 3
abline(a=2,b=1/10, col="yellow", lwd=2, lty=2)

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

Menambahkan tulisan

Akan ditambahkan tulisan “Titik potong”

X dan y mengatur koordinat text

labels mengatur isi text

cex untuk mengatur besar kecilnya tulisan

#PLOT AWAL
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))

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

#Menambah garis 1
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)

#Menambah garis 2
abline(h=mean(y), col="red", lwd=2.5)

#Menambah garis 3
abline(a=2,b=1/10, col="yellow", lwd=2, lty=2)

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

#Menambah tulisan

text(x=29,y=3.3, labels="Titik Potong",cex=0.7)

#PLOT AWAL
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))

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

#Menambah garis 1
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)

#Menambah garis 2
abline(h=mean(y), col="red", lwd=2.5)

#Menambah garis 3
abline(a=2,b=1/10, col="yellow", lwd=2, lty=2)

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

#Menambah tulisan 1

text(x=29,y=3.3, labels="Titik Potong",cex=0.7)

#Menambah tulisan 2
text(x=3,y=7.3, labels="Data Awal",cex=1)

#PLOT AWAL
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))

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

#Menambah garis 1
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)

#Menambah garis 2
abline(h=mean(y), col="red", lwd=2.5)

#Menambah garis 3
abline(a=2,b=1/10, col="yellow", lwd=2, lty=2)

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

#Menambah tulisan 1

text(x=29,y=3.3, labels="Titik Potong",cex=0.7)

#Menambah tulisan 2
text(x=3,y=7.3, labels="Data Awal",cex=1)

#Menambah tulisan 3
text(x=46,y=7.3, labels="Data Baru",cex=1.2)

LATIHAN 1

plot(sin,-pi,2*pi)

rpois membangkitkan bilangan acak poisson dengan n=100 dan lambda=5 -> fungsi massa peluang poisson dengan lambda=5

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

Latihan 2

a1 <- 1:25
a2 <- rnorm(25,4,2)

plot(a1,a2,pch="w", main="W")

Latihan 3

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

Latihan 4

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

par(mfrow=c(2,2))

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

#plot 2
plot(sin,-pi, 2*pi)

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

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

Latihan 5

Membuat 4 jenis grafik dalam 1 tampilan yang berisi 100 angka acak berdistribusi Normal (5,1) dengan format

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

## [1] 1 2 3 4

screen(3)
boxplot(yb)
title("BoxplotBilangan Acak Normal",cex.main=0.7)

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

screen(2)
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)
plot(xb,yb,pch=16,col=rainbow(100))
title("ScatterPlot Bilangan Acak Normal",cex.main=0.7)

Pembuatan Grafik

Secara umum pembuatan grafik pada package ggplot2 dapat dilakukan dengan menggunakan dua cara yaitu
1. qplot(x,y,geom=“xxx”)

2. ggplot()+geom_xxx(aes(x,y)).
   Sintaks ggplot()+geom_xxx() sebaiknya digunakan saat ingin membuat beberapa jenis grafik (misal scatter plot+line plot) dan saat ingin mengganti parameter-parameter dari plot (misal bentuk titik atau warna titik).

##menggunakan cara 1
set.seed(12)
qplot(x=rnorm(10),y=rnorm(10),geom = "point")

#menggunakan cara 2
ggplot()+geom_point(aes(x=rnorm(10),y=rnorm(10)))

***

Penggantian warna background dapat dilakukan dengan menggunakan fungsi theme_xxx.. Theme yang dapat digunakan antara lain :
* theme_bw()
* theme_dark()
* theme_classic()
* theme_minimal()
* theme_light()
* theme_linedraw()
* theme_gray()

ggplot()+geom_point(aes(x=rnorm(10),y=rnorm(10)))+theme_linedraw()

***

Tanda +

Tanda + sangat berperan penting dalam ggplot2. Tanda ini digunakan untuk menambahkan beberapa jenis plot yang berbeda dalam satu gambar. Selain itu tanda ini juga digunakan untuk mengakses fungsi dalam package ggplot2 yang dapat mengubah parameter grafik. Contoh :

set.seed(12)
ggplot()+geom_line(aes(x=1:10,y=rnorm(10)))+
  #manambah plot titik
  geom_point(aes(x=1:10,y=rnorm(10)))+
  #mengubah nama sumbu x
  xlab("Waktu")+
  #menambahkan nama grafik
  ggtitle("Grafik Percobaan")+
  #mengubah tema
  theme_classic()

Fungsi Aes Fungsi aes merupakan fungsi yang digunakan untuk mengatur bagaimana variabel-variabel dalam data dipetakan dalam plot. Misalnya saja, pada scatterplot variable apa yang menjadi sumbu x dan sumbu y, memberi warna scatterplot berdasarkan kategori tertentu. Berikut contohnya :

set.seed(12)
ggplot()+geom_point(aes(x=rnorm(2),y=rnorm(2),col=c('gen1','gen2')))

argumen col berperan sebagai pewarna titik berdasarkan kategori.

Visualisasi Data

Data yang akan digunakan adalah housing prices yang merupakan dataset bawaan dari R. Berikut merupakan isi dari dataset housing prices.

data_house<-read.csv("https://raw.githubusercontent.com/gerrydito/Sains-Data-S2/master/Praktikum/Visualisasi%20Data/house_price.csv", header = TRUE, sep=",")
head(data_house)

##                  date   price bedrooms bathrooms sqft_living sqft_lot floors
## 1 2014-05-02 00:00:00  313000        3      1.50        1340     7912    1.5
## 2 2014-05-02 00:00:00 2384000        5      2.50        3650     9050    2.0
## 3 2014-05-02 00:00:00  342000        3      2.00        1930    11947    1.0
## 4 2014-05-02 00:00:00  420000        3      2.25        2000     8030    1.0
## 5 2014-05-02 00:00:00  550000        4      2.50        1940    10500    1.0
## 6 2014-05-02 00:00:00  490000        2      1.00         880     6380    1.0
##   waterfront view condition sqft_above sqft_basement yr_built yr_renovated
## 1          0    0         3       1340             0     1955         2005
## 2          0    4         5       3370           280     1921            0
## 3          0    0         4       1930             0     1966            0
## 4          0    0         4       1000          1000     1963            0
## 5          0    0         4       1140           800     1976         1992
## 6          0    0         3        880             0     1938         1994
##                     street      city statezip country
## 1     18810 Densmore Ave N Shoreline WA 98133     USA
## 2          709 W Blaine St   Seattle WA 98119     USA
## 3 26206-26214 143rd Ave SE      Kent WA 98042     USA
## 4          857 170th Pl NE  Bellevue WA 98008     USA
## 5        9105 170th Ave NE   Redmond WA 98052     USA
## 6           522 NE 88th St   Seattle WA 98115     USA

Plot yang biasa digunakan untuk menggambar DISTRIBUSI DATA

Histogram

Misal kita ingin mengetahui distribusi harga rumah dari kolom price maka dapat dengan

ggplot(data_house) +
  geom_histogram(aes(x = price),fill="darkred",col="darkred",binwidth = 5000) +
  ggtitle("Histogram Harga Rumah") +
  ylab("Jumlah Rumah") +
  xlab("Harga Rumah") + 
  theme(plot.title = element_text(hjust = 0.5))

Keterangan :

Boxplot

Biasanya digunakan untuk melihat perbandingan distribusi.

Pada ilustrasi ini kita akan membandingkan distribusi harga rumah dari 10 kota. Oleh karena itu, langkah pertama yang kita lakukan adalah menyaring ke-10 kota tersbut menggukan fungsi filter.

data_house%>%
  filter(city %in% c("Shoreline","Seattle"
                     ,"Auburn","Duvall","Burien",
                     "Tukwila","Vashon","Yarrow Point",
                     "SeaTac","Medina")
         )%>%
  ggplot() +
  geom_boxplot(aes(x = city,y = price,fill=city),show.legend = F) +
  ggtitle("Sebaran Harga Rumah Setiap Kota") +
  ylab("Harga Rumah") +
  xlab("Kota") + 
  theme(plot.title = element_text(hjust = 0.5))+
  theme_bw()

Hal yang perlu diperhatikan :

1.Nama-nama peubah atau kelompok harus ada dalam satu kolom yang bertipe character atau factor.

2. Argumen x di dalam fungsi aes disi dengan nama kolom yang isinya nama-nama peubah atau kelompok tersebut.

3. Jika ingin memberikan warna yang beda-beda pada boxplot maka argumen fill harus diletakan dalam fungsi aes dan diisi dengan nama kolom yang isinya nama-nama peubah atau kelompok tersebut.

Ridgeline

Plot alternative lain untuk membandingkan distribusi data adalah plot Ridgeline, yang mirip seperti histogram. Plot ridgeline bisa dijalankan dalam R dengan menggunakan fungsi geom_density_ridges dari package ggridges. Isi argumenya mirip seperti pada saat pembuatan plot boxplot.

data_house%>%
  filter(city %in% c("Shoreline","Seattle"
                     ,"Auburn","Duvall","Burien",
                     "Tukwila","Vashon","Yarrow Point",
                     "SeaTac","Medina")
         )%>%
  ggplot() +
  geom_density_ridges(aes(y = city,x = price,fill=city), show.legend = F) +
  ggtitle("Sebaran Harga Rumah Setiap Kota") +
  ylab("Kota") +
  xlab("Harga Rumah") + 
  theme(plot.title = element_text(hjust = 0.5))+
  theme_bw()

## Picking joint bandwidth of 133000

***

Plot plot untuk VISUALISASI JUMLAH

Bar Plot

Dalam package ggplot2 bar plot ini bisa ditampilkan menggunakan geom_bar atau geom_col. Perbedaan mendasar dari kedua fungsi itu adalah data yang perlu disediakan.

Fungsi geom_bar hanya membutuhkan satu kolom saja (dengan tipe data character atau factor), sedangkan geom_col membutuhkan dua kolom data yaitu nama kategori dan jumlah setiap kategori.

data_house%>%count(city)%>%
  ggplot()+
  geom_col(aes(x=fct_reorder(as.factor(city),n),y=n), fill="steelblue", 
           width=0.4) +
    scale_y_continuous(expand = c(0,0))+
    coord_flip() +
    ggtitle("Jumlah Rumah Setiap Kota") +
    xlab("") +
    ylab("Jumlah Rumah") +
    theme(plot.title = element_text(hjust = 0.5))+
    theme_classic()

* fct_reorder : untuk mengurutkan kategori berdasarkan nilai tertentu.

• scale_y_continuous(expand = c(0,0)) digunakan agar barplot dimulai dari nol.

• Fungsi coord_flip : untuk menukar sumbu x dan sumbu y.

Lolipop Chart

Plot alternative dari bar chart adalah lolipop chart. Seperti namanya plot ini terinspirasi dari permen lolipop. Untuk membuat plot ini dibutuhkan dua fungsi yaitu geom_segment dan juga geom_point.

Fungsi geom_segment digunakan untuk menggambarkan garis sedangkan fungsi geom_point digunakan untuk menggambarkan titik.

data_house%>%count(city)%>%
  mutate(city=fct_reorder(as.factor(city),desc(n)))%>%
  ggplot()+
  geom_segment(aes(x=city,xend=city, y=0, yend=n), color="skyblue")+
  geom_point(aes(x=city,y=n),color="steelblue", size=2)+ 
    scale_y_continuous(expand = c(0,0))+
    coord_flip() +
    ggtitle("Jumlah Rumah Setiap Kota") +
    xlab("") +
    ylab("Jumlah Rumah") +
    theme(plot.title = element_text(hjust = 0.5))

> Plot untuk Visualisasi Korelasi

Scatter Plot Korelasi atau hubungan dari dua peubah bisa kita visualisasikan menggunakan scatterplot. Jika scatterplot membentuk pola garis maka bisa dikatakan bahwa kedua peubah tersebut memiliki korelasi yang kuat.

ggplot(data_house) +
  geom_point(aes(x = sqft_living,y = price),color="steelblue",size=2) +
  ggtitle("Scatter Plot Harga Rumah vs Luas Rumah") +
  ylab("Luas Rumah") +
  xlab("Harga Rumah") + 
  theme(plot.title = element_text(hjust = 0.5))

***

Correlogram

Penggunaan scatterplot memiliki keterbatasan jika korelasi peubah yang ingin dilihat ada banyak. Coreelogram bisa digunakan untuk mengatasi hal tersebut. Correlogram ini membuat grafik berdasarkan nilai koefisien korelasi yang dikonverssikan dalam bentuk warna

ggcorr(data_house, method = c("everything","pearson"),geom = "tile") 

## Warning in ggcorr(data_house, method = c("everything", "pearson"), geom =
## "tile"): data in column(s) 'date', 'street', 'city', 'statezip', 'country' are
## not numeric and were ignored

Plot Interaktif

Untuk membuat plot interaktif dari ggplot2 kita bisa menggunakan fungsi ggplotly yang berasal dari package plotly.

p1 <- ggplot(data_house) +
  geom_point(aes(x = sqft_living,y = price),color="steelblue",size=2) +
  ggtitle("Scatter Plot Harga Rumah vs Luas Rumah") +
  ylab("Luas Rumah") +
  xlab("Harga Rumah") + 
  theme(plot.title = element_text(hjust = 0.5))
ggplotly(p1)

Invers Transform Method

Fungsi sebaran kumulatif, nilainya 0 sampai 1. Kita bisa dekati dengan sebaran seragam (0,1) karena sebarannya sama.

Kelebihan Bisa diaplikasikan terhadap berbagai sebaran (kontinu dan diskret).

Kesulitan : memperoleh kebalikan dari fungsi sebaran kumulatif.

Algoritma

1. Tentukan fsk F(x) dari sebaran yang diinginkan
2. Tentukan invers F(x)
3. Bangkitkan U~Uniform(0,1)
4. Transformasi X= invers F(u)

Latihan 1

Bangkitkan bilangan acak yang berdistribusi eksponensial(3) dengan Inverse Transform Method, yang amatannya berjumlah
1000. Bandingkan hasilnya dengan fungsi bawaan R rexp dengan menggunakan histogram.

#Eksponensial (lambda=3)

eks <- function(n, lambda){
  U <- runif(n)
  X <- -log(1-U)/lambda
  return(X)
}

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

Agar hasilnya selalu sama setiap menjalankan program, bisa menambahkan fungsi set.seed(). Contoh set.seed(10)

Acceptance-Rejection Method

Membangkitkan (x,y) di daerah tersebut.

Kriteria :

Terima x jika (x,y) di dalam daerah kurva

Tolak x jika (x,y) di luar daerah kurva

Algoritma

1. Bangkitkan x~U(x0,x1)
2. Tentukan c sehingga f(x) <= c untuk x0 < x < x1
3. Bangkitkan y1 ~ U(0,c)
4. Tentukan y2 = f(x)
5. Jika y1 <= y2, terima x

Contoh 1

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

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)

Contoh 2

Bangkitkan peubah acak X ~ Beta(2,2) yang memiliki fungsi PDF f(x)=6x(1−x),0<x<1 dengan menggunakan Acceptance-Rejection Method sebanyak 1000 amatan. Kemudian gambarkan histrogramnya!

Dalam Kasus ini kita pilih Y∼U(0,1) karena memiliki domain yang sama dengan peubah acak X yaitu 0<x<1.

Pdf dari Y adalah g(y)=1,0≤y≤1

# banyaknya amatan
n=1000
set.seed(10)
y = runif(n,0,1)

#Bangkitkan U ~ uniform(0,1)
set.seed(30)
u = runif(n,0,1)

#Mencari nilai c dengan fungsi optimize di R
h <- function(y) 6*y*(1-y)
derivH <- optimize(f=h,interval=c(0,1),maximum=T)
derivH$objective

## [1] 1.5

#Kriteria
f <- function(y) 6*y*(1-y)
g <- function(y) 1
nilaic <- derivH$objective
#kriteria penerimaan
criteria <- u < f(y)/(nilaic*g(y))
head(criteria)

## [1] TRUE TRUE TRUE TRUE TRUE TRUE

#jadikan Y sebagai X
x <- y[criteria]
head(x)

## [1] 0.50747820 0.30676851 0.42690767 0.69310208 0.08513597 0.22543662

# banyaknya amatan
length(x)

## [1] 668

#Perbandingan histogram di R
z <- rbeta(n,shape1 = 2,shape2 = 2)
par(mfrow=c(1,2)) 
hist(x,main="beta(2,2) dari AR")
hist(z,main="beta(2,2) dari rbeta")

Contoh 3

Bangkitkan peubah acak X∼Gamma(32,1) yang memiliki fungsi PDF f(x)=1Γ(32)x12e−x,0<x<∞ dengan menggunakan Acceptance-Rejection Method sebanyak 1000 amatan. Kemudian gambarkan histrogramnya!

Dalam Kasus ini kita pilih Y∼exponential(3/2) karena memiliki domain yang sama dengan peubah acak X yaitu 0<x<∞ dan karena distribusi exponensial merupakan kasus khusus dari distribusi gamma. Pdf dari Y adalah g(y)=132e−y32,0≤y≤∞ atau g(y)=23e−23y,0≤y≤∞

# banyaknya amatan
n=1500
set.seed(10)
y = rexp(n,rate = 2/3)

set.seed(30)
u = runif(n,0,1)


#menghitung nilai c
h <- function(y) (-1) * (3/(2*gamma(3/2))) * (y)^(1/2) * exp(-y/3)
derivH <- optim(par = 0,fn=h,lower = 0,method = "L-BFGS-B")
derivH$value

## [1] -1.257317

#kriteria
f <- function(y) (3/(2*gamma(3/2))) * (y)^(1/2) * exp(-y)
g <- function(y) (2/3) * exp(-(2/3) * y)
nilaic <- -derivH$value #jangan lupa tambahkan negatif
#kriteria penerimaan
criteria <- u < f(y)/(nilaic*g(y))
head(criteria)

## [1] TRUE TRUE TRUE TRUE TRUE TRUE

#jadikan Y sebagai X
x <- y[criteria]
head(x)

## [1] 0.02243461 1.38033181 1.12823841 2.36256278 0.34748792 1.63000950

# banyaknya amatan
length(x)

## [1] 1393

#karena >1000 ambil 1000 saja
x <- x[1:1000]

#menampilkan histogram di R
z <- rgamma(n,shape = 3/2,scale = 1)
par(mfrow=c(1,2)) 
hist(x,main="gamma(3/2,1) dari AR")
hist(z,main="beta(3/2,1) dari rgamma")

Direct Transformation

Memanfaatkan beberapa fungsi transformasi dari berbagai sebaran yang ada. Harus mengetahui hubungan dari sebaran, misal : bernoulli dan binomial, eksponensial dan gamma, gamma dengan chi-square, dsb

Contoh 1

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)

model1 <- lm(Y~x1)
summary(model1)

## 
## Call:
## lm(formula = Y ~ x1)
## 
## 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

model2 <- lm(Y~x2)
summary(model2)

## 
## Call:
## lm(formula = Y ~ x2)
## 
## 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

model3 <- lm(Y~x1+x2)
summary(model3)

## 
## Call:
## lm(formula = Y ~ x1 + x2)
## 
## 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

model4 <- lm(Y~x1:x2)
summary(model4)

## 
## Call:
## lm(formula = Y ~ x1:x2)
## 
## 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

model5 <- lm(Y~x1*x2)
summary(model5)

## 
## Call:
## lm(formula = Y ~ x1 * x2)
## 
## 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

R2 <- matrix(c(summary(model1)$r.squared,
               summary(model1)$adj.r.squared,
               summary(model2)$r.squared,
               summary(model2)$adj.r.squared,
               summary(model3)$r.squared,
               summary(model3)$adj.r.squared,
               summary(model4)$r.squared,
               summary(model4)$adj.r.squared,
               summary(model5)$r.squared,
               summary(model5)$adj.r.squared), 5, byrow=T)

colnames(R2) <- c("R2", "R2.adj"); R2*100

##             R2     R2.adj
## [1,]  4.432592  0.2774876
## [2,] 60.699195 58.9904640
## [3,] 87.090357 85.9167529
## [4,] 68.036265 66.6465370
## [5,] 87.163648 85.3298839

coef(model3)

## (Intercept)          x1          x2 
##  -0.6288095   2.7295126   0.7245911

confint (model3)

##                   2.5 %     97.5 %
## (Intercept) -28.4234482 27.1658292
## x1            1.8854322  3.5735929
## x2            0.5979779  0.8512044

cbind(coef(model3), confint(model3))

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

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

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

Contoh 2

Bangkitkan peubah acak X∼t(2) sebanyak 1000 amatan yang didapatkan dari hasil transformasi langsung peubah acak Y∼χ2(2) dan Z∼N(0,1). Buatlah histogramnya dan bandingkan hasilnya dengan hasil bangkitan yang diperoleh dari fungsi rt.

# Membangkitkan distribusi t(2)

#banyaknya amatan
n <- 1000
set.seed(15)
z <- rnorm(n,mean = 0,sd = 1)
y <- rchisq(n,df = 2)

x <- z/sqrt(y/2)
head(x)

## [1]  0.4249926  2.2981978 -0.3956900  0.8938340  0.9930945 -1.4325087

#menampilkan histogram di R
set.seed(5)
v <- rt(n,df=2)
par(mfrow=c(1,2)) 
hist(x,main="t(2) dari Direct")
hist(v,main="t(2) dari rt")