Tugas Praktikum Kedua STA561 Pemrograman Statistika
Klik disini untuk ke halaman rpubs.
Plot (Praktikum Pertemuan Keempat)
Plot:
Generic function for plotting of R objects. General Form:
plot(x, y, .)
Possible type to be drawn:
"p" for points
"l" for lines
"b" for both
"c" for the lines part alone of "b"
"o" for overplotted
"h" for histogram like
"s" for stair steps
"S" for other steps
"n" for no plotting
Pilihan Warna Grafik:
Numeric 1-8 (recycle)
Character "colorname" (colors())
Fungsi rainbow(.), rgb(.)
Package "colorspace", "colourpicker", "RColorBrewer"
dll
Plot Dasar
#Points
x <- 1:40
y <- rnorm(40,5,1)
plot(x,y,type="p")
#Lines
x <- 1:40
y <- rnorm(40,5,1)
plot(x,y,type="o")
#tidak ada plot
x <- 1:40
y <- rnorm(40,5,1)
plot(x,y,type="n")
#Plot Points (Titik) dengan nama judul grafik, nama axis horizontal dan vertikal serta warna titik merah dengan karakter plotting nya no 16
x <- 1:40
y <- rnorm(40,5,1)
plot(x,y,type="p",xlab="Sumbu x",ylab="Sumbu y",main="Bilangan Acak Normal",col=2,pch=16)
#Plot Points (Titik) dengan nama judul grafik, nama axis horizontal dan vertikal serta warna titik rainbow dengan karakter plotting nya no 16. Sumbu X limit sampai dengan angka 50
x <- 1:40
y <- rnorm(40,5,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))
#menambahkan amatan
x1 <- 41:50
y1 <- rnorm(10,5,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))
points(x1,y1,cex=2)
#menambahkan garis
x <- 1:40
y <- rnorm(40,5,1)
x1 <- 41:50
y1 <- rnorm(10,5,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))
points(x1,y1,cex=2)
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
x <- 1:40
y <- rnorm(40,5,1)
x1 <- 41:50
y1 <- rnorm(10,5,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))
points(x1,y1,cex=2)
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
x <- 1:40
y <- rnorm(40,5,1)
x1 <- 41:50
y1 <- rnorm(10,5,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))
points(x1,y1,cex=2)
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)
arrows(x0=30,y0=3.5,x1=40,y1=mean(y)-.1,lwd=2)
#menambahkan tulisan
x <- 1:40
y <- rnorm(40,5,1)
x1 <- 41:50
y1 <- rnorm(10,5,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))
points(x1,y1,cex=2)
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)
arrows(x0=30,y0=3.5,x1=40,y1=mean(y)-.1,lwd=2)
text(x=29,y=3.3,labels="Titik potong",cex=0.7)
#menambahkan tulisan
x <- 1:40
y <- rnorm(40,5,1)
x1 <- 41:50
y1 <- rnorm(10,5,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))
points(x1,y1,cex=2)
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)
arrows(x0=30,y0=3.5,x1=40,y1=mean(y)-.1,lwd=2)
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
Could you explain what are these programs do for ?
plot(sin,-pi, 2*pi)
plot(table(rpois(100,5)),type="h",col="red",
lwd=1,main="rpois(100,lambda=5)")
Latihan 2
Create some programs to make the scatterplot with "w" as point:
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
#Create some programs to make the histogram, using 100 observation of X~cHI-sQUARE df 4
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)
par(mfrow=c(2,2))
plot(1:40,y,type="p",xlab="Sumbu x",ylab="Sumbu y",
main="Bilangan Acak Normal",col=2,pch=16)
plot(sin,-pi, 2*pi)
plot(table(rpois(100,5)),type="h",col="red",
lwd=1,main="rpois(100,lambda=5)")
plot(a1,a2,type="n",main="W")
text(a1,a2,labels=paste("w",1:25,sep=""),
col=rainbow(25),cex=0.8)
par(mfcol=c(2,2))
plot(1:40,y,type="p",xlab="Sumbu x",ylab="Sumbu y",
main="Bilangan Acak Normal",col=2,pch=16)
plot(sin,-pi, 2*pi)
plot(table(rpois(100,5)),type="h",col="red",lwd=1,
main="rpois(100,lambda=5)")
plot(a1,a2,type="n",main="W")
text(a1,a2,labels=paste("w",1:25,sep=""),
col=rainbow(25),cex=0.8)
Latihan 5
# Create 4 graph in one window from 100 random numbers which follow N(5,1). windows()
yb <- rnorm(100,5,1)
split.screen(c(2,2))## [1] 1 2 3 4screen(3)
boxplot(yb)
title("Boxplot Bilangan Acak Normal",cex.main=0.7)
screen(4)
xb <- 1:100
plot(xb,yb,type="l",lwd=2,col="blue")
title("Line Plot 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("Scatter Plot Bilangan Acak Normal",cex.main=0.7)
Pembangkitan Bilangan Acak (Praktikum Pertemuan Kelima)
Pembangkitan bil. acak merupakan alat yang diperlukan dalam komputasi statistik -> umumnya untuk simulasi. Bilangan acak yang dibangkitkan merupakan pseudorandom (acak semu) dan memenuhi sebaran statistik tertentu (pdf/pmf, cdf). Pada umumnya, semua metode pembangkitkan bil. acak tergantung dari pembangkitan bilangan acak uniform
Fungsi umumnya dimulai dengan huruf r diikuti dengan nama sebaran. Contoh:
Data sebaran pseudo normal: rnorm()
Data sebaran pseudo seragam: runif()
Fungsi density/mass (pdf/pmf) dimulai 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
#pembangkitan bilangan acak
x <- rnorm(10) #x~N(0,1)
x## [1] 0.87461571 -1.63715220 1.10821313 -1.38597386 0.76267354 0.11406511
## [7] 1.09850713 -0.09149863 0.87311930 0.21560128#pembangkitan bilangan acak
x1 <- rnorm(10,3,2) #x1~N(3,sd=2)
x1## [1] 4.63272039 0.02745841 0.25981921 2.94196937 6.47390326 3.48023029
## [7] 3.12626712 5.60120548 6.94705755 4.51634334#pembangkitan bilangan acak
x2 <- rbinom(10,1,0.4) #x2~bernoulli(0.4)
x2## [1] 1 0 0 1 0 1 0 0 0 0#mencari nilai peluang kumulatif peubah acak
p1 <- pnorm(1.645) #P(Z<1.645)=0.95
p1## [1] 0.9500151#mencari nilai peluang kumulatif peubah acak
p2 <- pnorm(1.96) #P(Z<1.975)=0.975
p2## [1] 0.9750021#mencari nilai peluang kumulatif peubah acak
p3 <- pnorm(-1.96)
p3## [1] 0.0249979#mencari nilai peluang kumulatif peubah acak
p4 <- pf(15,df1=10,df2=15)
p4## [1] 0.9999955#mencari nilai kuantil dari peluang peubah acak
q1 <- qnorm(0.975)
q1## [1] 1.959964#mencari nilai kuantil dari peluang peubah acak
q2 <- qnorm(0.95,2,1) #X~N(2,1), P(X<x)=0.95
q2## [1] 3.644854#mencari nilai density peubah acak
##plot density normal
a <- seq(-4,4,length=1000)
da <- dnorm(a)
plot(a,da)
tiga teknik umum dalam melakukan pembangkitan bilangan acak, yaitu:
Invers Transform Method
Acceptance-Rejection Method
Direct Transformation
Invers Transform Method
#Latihan 1
#Eksponensial(lambda=3)
eks <- function(n,lambda){
U <- runif(n)
X <- -log(1-U)/lambda
return(X)
}
yy1 <- eks(1000,3) #inverse transform method
yy2 <- rexp(1000,rate=3) #fungsi bawaan R
par(mfrow=c(1,2))
hist(yy1,main="Exp dari Inverse Transform")
hist(yy2,main="Exp dari fungsi rexp")
Acceptance-Rejection Method
#Latihan 2
#Bangkitkan bilangan acak yang memiliki fkp f(y)=3y^2; 0<y<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] 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
Memanfaatkan beberapa fungsi transformasi dari berbagai sebaran yang ada
#Jika X~Uniform(0,1) menggunakan Y=(b-a)X+a,maka Y menyebar Uniform (a,b)
#Membangkitkan bilangan acak Y = 2X + 3
x <- runif(100)
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.813428Pembangkitan Bilangan Acak untuk Model Analisis Regresi
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.3124set.seed(123)
X1 <- runif(25,10,25)
X2 <- runif(25,90,200)
Y <- 10 + 2.3*X1 + 0.7*X2 + rnorm(25,0,9)
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-06set.seed(123)
X1 <- runif(25,10,25)
X2 <- runif(25,90,200)
Y <- 10 + 2.3*X1 + 0.7*X2 + rnorm(25,0,9)
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-10set.seed(123)
X1 <- runif(25,10,25)
X2 <- runif(25,90,200)
Y <- 10 + 2.3*X1 + 0.7*X2 + rnorm(25,0,9)
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-07set.seed(123)
X1 <- runif(25,10,25)
X2 <- runif(25,90,200)
Y <- 10 + 2.3*X1 + 0.7*X2 + rnorm(25,0,9)
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-09R2 <- 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.3298839coef(model3)## (Intercept) X1 X2
## -0.6288095 2.7295126 0.7245911confint(model3)## 2.5 % 97.5 %
## (Intercept) -28.4234482 27.1658292
## X1 1.8854322 3.5735929
## X2 0.5979779 0.8512044cbind(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.8512044anova(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 ' ' 1par(mfrow=c(2,2)); plot(model3)
Eksplorasi Mandiri
Mencoba membangkitkan bilangan acak dari populasi normal serta simulasi penarikan sampel, dan membuat plot nya.
#Package
library(dplyr)
library(ggplot2)
library(gganimate)
library(ggpubr)
library(gridExtra)
library(reshape2)#Generate normal population dengan N = 10000 sd = 5 mu = 15
populasinormal =rnorm(10000,15,5)#Fungsi untuk menghasilkan plot untuk berbagai macam n sample
simulate_density = function(populasi , n){
sample_mean = replicate(100,sample(populasinormal,n)) %>% colMeans() %>% as.data.frame()
colnames(sample_mean) = "xbar"
plot_mnorm = ggplot(sample_mean, aes(x = xbar)) +
geom_density(col = "black",fill = "#56B4E9")+ ylim(c(0,.7))+
labs(title = paste("Density Plot Sample Mean dengan n =",n))+
theme(
plot.title = element_text(color = "black", size = 10, face = "bold")
)
plot_mnorm
}#Generate Density Plot dengan fungsi yang telah dibentuk
hist_5 = simulate_density(populasinormal,5)
hist_20 = simulate_density(populasinormal,20)
hist_50 = simulate_density(populasinormal,50)
hist_70 = simulate_density(populasinormal,70)
ggarrange(hist_5,hist_20,hist_50,hist_70,nrow = 2,ncol = 2)
#Generate Density Plot dengan fungsi yang telah dibentuk
hist_5 = simulate_density(populasinormal,5)
hist_20 = simulate_density(populasinormal,20)
hist_50 = simulate_density(populasinormal,50)
hist_70 = simulate_density(populasinormal,70)
ggarrange(hist_5,hist_20,hist_50,hist_70,nrow = 2,ncol = 2)
#Pembentukan sample mean untuk setiap n
for(n in 2:30) {
x = replicate(100,sample(populasinormal,n)) %>% colMeans()
if(n==2) sample_mean = c(dataset = n, xbar =x )
else {sample_mean = rbind(sample_mean,c(dataset = n, xbar =x ))}
}
sample_mean = as.data.frame(sample_mean)
sample_mean = sample_mean %>% melt(id=1) %>% arrange(dataset)#Plot animasi perubuhan density
dens = sample_mean %>% ggplot(aes(value)) +geom_density(alpha = 0.85, color = "#000000",fill = "#FF3838", size=1) +
transition_states(states = dataset,transition_length = 48 ,state_length = 2) +
ease_aes('cubic-in-out') +
labs(title = 'Plot Density Sample Mean dari Normal Distribution dengan n = {closest_state}')+
theme(plot.title = element_text(color = "black", size = 12, face = "bold"))
animate(dens,nframes = 175,end_pause=40)
Kesimpulan:
Dari simulasi diatas dapat disimpulkan secara empirik bahwa untuk populasi yang memiliki distribusi normal, sample mean akan menyebar normal untuk n berapapun.
TERIMA KASIH
Tugas Praktikum 2 STA561 Pemrograman Statistika, Reni Amelia, Mahasiswa Pascasarjana Statistika dan Sains Data, IPB University, reniamelia@apps.ipb.ac.id↩