————————– Dalil Limit Pusat ———————————–

Teorema limit pusat: Apapun sebaran populasi X yang memiliki nilai tengah μ dan ragam σ2 , jika diambil sampel secara acak dari populasi tersebut sebanyak n (dimana n cukup besar), maka x akan menyebar mendekati sebaran Normal dengan nilai tengah μ dan ragam σ2 /n

Dalil limit pusat sangat berguna sebagai dasar atau alasan mengapa kita sering menggunakan sebaran normal dalam inferensi statistika walaupun sebaran datanya tidak normal.

Algoritma yang akan digunakan untuk melihat bagaimana dalil limit pusat digunakan adalah sebagai berikut: a) Tentukan ukuran contoh (𝑛) b) Tentukan sebaran populasi data c) Ulang 𝑘 kali • Ambil 𝑛 contoh acak dari sebaran data yang sudah ditentukan • Hitung rataannya lalu simpan d) Periksa sebaran dari 𝑘 rataan

Pertama, akan disimulasikan berdasarkan populasi tak terhingga dengan sebaran normal baku dengan ukuran contoh (n) sebesar 10,30, dan 100

set.seed(1)
library(grDevices)
k<-1000
n<-10
x11<-matrix(rnorm(n*k),k) 
x11<-apply(x11,1,mean)
mean(x11); var(x11)
## [1] -0.006537039
## [1] 0.0979283
k<-1000
n<-30 
x12<-matrix(rnorm(n*k),k) 
x12<-apply(x12,1,mean) 
 mean(x12); var(x12)
## [1] -0.0007272005
## [1] 0.03231146
k<-1000
n<-100
x13<-matrix(rnorm(n*k),k)
x13<-apply(x13,1,mean) 
 mean(x13); var(x13)
## [1] 0.0006768572
## [1] 0.009937682
 par(mfrow=c(1,3))
hist(x11);hist(x12);hist(x13);

Dapat dilihat, nilai harapan dari ketiga plot tersebut adalah sekitar 0 dan ragamnya berbeda-beda untuk setiap ukuran contoh. Pada ukuran contoh 10,30 dan 100, ragamnya berturut-turut sekitar 0.0979283 (mendekati 1/10 = 0.1) , 0.03231146 (mendekati 1/30 = 0.033), dan 0.009937682 (mendekati 1/100 = 0.01) sehingga semakin besar ukuran contohnya makan akan didapatkan histogram yang semakin ramping

Selanjutnya, akan disimulasikan berdasarkan populasi tak terhingga dengan sebaran eksponensial dengan nilai harapan 1 dan ragam 1 dengan ukuran contoh (n) sebesar 10,30, dan 100

set.seed(1)
k<-1000
n<-10 
x21<-matrix(rexp(n*k),k) 
x21<-apply(x21,1,mean) 
mean(x21); var(x21)
## [1] 0.9983612
## [1] 0.09828422
k<-1000
n<-30
x22<-matrix(rexp(n*k),k)
x22<-apply(x22,1,mean)
mean(x22); var(x22)
## [1] 0.9978863
## [1] 0.03344943
k<-1000
n<-100
x23<-matrix(rexp(n*k),k)
x23<-apply(x23,1,mean)
mean(x23); var(x23)
## [1] 1.001066
## [1] 0.009842486
par(mfrow=c(1,3))
hist(x21);hist(x22);hist(x23);

Dapat dilihat, nilai harapan dari ketiga plot tersebut adalah sekitar 1 (walaupun tidak sama persis) dan ragamnya berbeda-beda untuk setiap ukuran contoh. Pada ukuran contoh 10,30 dan 100, ragamnya berturut-turut yaitu 0.09828422 (mendekati 0.1) , 0.03344943 (mendekati 0.033), dan 0.009842486 (mendekati 0.01) sehingga semakin besar ukuran contohnya makan akan didapatkan histogram yang semakin simetris

Selanjutnya, akan disimulasikan berdasarkan populasi tak terhingga dengan sebaran seragam(0,1) dengan nilai harapan 0.5 dan ragam 1/12 dengan ukuran contoh (n) sebesar 10,30, dan 100

