1. Dalil Limit Pusat

• Untuk sembarang peubah acak Xi yang saling bebas dan identik, Y= \(\sum_{i=1}^{n}\) Xi akan menyebar normal untuk n cukup besar.

• Untuk Ui ~ Uniform (0,1) \(\rightarrow\) X=\(\sum_{i=1}^{n}\) Ui ~ Normal (\(\mu,\sigma ^{2}\))

  • E(Ui) = \(\frac{1}{2}\) dan Var(Ui) = \(\frac{1}{12}\) , maka X=\(\sum_{i=1}^{12}\) Ui - 6 ~ Normal (0,1)

Pembuktian :

E(X) = E ( \(\sum_{i=1}^{12}\) Ui - 6 ) = E (U1) + E (U2) + ... + E (U12) - 6 = \(\frac{1}{2}\) + \(\frac{1}{2}\) + ... + \(\frac{1}{2}\) - 6 = 12 (\(\frac{1}{2}\)) - 6 = 0

Var(X) = V ( \(\sum_{i=1}^{12}\) Ui - 6 ) = Var (U1) + E (U2) + ... + E (U12) - 0 = \(\frac{1}{12}\) + \(\frac{1}{12}\) + ... + \(\frac{1}{12}\) = 12 (\(\frac{1}{12}\)) = 1

Bagaimana jika kita ingin membangkitkan peubah acak Normal(5,2)?

yaitu dengan cara menentukan nilai n dan k untuk persamaan X = \(\sum_{i=1}^{n}\) Ui - k ; X ~ Normal (5,2)

Cara memperoleh n dan k gunakan rumus berikut :

1. n = 12\(\sigma ^{2}\)

2. k = \(\frac{1}{2}\) n - \(\mu\)

sehingga, X = \(\sum_{i=1}^{24}\) Ui - 7 dimana X ~ Normal (5,2) dan Ui ~ Uniform (0,1)

#1. Membangkitkan Normal (0,1)
i  = 1000
mu = 0 
sigma2 = 1
n  = 12*sigma2
k  = (n/2) - mu 

set.seed(1234)
U  = runif(n*i) #membangkitkan U ~ uniform (0,1)
Um = matrix (U,i)
X  = apply(Um,1,sum)-k

hasil = data.frame( parameter = c("miu","sigma2"), 
                    Nilai= c(mean(X),var(X)))
hasil

  parameter      Nilai
1       miu 0.02175383
2    sigma2 0.97653380

hist(X,col = "Light Coral", main = paste ("Histogram of", "X ~ Normal (0,1)"))

#2. Membangkitkan Normal (5,2)
i  = 1000
mu = 5 
sigma2 = 2
n  = 12*sigma2
k  = (n/2) - mu 

set.seed(1234)
U  = runif(n*i) #membangkitkan U ~ uniform (0,1)
Um = matrix (U,i)
X  = apply(Um,1,sum)-k

hasil = data.frame( parameter = c("miu","sigma2"), 
                    Nilai= c(mean(X),var(X)))
hasil

  parameter    Nilai
1       miu 5.089300
2    sigma2 2.092465

hist(X,col = "Light Coral", main = paste ("Histogram of", "X ~ Normal (5,2)"))

2. Metode Box-Muller

Algoritma Box-Muller digunakan untuk pembangkitan bilangan acak yang berdistribusi gaussian dengan memberikan variasi pada pemilihan skenario yang dibuat untuk mensimulasikan skenario secara dinamik (Thomas et.al, 2007 dalam Ihsan, 2019)

Untuk dua peubah acak menyebar U(0,1) yang saling bebas, U1 dan U2 , maka :

N1 = (-2 \(\log_{e}\) U1 )\(^{1/2}\) cos (2\(\pi\)U2)

N2 = (-2 \(\log_{e}\) U1 )\(^{1/2}\) sin (2\(\pi\)U2)

dimana N1 dan N2 masing-masing menyebar Normal (0,1)

Fungsi di atas dapat dituliskan menjadi :

N1 = R cos \(\Theta\)

N2 = R sin \(\Theta\)

dengan

R = (-2 \(\log_{e}\) U1)\(^{1/2}\) ~ Eks (1/2)

\(\Theta\) = 2\(\pi\)U2 ~ U (0, 2\(\pi\))

(N1 , N2) \(\leftrightarrow\) (R , \(\pi\))

Koordinat Cartesius \(\leftrightarrow\) Koordinat Polar

i  = 1000
set.seed(123)
U1 = runif (i)
U2 = runif (i)
R  = sqrt (-2*log(U1))
Theta = 2*pi*U2
N1 = R*cos(Theta)
N2 = R*sin(Theta)

hasil = data.frame( parameter = c("miu N1","sigma2 N1","miu N2","sigma2 N2"), 
                    Nilai= c(mean(N1),var(N1),mean(N2),var(N2)))
hasil

  parameter       Nilai
1    miu N1 -0.05325227
2 sigma2 N1  1.03414570
3    miu N2 -0.02125216
4 sigma2 N2  0.97728585

