Monte Carlo

Alfa Nugraha¹

2023-12-05

Install library berikut jika belum ada

install.packages("ggplot2")
install.packages("mvtnorm")
install.packages("MBESS")

Latihan

Latihan 1

Seorang peneliti ingin melihat bagaimana kesalahan tipe 1 (\(\alpha\)) dari uji t (t test) dan uji peringkat bertanda Wilcoxon (Wilcoxoan signed-rank test) [https://sixsigmastudyguide.com/1-sample-wilcoxon-non-parametric-hypothesis-test/] untuk satu populasi pada berbagai ukuran sampel. Hipotesis yang akan diujikan adalah sebagai berikut:

\[ H_{0}: \mu=1500 \\ H_{1}: \mu \neq 1500 \] Bangkitan data berdasarkan distribusi \(N(\mu=1500,\sigma=200)\), \(\alpha=0.05\) dan ulangan sebanyak 10.000. Ukuran sampel yang diteliti adalah 5 sampai 50. Berilah kesimpulan atas hasilnya!

Jawaban 1

Langkah 1 mendefinisikan parameter-parameter yang digunakan untuk simulasi

# mendefinisikan ukuran sampel dari 5 saampai 50
n <- seq(5,50) 
# mendefinisikan miu distribusi normal
mu0 <- 1500
# mendefinisikan sigma distribusi normal
sigma0 <- 200
# mendefinisikan banyaknya ulangan
num_rep <- 10000
# mendefinisikan alpha
alpha <- 0.05

Langkah 2 mendefinisikan fungsi untuk mengextract p-value dari masing-masing uji (test).

t_test_pvalue <- function(x,alt="two.sided",mu=mu0){
  t.test(x, alternative = alt, mu = mu)$p.value  
}

wilcox_test_pvalue <- function(x,alt="two.sided",mu=mu0){
  wilcox.test(x, alternative = alt, mu = mu)$p.value  
}

Langkah 3 menjalankan simulasi

# membuat objek data.frame dummy untuk menyimpan hasil simulasi
result_sim1 <- data.frame(n=0,alpha_empiris=0,test_name="")
set.seed(2020)
# iterasi simulasi berdasarkan ukuran sampel
for(i in seq_along(n)){
  # membangkitan data berdistribusi normal dengan ulangan 10rb
  sim_mc_norm <- replicate(num_rep,                     rnorm(n[i],mu0,sigma0))
  # menghitung p value uji t dari 10rb data hasil bangkitan
  t_test_res <- apply(sim_mc_norm,2,t_test_pvalue)
  # menghitung p value uji wilcoxon dari 10rb data hasil bangkitan
  wilcox_test_res <- apply(sim_mc_norm,2,wilcox_test_pvalue)
  # menghitung alpha empiris dari uji t
  t_p_value_sim <- mean(t_test_res<alpha)
  # menghitung alpha empiris dari uji wilcoxon
  wilcox_p_value_sim <- mean(wilcox_test_res<alpha)
  # menyimpan hasil simulasi 
  result_sim1 <- rbind(
    result_sim1,list(n[i],t_p_value_sim,"t_test"),
    list(n[i],wilcox_p_value_sim,"wilcox_test")
  )
}
# membuang baris dummy
result_sim1 <- result_sim1[-1,]

Langkah 4 membuat grafik dari hasil simulasi

ggplot(data = result_sim1,aes(n,alpha_empiris)) +
  geom_point(aes(color=test_name)) + 
  geom_hline(yintercept = 0.05) + 
  scale_y_continuous(limits = c(0,0.1))

Latihan 2

Seorang peneliti ingin melihat bagaimana nilai koefisien korelasi dari metode korelasi pearson dan metode korelasi spearman untuk satu populasi pada berbagai ukuran sampel jika datanya terdapat pencilan. Bangkitan data pencilan dengan cara berikut:

Bangkitkan data berdasarkan distribusi multivariate Normal dengan vektor rataan \((10,12)\) dan korelasi 0.85 serta vektor simpangan baku \((2,3)\)
Bangkitan outlier pada kedua peubah dengan berdasarkan distribusi \(N(25,2)\) sebanyak 20% dari ukuran sampel pada langkah a.
Ganti amatan dari hasil bangkitan pada langkah a dengan hasil bangkitan pada langkah b secara acak.

Bangkitkan data dengan ulangan 10000 dan ukuran sampel yang diteliti adalah 10,20,30,40,50,60,70,80,90,100. Berikan kesimpulan dari hasil simulasi!

Jawaban 2

Langkah 1 abstraksi bangkitan data

langkah a

set.seed(1212)
# mendefinisikan vektor simpangan baku
sd0 <- c(2,3)
# mendefinisikan matrix korelasi
cor_mat <- matrix(c(1,0.85,0.85,1),nrow=2)
# mengubah matrix korelasi ke kovarians
cov_mat <- MBESS::cor2cov(cor_mat,sd0)
# mendefinisikan vektor rataan
mu_mvn <- c(10,12)
# membangkitkan data berdistribus multivariate normal
x_mvn <- rmvnorm(n=40, mean=c(10,12), sigma=cov_mat)

langkah b

set.seed(1212)
#membangkitkan data berdasarkan distribusi normal
x1_out <- rnorm(8,25,2)
y1_out <- rnorm(8,25,2)
x1y1_out <-cbind(x1_out,y1_out)

langkah c

# memilih secara acak amatan yang akan diganti dengan pencilan
index_out <- sample.int(nrow(x_mvn),8)
# mengganti amatan dari langkah a dengan amatan dari langkah b
x_mvn[index_out,] <- x1y1_out

qplot(x = V1,y=V2,data = as.data.frame(x_mvn))

## Warning: `qplot()` was deprecated in ggplot2 3.4.0.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Langkah 2 mendefinisikan parameter-parameter yang digunakan untuk simulasi

# mendefinisikan ukuran sampel yang dikaji
n <- seq(10,100,10)
# mendefinisikan jumlah ulangan
num_rep2 <- 10000
# mendefinisikan vektor simpangan baku
sd0 <- c(2,3)
# mendefinisikan matrix korelasi
cor_mat <- matrix(c(1,0.85,0.85,1),nrow=2)
# mengubah matrix korelasi ke kovarians
cov_mat <- MBESS::cor2cov(cor_mat,sd0)
# mendefinisikan vektor rataan
mu_mvn <- c(10,12)
# mendefinisikan mu dan sigma untuk data pencilan
mu1 <- 25
sigma1 <- 2

Langkah 3 mendefinisikan fungsi-fungsi yang dibutuhkan

# generate multivariate normal dengan pencilan
gen_mvn_out <- function(n0,mu0,sigma0,n1,mu1,sigma1){
x_mvn <- rmvnorm(n=n0, mean=mu0, sigma=sigma0)
x1_out <- rnorm(n1,mu1,sigma1)
y1_out <- rnorm(n1,mu1,sigma1)
x1y1_out <-cbind(x1_out,y1_out)
index_out <- sample.int(nrow(x_mvn),n1)
x_mvn[index_out,] <- x1y1_out
return(x_mvn)
}

# extract koefisien korelasi
cor_res <- function(x_mat,method){
  cor(x_mat[,1],x_mat[,2],method = method)
}

Langkah 4 menjalankan simulasi

set.seed(1212)
# membuat objek data.frame dummy untuk menyimpan hasil simulasi
result_sim2 <- data.frame(n=0,correlation=0,type="")
for(i in seq_along(n)){
  # membangkitan data berdistribusi multivariate normal dengan ulangan 10rb
  sim_mc_mvn <- replicate(num_rep2, 
                          gen_mvn_out(n0 = n[i], mu_mvn,cov_mat,n1=0.2*n[i],mu1,sigma1),
                          simplify = F)
  # menghitung keofisien korelasi pearson dari 10rb data hasil bangkitan
  pearson_res <- sapply(sim_mc_mvn,cor_res,method="pearson")
  # menghitung keofisien spearman pearson dari 10rb data hasil bangkitan
  spearman_res <- sapply(sim_mc_mvn,cor_res,method="spearman")
  # menghitung rata-rata koefisien
  pearson_sim <- mean(pearson_res)
  spearman_sim <- mean(spearman_res)
  
  # menyimpan hasil simulasi
  result_sim2 <- rbind(result_sim2, 
                       list(n[i],pearson_sim,"pearson"),
                       list(n[i],spearman_sim,"spearman")
  )
}
result_sim2 <- result_sim2[-1,]

Langkah 5 membuat grafik berdsarkan hasil simulasi

ggplot(data = result_sim2,aes(n,correlation)) +
  geom_point(aes(color=type)) +
  geom_hline(yintercept = 0.85) +
  scale_y_continuous(limits = c(0.5,1))

Latihan 3

Bangkitkan peubah \(Y, X_{1},X_{2}\) dan \(X_{3}\) berdasarkan model regresi linear berganda berikut ini:

\(Y= 10 + 3X_{1} + 5X_{2} + 7X_{3} + \epsilon\) dengan mengasumsikan bahwa \(\epsilon \sim N(0,1)\). Banyaknya amatan yang dibangkitkan adalah 1000

Jawaban 3

Bangkitkan residual berdasarkan distribusi normal \(\epsilon \sim N(0,\sigma)\)

#banyaknya amatan
n <- 1000
set.seed(12)
epsilon <- rnorm(n,mean = 0,sd = 1)

Bangkitkan peubah penjelas \(X_{i}\) yang saling bebas.

set.seed(10)
X <- replicate(3,rnorm(n,mean=0,sd=1.5))
head(X)

##             [,1]       [,2]        [,3]
## [1,]  0.02811926  1.5750205 -0.46179747
## [2,] -0.27637881  0.4291389  1.13712837
## [3,] -2.05699582  0.3608472 -0.86079512
## [4,] -0.89875157  1.2490578 -1.40811669
## [5,]  0.44181769 -0.3344748 -0.04154899
## [6,]  0.58469145  0.4325163 -1.59937298

Pembangkitan peubah penjelas \(X_{i}\) dapat dilakukan juga dengan fungsi mvnorm jika ingin melakukan simulasi multikolinearitas.

Bangkitkan peubah respon \(Y\) dengan menggunakan model regresi

betas <- matrix(c(10,3,5,7),nrow = 4,ncol = 1)
Xgab <- cbind(1,X)

Y <- Xgab%*%betas+epsilon

dataRegresi <- data.frame(Y,X)
colnames(dataRegresi) <- c("Y","X1","X2","X3")
head(dataRegresi)

Note : sintaks Xgab <- cbind(1,X) digunakan untuk menggabungkan peubah penjelas yang ada di matrik \(X\) dengan vektor yang isinya 1. Vektor ini digunakan sebagai pengali konstanta model (dalam hal ini bernilai 10).

Untuk membuktikan bahwa data bangkitan sudah tepat, maka data tersebut dimodelkan dengan fungsi lm yang ada di R.

modelRegresi <- lm(Y~.,data = dataRegresi)
summary(modelRegresi)

## 
## Call:
## lm(formula = Y ~ ., data = dataRegresi)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.0364 -0.6052 -0.0207  0.6087  3.1774 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  9.97484    0.03033   328.8   <2e-16 ***
## X1           2.97225    0.02042   145.6   <2e-16 ***
## X2           4.97282    0.01939   256.5   <2e-16 ***
## X3           7.00493    0.01989   352.2   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9589 on 996 degrees of freedom
## Multiple R-squared:  0.9955, Adjusted R-squared:  0.9954 
## F-statistic: 7.275e+04 on 3 and 996 DF,  p-value: < 2.2e-16

Berdasarkan output regresi diatas, nilai koefisien regresi yang dihasilkan mendekati nilai koefisien regresi yang telah dibangkitkan. Hal ini membuktikan bahwa data sudah sesuai dengan model regresi yang diinginkan.

Latihan 4

Berdasarkan Latihan 3 lakukan pengulangan sebanyak 100 kali kemudian hitung rata-rata dan simpangan baku dari masing-masing koefisien regresi.

Jawaban 4

#banyaknya amatan
n <- 1000
ulangan <- 100

# membuat tempat penyimpanan
# koefisien tiap ulangan
koefisien_ulangan <- NULL
X <- replicate(3,rnorm(n,mean=0,sd=1.5))
set.seed(13)
for(i in seq(ulangan)){
  epsilon <- rnorm(n,mean = 0,sd = 1)
  Xgab <- cbind(1,X)
  Y <- Xgab%*%betas+epsilon
  dataRegresi <- data.frame(Y,X)
  colnames(dataRegresi) <- c("Y","X1","X2","X3")
  modelRegresi <- lm(Y~.,data = dataRegresi)
  # rbind digunakan untuk menggabungkan baris
  koefisien_ulangan <- rbind(koefisien_ulangan,coef(modelRegresi))
}
head(koefisien_ulangan)

##      (Intercept)       X1       X2       X3
## [1,]    9.996860 3.013150 5.007435 7.008822
## [2,]   10.002799 3.011432 5.057188 6.976413
## [3,]   10.016592 2.986175 5.033262 6.959016
## [4,]    9.974942 3.025846 5.014508 7.017001
## [5,]    9.999086 2.987499 4.950832 6.984672
## [6,]   10.026188 2.991365 4.971202 6.996268

#menghitung mean dari setiap koefisien
colMeans(koefisien_ulangan)

## (Intercept)          X1          X2          X3 
##   10.002224    3.002941    5.002935    7.001979

#menghitung sd dari setiap koefisien
apply(koefisien_ulangan,2,sd)

## (Intercept)          X1          X2          X3 
##  0.03058885  0.01977378  0.02204635  0.02122027

Statistika dan Sains Data, IPB University, alfanugraha@apps.ipb.ac.id ↩︎