set.seed(1)
k<-1000
n<-10 
x31<-matrix(runif(n*k),k) 
x31<-apply(x31,1,mean)
 mean(x31); var(x31)
## [1] 0.500168
## [1] 0.008924686
k<-1000
n<-30
x32<-matrix(runif(n*k),k)
x32<-apply(x32,1,mean) 
mean(x32); var(x32)
## [1] 0.4990147
## [1] 0.002721224
k<-1000
n<-100
x33<-matrix(runif(n*k),k)
x33<-apply(x33,1,mean)
mean(x33); var(x33)
## [1] 0.4992437
## [1] 0.0008488931
par(mfrow=c(1,3))
hist(x31);hist(x32);hist(x33);

Dapat dilihat,nilai harapan dari ketiga plot tersebut adalah sekitar 0.5 (walaupun tidak sama persis) dan ragamnya berbeda-beda untuk setiap ukuran contoh. Pada ukuran contoh 10,30 dan 100, ragamnya berturut-turut yaitu 0.008924686 (mendekati 1/120=0.0833) , 0.002721224 (mendekati 1/360=0.00277), dan 0.0008488931 (mendekati 1/1200=0.000833) sehingga semakin besar ukuran contohnya makan akan didapatkan histogram yang semakin simetris

Agar lebih mudah melihat perbandingan untuk setiap sebaran seragam(0,1), normal baku, dan eksponensial (1), dapat disajikan sebagai berikut:

par(mfrow=c(3,3))
hist(x11);hist(x12);hist(x13);
hist(x21);hist(x22);hist(x23);
hist(x31);hist(x32);hist(x33);

Lalu untuk simulasi dari populasi terhingga dengan dibatasi jumlah populasi yaitu sebanyak 10000, akan disimulasikan berdasarkan sebaran normal baku, eksponensial dengan nilai harapan dan ragam 1 dan seragam(0,1) dengan ukuran contoh (n) sebesar 10,30, dan 100 sebagai berikut:

set.seed(1)
y1<-rnorm(10000) #normal baku
y2<-rexp(10000) #eksponensial nilai harapan dan ragam = 1
y3<-runif(10000) #seragam(0,1)
mean(y1); var(y1) #normal baku
## [1] -0.006537039
## [1] 1.024866
mean(y2); var(y2) #eksponensial nilai harapan dan ragam = 1
## [1] 1.000982
## [1] 0.9771294
mean(y3); var(y3)  #seragam(0,1)
## [1] 0.500708
## [1] 0.08282501
k<-1000
n<-10
z11<-matrix(sample(y1,
                   n*k),k)
z21<-matrix(sample(y2,
                   n*k),k) 
z31<-matrix(sample(y3,
                                              n*k),k) 
z11<-apply(z11,1,mean) 
z21<-apply(z21,1,mean) 
z31<-apply(z31,1,mean) 
mean(z11);var(z11) #normal baku ukuran contoh 10
## [1] -0.006537039
## [1] 0.09933862
mean(z21);var(z21) #eksponensial nilai harapan dan ragam 1 dan ukuran contoh 10
## [1] 1.000982
## [1] 0.09439403
mean(z31);var(z31)  #seragam(0,1) ukuran contoh 10
## [1] 0.500708
## [1] 0.008267329
par(mfrow=c(1,3))
hist(x11);hist(x21);hist(x31);

k<-1000
n<-30
z12<-matrix(sample(y1,
                   n*k,rep=T),k)
z22<-matrix(sample(y2,
                   n*k,rep=T),k) 
z32<-matrix(sample(y3,
                                              n*k,rep=T),k) 
z12<-apply(z12,1,mean) 
z22<-apply(z22,1,mean)
z32<-apply(z32,1,mean)
mean(z12);var(z12) #normal baku ukuran contoh 30
## [1] -0.0006392908
## [1] 0.03037395
mean(z22);var(z22) #eksponensial nilai harapan dan ragam 1 dan ukuran contoh 30
## [1] 1.007364
## [1] 0.03325349
mean(z32);var(z32) #seragam(0,1) ukuran contoh 30
## [1] 0.5007944
## [1] 0.002800989
par(mfrow=c(1,3))
hist(x12);hist(x22);hist(x32);