par(mfrow = c(1,2))
hist(N1,col = "Misty Rose", main = paste ("Histogram of", "N1 ~ Normal (0,1)"))
hist(N2,col = "Misty Rose", main = paste ("Histogram of", "N2 ~ Normal (0,1)"))

3. Metode Polar-Marsaglia

Metode pembangkit bilangan acak pseudo untuk menghasilkan sepasang independen peubah acak normal baku.Perbedaan dengan metode Box-Muller yaitu terletak pada penggunaan komputasi sinus dan cosinus. Pada Metode Polar Marsaglia tidak digunakan komputasi sinus dan cosinus.

Dari persamaan pada metode Box-Muller :

N1 = (-2 \(\log_{e}\) U1 )\(^{1/2}\) cos (2\(\pi\)U2)

N2 = (-2 \(\log_{e}\) U1 )\(^{1/2}\) sin (2\(\pi\)U2)

menjadi,

N1 = (-2 \(\log_{e}\) R\(^2\) )\(^{1/2}\) V1 (V1\(^2\) + V2\(^2\))\(^{-1/2}\)

N2 = (-2 \(\log_{e}\) R\(^2\) )\(^{1/2}\) V2 (V1\(^2\) + V2\(^2\))\(^{-1/2}\)

menjadi,

N1 = (-2 \(\log_{e}\) (V1\(^2\) + V2\(^2\)) )\(^{1/2}\) V1 (V1\(^2\) + V2\(^2\))\(^{-1/2}\)

N2 = (-2 \(\log_{e}\) (V1\(^2\) + V2\(^2\)) )\(^{1/2}\) V2 (V1\(^2\) + V2\(^2\))\(^{-1/2}\)

sehingga,

N1 = V1 {(-2 \(\log_{e}\) W)/W }\(^{1/2}\)

N2 = V2 {(-2 \(\log_{e}\) W)/W }\(^{1/2}\)

dimana W = V1\(^2\) + V2\(^2\)

i  = 1000
set.seed(123)
U1 = runif (i)
U2 = runif (i)
V1 = 2*U1 - 1
V2 = 2*U2 - 1
W  = V1^2 + V2^2

W1 = W [W<1]
V1 = V1 [W<1]
V2 = V2 [W<1]

i_1 = i-length(W1)
while (i_1>0) {
  U_1 = runif(i_1)
  U_2 = runif(i_1)
  V_1 = 2*U_1-1
  V_2 = 2*U_2-1
  W_1 = V_1^2 + V_2^2
  W1  = c(W1,W_1[W_1<1])
  V1  = c(V1,V_1[W_1<1])
  V2  = c(V2,V_2[W_1<1])
  i_1 = i-length(W1)
}

N1 = V1*sqrt(-2*log(W1)/W1)
N2 = V2*sqrt(-2*log(W1)/W1)

hasil = data.frame( parameter = c("miu N1","sigma2 N1","miu N2","sigma2 N2"), 
                    Nilai= c(mean(N1),var(N1),mean(N2),var(N2)))
hasil

  parameter        Nilai
1    miu N1  0.008316551
2 sigma2 N1  1.062962164
3    miu N2 -0.012600728
4 sigma2 N2  0.979850242

par(mfrow = c(1,2))
hist(N1,col = "Thistle", main = paste ("Histogram of", "N1 ~ Normal (0,1)"))
hist(N2,col = "Thistle", main = paste ("Histogram of", "N2 ~ Normal (0,1)"))

Ketakbiasan (Penduga Tak Bias)

\(\overline{x}\) adalah penduga tak bias bagi \(\mu\) , jika E ( \(\overline{x}\) ) = \(\mu\)

\(s^2\) adalah penduha tak bisa bagi \(\sigma ^{2}\) , jika E ( \(s^2\) ) = \(\sigma ^{2}\)

Maka dalam hal ini kita akan membuktikan apakah benar nilai harapan penduga parameter sama dengan nilai parameternya

Algoritme :

• Tentukan sebaran yang akan digunakan

• Ulangi sebanyak k kali

  • Bangkitkan n buah data dari sebaran yang sudah ditentukan

  • Hitung nilai \(\overline{x}\) dan \(s^2\)

• Hitung rata-rata dari \(\overline{x}\) dan \(s^2\) , kemudian bandingkan dengan \(\mu\) dan \(\sigma ^{2}\)

n = 10
k = 1000 #ulangan
populasi = rnorm(100,10,sqrt(5)) # populasi menyebar normal
miu      = mean(populasi)
sigma2   = var(populasi)*(100-1)/100 #fungsi var di R adalah varian contoh sehingga perlu disesuaikan 
sampel   = matrix(NA,k,n)
for (i in 1:k) sampel[i,] = sample (populasi,n)
xbar     = apply (sampel,1,mean)
s2       = apply(sampel,1,var)
E.xbar   = mean(xbar)
E.s2     = mean(s2)

hasil = data.frame( "."  = c("populasi","E(sampel)"), 
                    mean   = c(miu, E.xbar), 
                    varian = c(sigma2, E.s2))
hasil

          .     mean   varian
