Pembangkitan Bilangan Acak Sebaran Normal

Sebaran Normal dapat di bangkitkan melalui beberapa metode berikut.

Dalil Limit Pusat

dengan n dan k didapatkan dari rumus berikut.

Berikut adalah contoh pembangkitan sebaran normal baku.

i <- 1000
mu <- 0
sigma2 <- 1
n <- 12*sigma2
k <- mu-n/2
U <- runif(n*i)
Um <- matrix(U,i)
N <- apply(Um,1,sum)+k
mean(N)
## [1] -0.04647277

Bisa dilihat rata-rata dari data yang telah dibangkitkan mendekati 0.

var(N)
## [1] 0.9920685

Bisa dilihat ragam dari data yang telah dibangkitkan mendekati 1.

hist(N)

Bisa dilihat histogram dari data yang dibangkitkan memiliki sebaran yang simetrik sesuai dengan karakteristik sebaran normal.

Berikut adalah contoh pembangkitan sebaran normal dengan nilai tengah sebesar 4 dan ragam sebesar 16.

i <- 1000
mu <- 4
sigma2 <- 16
n <- 12*sigma2
k <- mu-n/2
U <- runif(n*i)
Um <- matrix(U,i)
M <- apply(Um,1,sum)+k
mean(M)
## [1] 3.776246

Bisa dilihat rata-rata dari data yang dibangkitkan mendekati 4.

var(M)
## [1] 14.83512
hist(M)

Bisa dilihat histogram dari data yang dibangkitkan memiliki sebaran yang simetrik sesuai dengan karakteristik sebaran normal.

Metode Box-Muller

Pendekatan ini mengubah data dengan koordinat kartesius menjadi koordinat polar dengan langkah berikut.

Berikut adalah contoh pembangkitan sebaran normal dengan metode Box-Muller.

i <- 1000
U1 <- runif(i)
U2 <- runif(i)
R <- sqrt(-2*log(U1))
Theta <- 2*pi*U2
N1 <- R*cos(Theta)
N2 <- R*sin(Theta)
mean(N1) 
## [1] 0.02818737
var(N1)
## [1] 0.9671955
mean(N2)
## [1] -0.005411695
var(N2)
## [1] 1.055862

dapat dilihat bahwa nilai rata-rata baik N1 maupun N2 mendekati 0 dan ragamnya mendekati 1.

hist(N1)

hist(N2)

Kedua histogram di atas memiliki sebaran yang simetrik (data menyebar normal).

Metode Polar Marsaglia

Merupakan pengembangan dari metode box muller dengan cara sebagai berikut.

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

M1 <- V1*sqrt(-2*log(W1)/W1)
M2 <- V2*sqrt(-2*log(W1)/W1)
mean(M1) 
## [1] 0.006909552
var(M1)
## [1] 1.078456
mean(M2)
## [1] -0.02891576
var(M2)
## [1] 0.9588654

dapat dilihat bahwa nilai rata-rata baik M1 maupun M2 mendekati 0 dan ragamnya mendekati 1.

hist(M1)

hist(M2)

Kedua histogram di atas memiliki sebaran yang simetrik (data menyebar normal).

Sebaran Percontohan & Dalil Limit Pusat

Sebaran percontohan ialah sebaran peluang dari suatu statistik tertentu. Dalam materi ini akan lebih didalami terkait sebaran percontohan untuk rataan dengan menggunakan dalil limit pusat.

Simulasi dapat dilakukan dengan menggunakan sebaran apapun. dalam contoh ini akan digunakan sebaran normal, eksponensial dan seragam.

Berikut adalah simulasi pada populasi tak terhingga dengan jumlah contoh n ={10,30,100}.