k<-1000
n<-100
z13<-matrix(sample(y1,
                   n*k,rep=T),k)
z23<-matrix(sample(y2,
                   n*k,rep=T),k) 
z33<-matrix(sample(y3,
                                              n*k,rep=T),k) 
z13<-apply(z13,1,mean) 
z23<-apply(z23,1,mean) 
z33<-apply(z33,1,mean)
mean(z13);var(z13) #normal baku ukuran contoh 100
## [1] -0.004304931
## [1] 0.01102185
mean(z23);var(z23) #eksponensial nilai harapan dan ragam 1 dan ukuran contoh 100
## [1] 0.9953042
## [1] 0.009878623
mean(z33);var(z33) #seragam(0,1) ukuran contoh 100
## [1] 0.5011831
## [1] 0.0008857218
par(mfrow=c(1,3))
hist(x13);hist(x23);hist(x33);

Selanjutnya untuk lebih mudah melihat perbandingannya, dapat disajikan sebagai berikut:

par(mfrow=c(3,3))
hist(z11);hist(z12);hist(z13);
hist(z21);hist(z22);hist(z23);
hist(z31);hist(z32);hist(z33);

Dapat dilihat,kurang lebih sama dengan populasi tak terhingga, karna ulangan yang digunakan pada populasi terhingga besar yaitu k=1000, sehingga semakin besar ukuran contoh, maka histogram yang terbentuk juga semakin simetrik untuk eksponensial nilai harapan dan ragam 1, normal baku maupun sebaran seragam(0,1)

————– Ketakbiasan dan Selang Kepercayaan ——————————–

Pertama akan dilihat ketakbiasan 𝒙 dan s2 sebagai Penduga μ dan σ2 • Algoritme sebagai berikut: – Tentukan sebaran yang akan digunakan, termasuk μ dan σ2 – Ulangi k kali • Bangkitkan n buah data dari sebaran yang ditentukan • Hitung 𝑥 dan s2 – Hitung rata-rata dari 𝑥 dan s2, bandingkan dengan μ dan σ

Misal akan dibangkitkan sebaran normal(10,3), dan akan dilihat ketidakbiasan 𝒙 dan s2 sebagai Penduga μ dan σ2 sebagai berikut

set.seed(1)
n<-10
k<-1000 #ulangan
pop<-rnorm(100,10,sqrt(3)) #populasi terhingga
mean(pop)
## [1] 10.1886
var(pop)*(100-1)/100 #fungsi var di R adalah ragam contoh 
## [1] 2.396083
cth<-matrix(NA,k,n)
for (i in 1:k) cth[i,]<-sample(pop,n) #cth<-matrix(rnorm(n*k,10,sqrt(5)),k) #populasi takhingga
xbar<-apply(cth,1,mean)
dev<-apply(cth,1,var)
mean(xbar)
## [1] 10.17729
mean(dev)
## [1] 2.412825

Dapat dilihat, nilai harapan dan ragam contoh yaitu sebesar berbeda tipis dengan nilai harapan dan ragam populalsi, μ dan σ2 yaitu 10 dan 3.

Selanjutnya akan dilihat Ketakbiasan Penduga Least-Square dalam Regresi Linier, dengan algoritma sebagai berikut: – Tentukan β0 dan β1 – Ulangi k kali • Bangkitkan n buah data X dari sebaran yang tertentu (misal: seragam) • Bangkitkan ε yang menyebar normal dengan nilai tengah 0 dan ragam 𝜎2 • Hitung b0 dan b1 – Hitung rata-rata dari b0 dan b1, bandingkan dengan β0 dan β1

Misal akan dilihat Ketakbiasan Penduga Least-Square dalam Regresi Linier dengan parameter yang ditentukan di awal beta0 yaitu 3 dan beta1 yaitu 20, maka dengan sintaks

set.seed(1)
n<-20 #ukuran contoh
ulangan<-100
beta<-c(3,20) #beta0=3 dan beta1=20 
sigmaerror<-3
betaduga<-matrix(NA,ulangan,2) 
for (i in 1:ulangan){
  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)
}
rataanbetaduga<-apply(betaduga,2,mean)
names(rataanbetaduga)<-c("b0","b1")
rataanbetaduga
##        b0        b1 
##  2.980584 20.000980