1  populasi 9.882038 4.741754
2 E(sampel) 9.901551 4.856654

n = 50 #ukuran contoh
k = 1000 #ulangan
beta       = c(8,20) # beta_0 = 8 dan beta_1 = 20
sigmaerror = 3

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

rataanbetaduga = apply(betaduga,2,mean)
names(rataanbetaduga) = c("bo","b1")
rataanbetaduga

       bo        b1 
 7.970786 20.005076 

Selang Kepercayaan

Apa arti dari SK 95% ?

• SK 95% bagi \(\theta\) : Kita percaya 95% bahwa selang a sampai b memuat nilai parameter \(\theta\) yang sebenarnya.

• SK 95% : Jika kita melakukan 100 kali percontohan acak dan setiap percontohan acak dibuat selang kepercayaannya, maka dari 100 SK yang terbentuk, ada 95 SK yang mencakup parameter sedangkan sisanya sebanyak 5 SK tidak mencakup parameter.

Algoritme :

• Tentukan sebaran yang akan digunakan

• Ulangi sebanyak k kali

  • Bangkitkan n buah data dari sebaran yang sudah ditentukan

  • Hitung nilai \(\overline{x}\) dan \(s^2\)

  • Hitung \(\sigma^{2}_\overline{x}\) dan buat selang kepercayaann (1- \(\alpha\) )%

• Hitung proporsi banyaknya selang kepercayaan yang memuat \(\mu\), bandingkan dengan (1-\(\alpha\))

n    = 50
k    = 100 #ulangan
alpha  = 0.05
mu   = 50
std  = 10

set.seed(503)
sampel = matrix(rnorm(n*k,mu,std),k)
xbar = apply(sampel,1,mean)
s    = apply(sampel,1,sd)
SE   = s/sqrt(n)
z    = qnorm(1-alpha/2)
SK   = (xbar-z*SE < mu & mu < xbar+z*SE)

sum(SK)/k #proporsi banyaknya SK yang memuat mu

[1] 0.95

matplot(rbind (xbar-z*SE, xbar+z*SE), rbind(1:k,1:k), col=ifelse(SK,"blue","red"), type = "l", lty = 1,main='Selang Kepercayaan', xlab='SK', ylab='banyak ulangan')
abline(v=mu)

1. Bootstrap

Ukuran n kecil mengakibatkan beberapa teori peluang dan statistika inferensia tidak terpenuhi. Selain itu, ukuran n kecil juga menjadi masalah dalam suatu pendugaan parameter yang berkaitan dengan tingkat keakuratan dan ketelitian.

• Ketika n kecil dan tidak ada informasi sebaran populasi dapat diatasi dengan Metode Resampling Bootstap

• Ide dasar dari bootstrap adalah membangun data semu dengan menggunakan informasi data asli

• Efisien untuk mendapatkan hasil yang robust

Ilustrasi :

Dugalah rasio item1 dengan total item menggunakan bootstrap !

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

[1] 0.5563869

se_b

[1] 0.03749447

2. Monte Carlo

• Simulasi yang memanfaatkan informasi mengenai sebaran data yang diketahui (dihipotesiskan, dianggap tahu) dengan pasti

• Simulasi didasarkan pada pembangkitan bilangan acak dari sebaran hipotetik

• Perilaku yang ingin dipelajari dapat diketahui jika proses diulang berkali-kali

• Banyak digunakan untuk "mengetahui" sebaran dari statistik

Ilustrasi :

1. Misal X ~ Normal (0,1), jika diketahui Y = 2X, berapakah P(Y>2) ?

Algoritme :

1. Bangkitkan X ~ N(0,1)

2. Hitung Y = 2X

3. Ulangi langkah 1 dan 2 sebanyak 100000 kali

4. Hitung presentasi nilai Y > 2

n  = 100000
x  = rnorm(n)
y  = 2*x
y1 = ifelse (y > 2, 1, 0)
mean(y1) #P(Y>2) secara empirik

[1] 0.16004

pnorm(2,0,2,lower.tail=F) #P(Y>2 dari sebaran hipotetik)

[1] 0.1586553

lambda = 1
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

[1] 0.7805

kesimpulan

[1] "Tak Tolak H0"

3. Jacknife

• Misalkan terdapat sampel \(y = ( y_1, y_2, ... , y_n )\) untuk menduga \(\theta\) dengan penduga \(\widehat{\theta} = f (y)\)

• Tujuannya adalah menduga SE (\(\widehat{\theta}\))

• Sampel observasi dari \(y_i = (y_i,...,y_{i-1} , y_{i+1},.. ,y_n )\) untuk i = 1 ,.., n disebut Jacknife Samples

Ilustrasi :

Dugalah rasio item1 dengan total item menggunakan Jacknife !

Algoritme :

1. Hitung \(\widehat{\theta}\) = r = \(\frac{\overline{y}}{\overline{x} + \overline{y}}\)

2. Tentukan sampel janknife \(x_(i)\) dan \(y_(i)\) ; i = 1, ... ,n

3. Hitung penduga jacknife :

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

[1] 0.5584068

se_jk

[1] 0.006207102