#Normal n=10
k<-1000
n<-10
x11<-matrix(rnorm(n*k),k)
x11<-apply(x11,1,mean)
mean(x11);var(x11)
## [1] 0.006389233
## [1] 0.1000967
#Normal n=30
k<-1000
n<-30
x12<-matrix(rnorm(n*k),k)
x12<-apply(x12,1,mean)
mean(x12);var(x12)
## [1] -0.001389417
## [1] 0.03297246
#Normal n=100
k<-1000
n<-100
x13<-matrix(rnorm(n*k),k)
x13<-apply(x13,1,mean)
mean(x13);var(x13)
## [1] 0.002965743
## [1] 0.009353753
#Eksponensial n=10
k<-1000
n<-10
x21<-matrix(rexp(n*k),k)
x21<-apply(x21,1,mean)
mean(x21);var(x21)
## [1] 0.9950717
## [1] 0.1036051
#Eksponensial n=30
k<-1000
n<-30
x22<-matrix(rexp(n*k),k)
x22<-apply(x22,1,mean)
mean(x22);var(x22)
## [1] 0.9999216
## [1] 0.0316938
#Eksponensial n=100
k<-1000
n<-100
x23<-matrix(rexp(n*k),k)
x23<-apply(x23,1,mean)
mean(x23);var(x23)
## [1] 1.002311
## [1] 0.009492788
#Seragam n=10
k<-1000
n<-10
x31<-matrix(runif(n*k),k)
x31<-apply(x31,1,mean)
mean(x31);var(x31)
## [1] 0.4991981
## [1] 0.008488607
#Seragam n=30
k<-1000
n<-30
x32<-matrix(runif(n*k),k)
x32<-apply(x32,1,mean)
mean(x32);var(x32)
## [1] 0.4999294
## [1] 0.002482794
#Seragam n=100
k<-1000
n<-100
x33<-matrix(runif(n*k),k)
x33<-apply(x33,1,mean)
mean(x33);var(x33)
## [1] 0.4991775
## [1] 0.0008637503

Berikut adalah histrogram hasil simulasi tersebut.

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

terlihat bahwa pada sebaran dengan n kecil, sebaran contoh cenderung mengikuti sebaran populasinya sedangkan sebaran contoh dengan n besar cenderung simetrik.

Ketakbiasan dan Selang Kepercayaan

Ketakbiasan Xbar dan S sebagai penduga miu dan sigma

n<-10
k<-1000 
pop<-rnorm(100,10, sqrt(5)) #pop terhingga
mean(pop)
## [1] 10.15971
var(pop) *(100-1)/100  #var=ragam contoh
## [1] 5.870086
cth<-matrix(NA,k,n)
for(i in 1:k){
  cth[i,]<-sample(pop,n)
}

xbar<-apply(cth,1,mean)
s2<-apply(cth,1,var)
mean(xbar)
## [1] 10.17534
mean(s2)
## [1] 5.953711

Bisa dilihat bahwa penduganya mendekati nilai parameter yang sebenarnya.

Ketakbiasan penduga OLS dalam regresi linear

n<-20
ulangan<-100
beta<-c(8,20)
sigma2error<-3 

betaduga<-matrix(NA,ulangan,2)
for(i in 1:ulangan){
  x<-runif(n)*10
  epsilon<-rnorm(n,0,sqrt(sigma2error))
  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 
##  8.046486 20.005904

bisa dilihat nilai penduganya mendekati nilai sebenarnya, yaitu b0 = 8 dan b1 = 20.

Selang Kepercayaan

n<-50
ulangan<-50
alpha<-0.05
mu<-50
std<-10
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.96
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)

plot di atas merupakan gambaran dari selang kepercayaan yang mengandung nilai parameter sebenarnya sementara yang berwarna merah adalah selang kepercayaan yang tidak megandung nilai parameter yaitu miu = 50.

Resampling

Simulasi Monte Carlo

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

output di atas menggambarkan kemiripan hasil dari peluang secara empirik (atas) dengan melalui sebaran hipotetik (bawah).

lambda<-5
k<-10000
data1<-c(3,8,4,2,5,5,6,7,4,1,3,3,4,4,3)
n<-length(data1)
m<-mean(data1)
s<-sd(data1)
tdata<-abs((m-lambda)/(s/sqrt(n))) 

data2<-matrix(rpois(n*k,lambda),k)
m1<-apply(data2,1,mean)
s1<-apply(data2,1,sd)
tdist<-abs((m1-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.0932
## [1] "Tak Tolak H0"

Simulasi 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;se_jk
## [1] 0.5584068
## [1] 0.006207102

Output di atas merupakan hasil pendugaan terhadap teta dan standard error menggunakan metode jacknife.

Simulasi 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
x<-matrix(sample(data1$item1,n*b,replace=T),b)
y<-matrix(sample(data1$item2,n*b,replace=T),b)
xbar<-apply(x,1,mean)
ybar<-apply(y,1,mean)
teta_i<-xbar/(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; se_b
## [1] 0.5578005
## [1] 0.03676413

Output di atas merupakan hasil pendugaan terhadap teta dan standard error menggunakan metode bootstrap.