Diperoleh beta duga yaitu sebesar berbeda tipis dengan parameter yang ditentukan diawal beta0 sebesar 3 dan beta1 sebesar 20

Selanjutnya untuk selang kepercayaan, seperti yang kita ketahui arti dari SK 95% yaitu SK 95% bagi θ: Kita percaya 95% bahwa selang a sampai b memuat nilai parameter θ yang sebenarnya dan SK 95%: Jika kita melakukan 100 kali percontohan acak dan setiap percontohan acak dibuat SK-nya, maka dari 100 SK yang terbentuk, ada 95 SK yang mencakup parameter, sisanya sebanyak 5 SK meleset tidak mencakup parameter

Maka dengan algoritma berikut:

– Tentukan sebaran yang akan digunakan – Ulangi k kali • Bangkitkan n buah data dari sebaran yang ditentukan • Hitung 𝑥 dan s2 • Hitung 𝜎2 dan buat selang kepercayaan (1-α)% 𝑥 – Hitung proporsi banyaknya selang kepercayaan yang memuat μ, bandingkan dengan (1-α)

Misal akan dibangkitkan selang kepercayaan 95% dari 100 ulangan dengan sebaran yang akan digunakan normal(50,3), maka dengan sintaks:

set.seed(131)
n<-50
ulangan<-100
alpha<-0.05
mu<-50; std<-3
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 #proporsi banyaknya SK yang memuat mu 
## [1] 0.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)

Diperoleh selang kepercayaan memuat myu sebanyak 95 dan 5 yang meleset sesuai dengan SK 95% yang telah ditentukan di awal.

————————————— RESAMPLING ———————————————

Simulasi Monte Carlo yaitu Simulasi yang memanfaatkan informasi mengenai sebaran data yang diketahui (dihipotesiskan, dianggap tahu) dengan pasti. Untuk Pendugaan Peluang misal P(Y > a) atau P(Y < a) maka dengan Y ~ g(x), sehingga X ~ f(x)

Misal X~N(0,1), jika diketahui Y=2X berapa P(Y>2)? diketahui algoritme : 1. Bangkitkan 𝑋~𝑁(0,1) 2. Hitung Y=3X 3. Ulangi langkah 1 dan 2 sebanyak 100000 kali 4. Hitung presentasi nilai 𝑌>2

Maka dengan sintaks:

set.seed(1)
n<-100000
x<-rnorm(n)
y<-2*x
y1<-ifelse(y>2,1,0)
mean(y1) 
## [1] 0.15846
pnorm(2,0,2,lower.tail=F) 
## [1] 0.1586553

diperoleh secara empirik P(Y>2) hampir sama dengan P(Y>2) yang diperoleh dari sebaran hipotetik

Selanjutnya, untuk Simulasi Monte Carlo dalam Pengujian Hipotesis:

Misal peubah acak 𝑋~exp(𝜆), kita ingin menguji hipotesis 𝐻0 : 𝜆 = 2 H1: 𝜆= 3 Jika kita menolak 𝐻0 ketika 𝑥 > 1. Hitung 𝛼 dan kuasa ujinya !

set.seed(1)
lambda<-2
N<-100000
k<-1
x<-rexp(N,lambda)
y<-ifelse(x>k,1,0)
hasil<-mean(y)
hasil#alpha
## [1] 0.13635
lambda<-3
N<-100000
k<-1
x<-rexp(N,lambda) 
y<-ifelse(x>k,1,0)
hasil<-mean(y) #kuasa uji
hasil
## [1] 0.04957

Sehingga, diperoleh hasil kuasa uji adalah sebesar 0.04957 dan alpha sebesar 0.13635

Selanjutnya,misal tersedia data dari sebaran eksponensial: 1.68760 0.03120 0.03037 0.61068 0.91673 0.63169 1.34939 2.99986 0.08164 2.70955

• 𝐻0: 𝜆 = 3

• Apa kesimpulan yang bisa diambil?

