G14170045 - Laporan Praktikum Pemrograman Statistika 2
Pemrograman Grafik
Misal akan dibuat plot dengan x adalah nilai 1 sampai 25 dan y adalah bilangan acak menyebar poisson dengan nilai lamba 3. Plot dengan type=‘b’ berarti grafik berupa titik dan garis sedangkan untuk type=‘o’ berarti grafik berupa titk dan garis yang saling tumpang tindih dan type=‘S’ berarti grafik berupa tangga dengan segmen pertama vertikal. Seperti yang diperlihatkan pada sintaks berikut:
x <- 1:25
y <- rpois(25,3)
plot(x,y,type="b")
plot(x,y,type="o")
plot(x,y,type="S")
Apabila grafik ingin berupa titik maka gunakan type=‘p’ dengan warna untuk col=1 berarti warna hitam, col=2 berarti warna merah, col=3 berarti warna hijau, dan col=4 berarti warna biru,dst. Selain itu, jika ingin diberi batas pada sumbu x dan y gunakan fungsi xlim dan ylim. Untuk warna raainbow dapat digunakan col=rainbow() dan diisi bilangan pada fungsi rainbow tsb. pch adalah plotting character untuk simbol dari titiknya, dan cex adalah ukuran dari pch. main digunakan untuk menambahkan judul.
plot(x,y,type="p",xlab="Sumbu x",ylab="Sumbu y", main="Bilangan Acak Poisson",col=3,cex=2,pch=16,xlim=c(0,30),ylim=c(0,12))
#menambahkan amatan
x1 <- 26:35
y1 <- rpois(10,3)
points(x1,y1,cex=1)
#menambahkan garis
x2 <- rep(mean(x),20)
y2 <- seq(min(c(y,y1)),13,length=20)
lines(x2,y2,col=4,lwd=2,lty=2)
x3 <- rep(min(x1)-0.5,20)
y3 <- seq(min(c(y,y1)),13,length=20)
lines(x3,y3,col=2,lwd=2,lty=2)
text(x=3,y=8,labels="Data 1 dibawah rata-rata x",cex=0.7)
text(x=20,y=7.3,labels="Data 1 diatas rata-rata x",cex=0.7)
text(x=28,y=7.3,labels="Data 2",cex=0.7)
Jika grafik berupa garis maka dapat digunakan fungsi lty untuk type dari garis dan lwd untuk ketebalan dari garis. sub digunakan untuk menambahkan judul bagian bawah. xlab dan ylab digunakan untuk menamakan masing-masing sumbu x dan sumbu y Seperti pada sintaks berikut:
plot(x, y, type="l",lty=3, lwd=2, sub='plot ke 2', main="Bilangan Acak Poisson",xlab= "X Axis", ylab='Y Axis')
Selanjutnya, untuk membuat histogram yang menyebar normal(3,3) dan bentuk density nya dapat digunakan sintaks berikut:
x3 <- rnorm(100,3,3)
hist(x3,freq=FALSE,ylim=c(0,0.2))
curve(dnorm(x,3,3),col=2,lty=2,lwd=2,add=TRUE)
Apabila ingin menggambar banyak grafik dalam satu device dapat gunakan par(), split.screen() dan layout(). Disini akan dicontohkan dengan par() dan layout()
par(mfrow=c(2,1))
plot(x,y,type="p",xlab="ini sumbu x",ylab="ini sumbu y", main="Bilangan Acak Poisson",col=3,cex=2,pch=16,xlim=c(0,30),ylim=c(0,12))
plot(x, y, type="l",lty=3, lwd=2, main="Bilangan Acak Poisson",xlab= "X Axis", ylab='Y Axis')
par(mfrow=c(1,1))untuk layout() adalah sebagai berikut:
set <- matrix(c
(1,1,3,1,1,3,2,2,4,2,2,4)
, nrow=3)
tempat <- layout(set)
plot(x,y,main='plot x dan y')
boxplot(y1,horizontal =
TRUE,main='boxplot y1')
boxplot(y2,horizontal =
TRUE,main='boxplot y2')
boxplot(x,horizontal =
TRUE,main='boxplot x')
Pembangkitan Bilangan Acak
#pembangkitan bilangan acak normal
x <- rnorm(3) #x~N(0,1) sebanyak 3 bilangan
x1 <- rnorm(10,3,3) #x1~N(3,sd=3) sebanyak 10 bilangan
x2 <- rbinom(10,1,0.3) #x2 menyebar bernoulli(0.3) sebanyak 10 bilangan
x3 <- rbinom(10,2,0.5) #x3 menyebar binomial(2,0.5) sebanyak 10 bilangan
x4 <- rpois(10,3) #x4 menyebar poisson dengan lambda 3 sebanyak 10 bilangan
#mencari nilai peluang kumulatif peubah acak p1 <- pnorm(1.645) #P(Z<1.645)=0.95
p2 <- pnorm(1.64) #P(Z<1.64)=0.95
p3 <- pnorm(-1.64)
#mencari nilai kuantil dari peluang peubah acak
q1 <- qnorm(0.90)
#mencari nilai density peubah acak
##plot density normal
x<- seq(-3,3,length=1000)
y <- dnorm(x)
plot(x,y)
Inverse-Transform
Metode ini menerapkan transformasi integral peluang. Untuk Inverse-Transform Method terdapat kesulitan dan keunggulan. • Kesulitan utama: memperoleh kebalikan dari fungsi sebaran kumulatif • Keunggulan: bisa digunakan untuk berbagai sebaran (termasuk sebaran diskret)
Misal ingin dibangkitkan Eksponensial(9), maka dengan sintaks:
i<-1000
lambda<-9
U<-runif(i)
X<--log(U)/lambda
hist(X)
Misal ingin dibangkitkan Binomial(3,0.33), maka dengan sintaks:
#Inverse-Transform
n=3
i<-1000
p<-.35
X<-runif(i)
Y<-NULL
for (z in 1:i)
ifelse (X[z]<=p,Y[z]<-1,Y[z]<-0)
(tabel<-table(Y)/length(Y)) ## Y
## 0 1
## 0.653 0.347barplot(tabel,main="Bernoulli")
table(Y)## Y
## 0 1
## 653 347misal diinginkan untuk bangkitkan bilangan yang menyebar dengan fungsi f(x)=4x^3, 0<x<1, sehingga F(x) adalah x^4 dan F^-1(u)) adalah u^1/4 maka:
x4 <- function(n){
U <- runif(n)
X <- U^(1/4)
return(X)
}
inverse_transform <- x4(1000) #inverse transform method
hist(inverse_transform,main="Histogram Hasil Metode Inverse Transform")
Acceptance-Rejection Method
Metode ini digunakan dengan cara menolak atau menerima bilangan yang dibangkitkan berdasarkan kriteria tertentu, bilangan yang diterima berarti berada didalam daerah yang diinginkan. Sebagai contoh:
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] 322par(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)
Direct Transformation
Beberapa transformasi sebaran dari sebaran lain dapat digunakan untuk membangkitkan bilangan acak. Misal akan dibangkitkan bilangan acak chisquare dari bilangan acak normal sebagai berikut:
y <- rnorm (1000)
x <- y^2
hist(x,prob=T,main='pembangkitan chisquare')
xx <- seq (0 ,20 ,0.01)
lines(xx,dchisq(xx,df=1),col="red")
catatan untuk direct transformation: • Yi ~ Eksponensial(l) maka untuk membangkitkan gamma G = S Yi ~Gamma(n,l) • Gamma(m/2, 2) = Chi-square(m) • Chi-square(db = 1) = Kuadrat dari Normal(0, 1) • Jika X1~Chi-square(m), X2~Chi-square(n), maka X3=X1+X2~Chi-square(m+n)
Pembangkitan bilangan acak untuk regresi
Misalkan akan dibangun model regresi linear sederhana sebagai berikut: 𝑌=20+5𝑥+𝜖
set.seed(15)
#membangkitkan data
x<-runif(n = 100,min = 0,max = 1)
err<-rnorm(n = 100,mean = 0,sd = 3)
y<-20+5*x+err
data<-as.data.frame(cbind(x,y))
head(data) model_regresi<-function(n=100,b0=20,b1=5,sigma=3){
x<-runif(n = 100,min = 0,max = 1)
err<-rnorm(n,mean = 0,sd = sigma)
y<-b0+b1*x+err
data<-as.data.frame(cbind(x,y))
model<-lm(y~x,data=data) }
#Lalu diulang sebanyak 1000 kali:
sim<-replicate(n=1000,model_regresi(),simplify = F)
library(purrr)
library(broom)
library(dplyr) ##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionsim%>% map_df(tidy)%>% head() #Hasil dari rataan 1000 nilai dugaan parameter tidak berbeda jauh dengan parameternya (β0=20 dan β1=5).