Laporan Individu 3 Simulasi Statistika
Annissa Nur Fitria Fathina (G14180031)
12/6/2021
Penerapan Metode Simulasi
Berikut adalah beberapa contoh dari penerapan metode simulasi pada teori-teori statistika.
Teori Limit Pusat
“Rata-rata dari contoh acak yang berasal dari sebaran apapun memiliki sebaran normal jika ukuran contohnya sangat besar”. Sehingga semakin besar n maka rata-rata contoh semakin mendekati sebaran normal dengan nilai tengah miu dan ragam sigma^2/n. Skenario simulasi: membangkitkan populasi geometrik, eksponensial, dan seragam lalu mengambil contoh berukurann 10,30,dan 100, kemudian hitung rata-rata, median, dan ragam contoh. Pengambilan contoh dilakukan sebanyak 1000 kali
Sebaran Populasi Geometrik
Prosedur simulasi: membangkitkan populasi dengan sebaran geometrik (0.1) dengan jumlah populasi tak terhingga, lalu diambil contoh berukuran 10, 30, dan 100 sebanyak 1000 kali ulangan dan hitung rata-rata, median, serta varian dari contoh, lalu buat histogram dan dilihat sebaran contohnya
Simulasi
k <- 1000
n1 <- 10
n2 <- 30
n3 <- 100
x11 <- matrix(rgeom(n1*k, 0.1),k)
x11 <- apply(x11,1,mean)
x12 <- matrix(rgeom(n2*k, 0.1),k)
x12 <- apply(x12,1,mean)
x13 <- matrix(rgeom(n3*k, 0.1),k)
x13 <- apply(x13,1,mean)
ukuran.contoh <- c("n=10","n=30","n=100")
ratarata1 <- c(mean(x11),mean(x12),mean(x13))
median1 <- c(median(x11), median(x12), median(x13))
ragam1 <- c(var(x11),var(x12),var(x13))
hasil1 <- cbind(ukuran.contoh,ratarata1,median1,ragam1)
as.data.frame(hasil1)## ukuran.contoh ratarata1 median1 ragam1
## 1 n=10 9.0494 8.6 8.71743707707708
## 2 n=30 8.95676666666667 8.81666666666667 3.05615036258481
## 3 n=100 9.01564 9.01 0.898393583983984Histogram Sebaran Contoh
par(mfrow=c(3,1))
hist(x11)
hist(x12)
hist(x13)
Sebaran Populasi Eksponensial
Prosedur simulasi: membangkitkan populasi dengan sebaran eksponensial (1) dengan jumlah populasi tak terhingga, lalu diambil contoh berukuran 10, 30, dan 100 sebanyak 1000 kali ulangan dan hitung rata-rata, median, serta varian dari contoh, lalu buat histogram dan dilihat sebaran contohnya
Simulasi
x21 <- matrix(rexp(n1*k),k)
x21 <- apply(x21,1,mean)
x22 <- matrix(rexp(n2*k),k)
x22 <- apply(x22,1,mean)
x23 <- matrix(rexp(n3*k),k)
x23 <- apply(x23,1,mean)
ratarata2 <- c(mean(x21),mean(x22),mean(x23))
median2 <- c(median(x21), median(x22), median(x23))
ragam2 <- c(var(x21),var(x22),var(x23))
hasil2 <- cbind(ukuran.contoh,ratarata2,median2,ragam2)
as.data.frame(hasil2)## ukuran.contoh ratarata2 median2 ragam2
## 1 n=10 1.00346150767869 0.977232425268395 0.0943444882941472
## 2 n=30 1.00443098425896 0.991274617706445 0.031982274217917
## 3 n=100 0.998307172296098 0.998402171316974 0.00975601667512815Histogram Sebaran Contoh
par(mfrow=c(3,1))
hist(x21)
hist(x22)
hist(x23)
Sebaran Populasi Seragam
Prosedur simulasi: membangkitkan populasi dengan sebaran seragam (0,1) dengan jumlah populasi tak terhingga, lalu diambil contoh berukuran 10, 30, dan 100 sebanyak 1000 kali ulangan dan hitung rata-rata, median, serta varian dari contoh, lalu buat histogram dan dilihat sebaran contohnya
Simulasi
x31 <- matrix(runif(n1*k),k)
x31 <- apply(x31,1,mean)
x32 <- matrix(runif(n2*k),k)
x32 <- apply(x32,1,mean)
x33 <- matrix(runif(n3*k),k)
x33 <- apply(x33,1,mean)
ratarata3 <- c(mean(x31),mean(x32),mean(x33))
median3 <- c(median(x31), median(x32), median(x33))
ragam3 <- c(var(x31),var(x32),var(x33))
hasil3 <- cbind(ukuran.contoh,ratarata3,median3,ragam3)
as.data.frame(hasil3)## ukuran.contoh ratarata3 median3 ragam3
## 1 n=10 0.500000642781169 0.49894976192154 0.00840068407309818
## 2 n=30 0.499831952581593 0.501293331582565 0.00277544315529291
## 3 n=100 0.500374405886808 0.499202231038362 0.000838587420821748Histogram Sebaran Contoh
par(mfrow=c(3,1))
hist(x31)
hist(x32)
hist(x33)
Kesimpulan Simulasi Teorema Limit Pusat
Berdasarkan output diatas, terlihat bahwa semakin besar ukuran contoh, nilai mean dan median nya memiliki nilai yang relatif sama, ini menunjukkan bentuk kurva yang semakin simetris. Dilihat dari histogramnya pun, semakin besar ukuran contoh, kurva mendekati sebaran normal.
Ketakbiasan Penduga Rata-rata dan Ragam
Data memiliki ukuran pemusatan diantaranya rata-rata dan median serta ukuran penyebaran diantaranya ragam. Berikut kita akan menguji ketakbiasan penduga ukuran pemusatan dan penyebaran dengan simulasi.
Simulasi
Skenario simulasi: Membangkitkan populasi dari sebaran normal baku (0,1) sejumlah 20 angka. Selanjutnya hitung parameter populasi berupa rata-rata, median, dan ragam. Lalu diambil contoh sebesar 4 dan dicari nilai statistik contoh berupa rata-rata, median, ragam (pembagi: n-1), dan ragam (pembagi: n) dari semua kemungkinan contohnya dan bandingkan nilai harapan statistik contoh dengan parameter populasinya.
Populasi
set.seed(1234)
populasi <- rnorm(20)
#Rata2 Populasi
mean <- mean(populasi)
#Ragam Populasi
varians = function(populasi)
{n = length(populasi)
sum = 0
summ = 0
for(i in 1:n)
{sum = sum+populasi[i]}
rata = sum/n
for(i in 1:n)
{sel = populasi[i]-rata
pang=sel^2
summ=summ+pang}
var=summ/n
print(var)}
sigma <- varians(populasi)## [1] 0.9764156#Median Populasi
median <- median(populasi)
parameter <- c(mean, median, sigma, sigma)Contoh
#CONTOH
choose(20,4)## [1] 4845library(prob)## Warning: package 'prob' was built under R version 4.0.5## Loading required package: combinat##
## Attaching package: 'combinat'## The following object is masked from 'package:utils':
##
## combn## Loading required package: fAsianOptions## Warning: package 'fAsianOptions' was built under R version 4.0.5## Loading required package: timeDate## Loading required package: timeSeries## Warning: package 'timeSeries' was built under R version 4.0.5## Loading required package: fBasics## Warning: package 'fBasics' was built under R version 4.0.5## Loading required package: fOptions## Warning: package 'fOptions' was built under R version 4.0.5##
## Attaching package: 'prob'## The following objects are masked from 'package:base':
##
## intersect, setdiff, unioncontoh <- urnsamples(populasi, size = 4, replace = F, ordered = F)
#1. Menghitung rata-rata contoh
ratarata <- rowMeans(contoh)
ratarata <- format(ratarata, digits = 4)
ratarata <- as.numeric(ratarata)
xbar <- mean(ratarata)
#2. Menghitung ragam contoh
standev <- rowStdevs(contoh)
standev <- format(standev, digits = 4)
standev <- as.numeric(standev)
#Ragam contoh dengan pembagi (n-1)
ragam.n1 <- standev^2
ragam_n1 <- mean(ragam.n1)
#Ragam contoh dengan pembagi (n)
ragam.n <- ragam.n1*3/4
ragam_n <- mean(ragam.n)
#3. Menghitung median contoh
library(psych)## Warning: package 'psych' was built under R version 4.0.5##
## Attaching package: 'psych'## The following object is masked from 'package:prob':
##
## sim## The following object is masked from 'package:fBasics':
##
## tr## The following object is masked from 'package:timeSeries':
##
## outlierringkasan <- describe(contoh)
median.s <- mean(ringkasan$median)
statistik <- c(xbar, median.s, ragam_n1, ragam_n)Perbandingan parameter dan statistik
ukuran <- c("rata-rata","median","ragam(/n-1)","ragam(/n)")
perbandingan <- cbind(ukuran, parameter, statistik)
as.data.frame(perbandingan)## ukuran parameter statistik
## 1 rata-rata -0.250664058895984 -0.250663764336429
## 2 median -0.528820680795415 -0.520366479282141
## 3 ragam(/n-1) 0.976415551662844 1.02780569486431
## 4 ragam(/n) 0.976415551662844 0.770854271148235Kesimpulan Simulasi Ketakbiasan Penduga Rata-rata dan Ragam
Nilai harapan dari rata-rata contoh = rata-rata populasi, sehingga xbar merupakan penduga yang tak bias bagi rata-rata populasi. Nilai harapan dari median contoh tidak sama dengan median populasi, sehingga median merupakan penduga yang berbias bagi nilai tengah populasi. Nilai harapan dari ragam contoh (pembagi n-1) mendekati ragam populasi, sehingga ragam contoh (pembagi n-1) merupakan penduga tak bias bagi ragam populasi. Nilai harapan dari ragam contoh (pembagi n) tidak mendekati ragam populasi, sehingga ragam contoh (pembagi n) merupakan penduga yang berbias bagi ragam populasi.
Ketakbiasan Penduga Koefisien Regresi dengan Teknik Least Square
Penduga koefisien regresi dikatakan tidak bias apabila nilai harapan dari beta duga = beta (nilai parameternya).
Simulasi
Skenario simulasi:menentukan beta0=8 dan beta1=16.Lalu membangkitkan 30 data dari sebaran seragam (0,1) sebagai peubah x. Lalu membangkitkan epsilon dari sebaran normal (0,4). Lalu hitung peubah y yang merupakan fungsi dari x, beta, dan epsilon, yaitu y=X*beta + epsilon.Lalu menghitung b0 dan b1 (b: penduga beta). Pembangkitan data hingga pendugaan beta dilakukan berulang sebanyak 1000 kali. Selanjutnya dihitung nilai harapan dari b0 dan b1, serta bandingkan nilainya dengan beta0 dan beta1.
n<-30
k<-1000
sigmaerror<-4
beta<-c(8,16)
betaduga<-matrix(NA,k,2)
for (i in 1:k){
x<-runif(n)*10
epsilon<- rnorm(n,0,sqrt(sigmaerror))
X<-cbind(1,x)
y<-X %*% beta + epsilon
betaduga[i,]<-solve(t(X) %*% X) %*% (t(X) %*% y)
}
beta.duga <- apply(betaduga,2,mean)
beta.duga <- round(beta.duga)
koefisien <- c("beta0","beta1")
perbandingan.beta <- cbind(koefisien, beta, beta.duga)
as.data.frame(perbandingan.beta)## koefisien beta beta.duga
## 1 beta0 8 8
## 2 beta1 16 16Nilai harapan dari beta duga = beta, sehingga beta duga merupakan penduga tak bias bagi beta (parameter regresi).
Selang Kepercayaan
Selang kepercayaan 95% bagi parameter berarti kita percaya 95% bahwa selang a hingga b mengandung nilai parameter yang sebenernya.
Simulasi
Skenario simulasi:menentukan sebaran yang akan digunakan, yaitu sebaran normal (40,16), serta menentukan tingkat kepercayaan yaitu 95% sehingga alpha=5%. Lalu bangkitkan data sebanyak 40 buah dari sebaran normal(40,16), serta hitung xbar, s^2, standar error, dan buat selang kepercayaan. Pembangkitan data hingga pembuatan selang kepercayaan diulang sebanyak 100 kali. Lalu dihitung proporsi selang kepercayaan yang memuat mu = 40, bandingkan dengan 95% (tingkat kepercayaannya).
n<-40
ulangan<-100
alpha<-0.05
mu<-40
std<-4
set.seed(124)
sampel<-matrix(rnorm(n*ulangan,mu,std),ulangan)
xbar<-apply(sampel,1,mean)
dev<-apply(sampel,1,sd)
SE<-dev/sqrt(n)
z<-qnorm(1-alpha/2)
sk<-(xbar-z*SE<mu & mu<xbar+z*SE)
sum(sk)/ulangan## [1] 0.95Proporsi yang didapatkan sama dengan tingkat kepercayaan yang telah ditentukan di awal. Sehingga SK 95% dapat diartikan bila kita melakukan 100 kali percontohan acak dan dibuat SKnya, maka dari 100 SK yang terbentuk, ada 95 SK yang memuat parameter (mu=40), dan 5 SK sisanya kemungkinan meleset atau tidak mencakup parameter. Berikut adalah penggambaran dari SK 95%.
matplot(rbind(xbar-z*SE,xbar+z*SE),
rbind(1:ulangan,1:ulangan),
col=ifelse(sk,"blue","red"),
type="l",lty=1)
abline(v=mu)
Dari output tersebut terlihat bahwa terdapat 5 selang yang tidak mencakup nilai parameter (mu=40) diberi warna merah, sementara 95 selang mencakup nilai parameter diberi warna biru.
Simulasi Monte Carlo untuk Pengujian Hipotesis Tentang Parameter
Skenario simulasi: Pengujian hipotesis tentang parameter dibatasi yaitu nilai tengah dan proporsi dari 1 populasi serta beda nilai tengah dan beda proporsi dari 2 populasi pada sebaran tidak normal yaitu sebaran poisson dengan jumlah populasi tak terhingga
Pengujian Hipotesis Nilai Tengah
Prosedur simulasi: mengambil 16 contoh acak yang menyebar poisson(4), lalu hitung statistik uji t. Selanjutnya membangkitkan contoh acak berukuran 16 dari sebaran poisson(4), hitung nilai t sebagai tdist, ulangi pengambilan contoh dan hitung tdist sebanyak 10000 kali. lalu hitung proporsi nilai tdist > thit, nilai ini disebut p-value
Hipotesis:
H0: lambda = 4 vs H1: lambda != 4 (lambda tidak sama dengan 4)
Statistik uji
n = 16
lambda = 4
set.seed(1234)
contoh <- rpois(n,lambda)
xbar <- mean(contoh)
stdev <- sd(contoh)
thit <- abs((xbar-4)/(stdev/sqrt(n)))
xbar;stdev;thit ## [1] 4## [1] 1.75119## [1] 0Titik kritis (p-value)
k = 10000
set.seed(1234)
distribusi <- matrix(rpois(n*k,lambda),k)
m <- apply(distribusi,1,mean)
s <- apply(distribusi,1,sd)
tdist <- abs((m-4)/(s/sqrt(n)))
y <- ifelse(tdist>thit,1,0)
pvalue <- mean(y)
pvalue## [1] 0.9518Kriteria Penolakan
H0 ditolak jika P-value < alpha. alpha yang digunakan sebesar 5%.
Kesimpulan
kesimpulan <- ifelse(pvalue<0.05, "Tolak H0", "Tak Tolak H0")
kesimpulan## [1] "Tak Tolak H0"Kesimpulannya cukup bukti untuk menyatakan nilai tengah populasi sama dengan 4 pada taraf nyata 5%.
Pengujian Hipotesis Beda Nilai Tengah 2 Populasi
Prosedur simulasi: mengambil 16 contoh acak yang menyebar poisson(4) dan poisson(8) dengan set.seed yang berbeda, lalu hitung statistik uji t. Selanjutnya membangkitkan contoh acak berukuran 16 dari sebaran poisson(4) dan poisson(8) dengan set.seed yang sama dengan contoh, hitung nilai t sebagai tdist dengan asumsi ragam populasi sama, ulangi pengambilan contoh dan hitung tdist sebanyak 10000 kali. lalu hitung proporsi nilai tdist > thit, nilai ini disebut p-value
Hipotesis:
H0: lambda1 = lambda2 vs H1: lambda1 != lambda2 (lambda1 tidak sama dengan lambda2)
Statistik uji
n = 16
lambda1 = 4
lambda2 = 5
set.seed(1234)
contoh1 <- rpois(n,lambda1)
xbar1 <- mean(contoh1)
stdev1 <- sd(contoh1)
set.seed(4567)
contoh2 <- rpois(n,lambda2)
xbar2 <- mean(contoh2)
stdev2 <- sd(contoh2)
sgab <- sqrt(((n-1)*stdev1^2+(n-1)*stdev2^2)/(n+n-2))
thit <- abs((xbar1-xbar2)/(sgab*sqrt(1/n+1/n)))
thit ## [1] 0.9745586Titik kritis (p-value)
k = 10000
set.seed(1234)
distribusi1 <- matrix(rpois(n*k,lambda1),k)
m1 <- apply(distribusi1,1,mean)
s1 <- apply(distribusi1,1,sd)
set.seed(4567)
distribusi2 <- matrix(rpois(n*k,lambda2),k)
m2 <- apply(distribusi2,1,mean)
s2 <- apply(distribusi2,1,sd)
sgab <- sqrt(((n-1)*s1^2+(n-1)*s2^2)/(n+n-2))
tdist <- abs((m1-m2)/(sgab*sqrt(1/n+1/n)))
y <- ifelse(tdist>thit,1,0)
pvalue <- mean(y)
pvalue## [1] 0.6607Kriteria Penolakan
H0 ditolak jika P-value < alpha. alpha yang digunakan sebesar 5%.
Kesimpulan
kesimpulan <- ifelse(pvalue<0.05, "Tolak H0", "Tak Tolak H0")
kesimpulan## [1] "Tak Tolak H0"Kesimpulannya tidak cukup bukti untuk menyatakan nilai tengah kedua populasi berbeda pada taraf nyata 5%.
Pengujian Hipotesis Proporsi untuk Contoh Berukuran Kecil
Prosedur simulasi: mengambil 16 contoh acak yang menyebar poisson(4), bila bilangan acak > 1 maka diberi label = 1, lalu hitung statistik uji zhit. Selanjutnya membangkitkan contoh acak berukuran 16 dari sebaran binomial(0.4), hitung nilai z sebagai zdist, ulangi pengambilan contoh dan hitung zdist sebanyak 10000 kali. lalu hitung proporsi nilai zdist > zhit, nilai ini disebut p-value
Hipotesis:
H0: P = 0.4 vs H1: P!= 0.4 (proporsi populasi tidak sama dengan 0.4)
Statistik uji
n = 16
lambda = 4
p0 = 0.4
set.seed(1234)
contoh <- rpois(n,lambda)
contoh <- ifelse(contoh>1,1,0)
pduga <- sum(contoh)/n
zhit <- abs((pduga-p0)/sqrt(p0*(1-p0)/n))
pduga;zhit## [1] 0.9375## [1] 4.388669Titik kritis (p-value)
k = 10000
set.seed(1234)
distribusi <- matrix(rbinom(n*k,1,p0),n,k)
p <- apply(distribusi,2,sum)/n
zdist <- abs((p-p0)/sqrt(p0*(1-p0)/n))
y <- ifelse(zdist>zhit,1,0)
pvalue <- mean(y)
pvalue## [1] 0Kriteria Penolakan
H0 ditolak jika P-value < alpha. alpha yang digunakan sebesar 5%.
Kesimpulan
kesimpulan <- ifelse(pvalue<0.05, "Tolak H0", "Tak Tolak H0")
kesimpulan## [1] "Tolak H0"Kesimpulannya cukup bukti untuk menyatakan proporsi populasi tidak sama dengan 0.4 pada taraf nyata 5%.