dengan algoritma: 1. Hitung tdata 2. Bangkitkan 10 contoh acak~exp(𝜆 =3) 3. Hitung tdist dari sebaran hipotetik 4. Ulangi langkah 2 dan 3 sebanyak 10000 kali 5. Hitung nilai-p=𝑃(𝑡𝑑𝑖𝑠𝑡>𝑡𝑑𝑎𝑡𝑎)

Maka dengan sintaks:

lambda<-3
k<-10000 
data1<-c(1.6876,0.03037,0.91673,1.34939,0.08164,
                  0.0312,0.61068,0.63169,2.99986,2.70955)
n<-length(data1)
m<-mean(data1)
s<-sd(data1) 
tdata<-abs((m-1/lambda)/(s/sqrt(n)))
data2<-matrix(rexp(n*k,lambda),k)
m1<-apply(data2,1,mean)
s1<-apply(data2,1,sd)
tdist<-abs((m1-1/lambda)/(s1/sqrt(n)))
y<-ifelse(tdist>tdata,1,0)
pvalue<-mean(y)
kesimpulan<-ifelse(pvalue<0.05,"Tolak H0","Tak Tolak H0")
pvalue;kesimpulan
## [1] 0.0962
## [1] "Tak Tolak H0"
 #atau bisa dengan cara lain 
y1<-c(tdata,tdist[-k])
y1<-sort(y1)
pvalue1<-1-which(y1==tdata)/k
kesimpulan1<-ifelse(pvalue1<0.05,"Tolak H0",
                    "Tak Tolak H0")
pvalue1;kesimpulan1
## [1] 0.0962
## [1] "Tak Tolak H0"

Maka kesimpulannya adalah Tak Tolak H0 dengan p-value sebesar 0.0962

Menggunakan resampling metode Jacknife: Selanjutnya adalah menduga galat baku dari rasio menggunakan jacknife, dengan algoritma: hitung teta, ambil sample dari x(i) dan y(i) dan hitung penduga jacknife. Sebagai contoh:

store<-LETTERS[1:6]
item1<-c(1363,670,761,746,991,798)
item2<-c(1087,571,518,612,770,655)
data1<-data.frame(store,item1,item2)
n<-nrow(data1)
teta_hat<-mean(data1$item1)/
  (mean(data1$item1)+mean(data1$item2))
y<-matrix(NA,n,n-1)
x<-matrix(NA,n,n-1)
for (i in 1:n) {
  y[i,]<-data1$item1[-i]
  x[i,]<-data1$item2[-i] }
ybar<-apply(y,1,mean)
xbar<-apply(x,1,mean)
teta_i<-ybar/(xbar+ybar)
teta_jk<-n*teta_hat-(n-1)*mean(teta_i) 
var_jk<-(n-1)/n*(sum(teta_i^2)-(n*mean(teta_i)^2)) 
se_jk<-sqrt(var_jk)
teta_jk; var_jk;
## [1] 0.5584068
## [1] 3.852812e-05

diperoleh teta duga sebesar 0.5584068 dan penduga ragam jacknife sebesar 3.852812e-05

Menggunakan resampling metode Bootstrap: dengan algoritma -> Secara berulang bangkitkan ukuran contoh dari y, hitung statistik yang diminati, lalu pelajari perilaku statistik selama pengulangan N, sebagai contoh:

store<-LETTERS[1:6]
item1<-c(1363,670,761,746,991,798)
item2<-c(1087,571,518,612,770,655)
data1<-data.frame(store,item1,item2)
n<-nrow(data1)
b<-1000
y<-matrix(sample(data1$item1,n*b,replace=T),b) 
x<-matrix(sample(data1$item2,n*b,replace=T),b) 
ybar<-apply(y,1,mean)
xbar<-apply(x,1,mean)
teta_i<-ybar/(xbar+ybar)
teta_b<-mean(teta_i)
var_b<-(sum(teta_i^2)-(b*teta_b^2))/(b-1)
se_b<-sqrt(var_b)
teta_b;var_b
## [1] 0.5565476
## [1] 0.001410045

diperoleh teta duga sebesar 0.5565476 dan penduga ragam bootstrap sebesar 0.001410045