Terdapat 2 problem yang harus diselesaikan oleh mahasiswa. Pada problem #1, mahasiswa diminta untuk membangkitkan data dari dua gugus yang berbeda, yakni data dengan sebaran simetris dan sebaran kombinasi chisq + normal. Lalu, coba identifikasi apakah data simetris cenderung lebih cepat mendekati sebaran normal? Resampling ini tentunya dilakukan dengan menggunakan jumlah sampel yang berbeda.

Pada problem #2, diberikan soal sedikit lebih kompleks. Berikut studi kasus yang harus dipecahkan melalui soal ini

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 1000 dan ukuran sampel yang diteliti adalah 10,20,30,40,50,60,70,80,90,100. Metode mana yang lebih baik?

Untuk menjawab problem #1, disajikan gambaran besar grafik sebaran yang akan digunakan

# sumber gugus data bangkitan, populasi sebesar 1jt
set.seed(121)
pop.norm = rnorm(1000000) #Norm(0,1)
pop.chisq = rchisq(1000000,5) #Chisq(db= 5)

par(mfrow= c(1,2))
# kurva normal pada pop.norm
hist(pop.norm, prob=TRUE, main="Normal Distribution", xlab="Value", col="lightblue")
curve(dnorm(x, mean=mean(pop.norm), sd=sd(pop.norm)), add=TRUE, col="red", lwd=2)

# kurva chisq untuk pop.chisq
hist(pop.chisq, prob=TRUE, main="Chi-Square Dist (df=5)", xlab="Value", col="lightgreen")
curve(dchisq(x, df=5), add=TRUE, col="blue", lwd=2)

Problem 1

a. Bangkitan data sebaran normal dengan μ = 0 dan σ =1.

set.seed(121)

par(mfrow=c(2,2))
generate_normal_hist <- function(n, mean = 0, sd = 1) {
  pop <- rnorm(1000000, mean, sd)
  samp <- sample(pop, n)
  
  hist(samp, freq = FALSE, col = "lightblue", main = paste0("Sample size n = ", n),
       xlab = "Value", ylab = "Density")
  
  curve(dnorm(x, mean, sd), add = TRUE, col = "red", lwd = 2)
}

# Contoh penggunaan fungsi untuk menghasilkan plot histogram dengan garis kurva
generate_normal_hist(4); generate_normal_hist(10); generate_normal_hist(12); generate_normal_hist(100)

Pada gugus data campuran, gunakan fungsi sebagai berikut;

generate_data_normal_chisq <- function(n, mu=0, sigma=1, df=1) {
  n_norm <- floor(n/2)
  n_chisq <- n - n_norm
  norm_data <- rnorm(n_norm, mean=mu, sd=sigma)
  chisq_data <- rchisq(n_chisq, df=df)
  data <- c(norm_data, chisq_data)
  return(data)
}

generate_data_chisq_ab <- function(n, a=1, b=5) {
  n_half <- floor(n/2)
  chisq_data_a <- rchisq(n_half, df=a)
  chisq_data_b <- rchisq(n-n_half, df=b)
  data <- c(chisq_data_a, chisq_data_b)
  return(data)
}

generate_data_mixed <- function(n, a=1, b=5, mu=0, sigma=1) {
  n_chisq_a <- floor(n/4)
  n_chisq_b <- floor(n/4)
  n_norm_a <- floor(n/4)
  n_norm_b <- n - n_chisq_a - n_chisq_b - n_norm_a
  chisq_data_a <- rchisq(n_chisq_a, df=a)
  chisq_data_b <- rchisq(n_chisq_b, df=b)
  norm_data_a <- rnorm(n_norm_a, mean=mu, sd=sigma)
  norm_data_b <- rnorm(n_norm_b, mean=mu, sd=sigma)
  data <- c(chisq_data_a, chisq_data_b, norm_data_a, norm_data_b)
  return(data)
}

b. Bangkitkan gugus data campuran dengan menggunakan fungsi di atas

Adapun komposisi dari ketiga data yang ingin dibangkitkan adalah sebagai berikut:

Data Normal - Chisq: 50% menyebar normal + 50% menyebar chisquare
Data Chisq (parameter a dan b): 50% sebaran chisquare dengan parameter a + 50% sebaran chisquare dengan parameter b
Data Mixed: 25% sebaran chisquare dengan parameter a + 25% sebaran chisquare dengan parameter b + 25% sebaran normal dengan parameter a + 25% sebaran normal dengan parameter b

set.seed(121)
n_samples <- c(4, 12, 60, 100)

# Data Normal - Chisq
data_norm_chisq <- lapply(n_samples, function(n) generate_data_normal_chisq(n))

# Data Chisq (parameter a dan b)
data_chisq_ab <- lapply(n_samples, function(n) generate_data_chisq_ab(n))

# Data Mized
data_mixed <- lapply(n_samples, function(n) generate_data_mixed(n))

par(mfrow=c(2,2))
for (i in seq_along(n_samples)) {
  hist(data_norm_chisq[[i]], main = paste0("Normal and Chi-square (n = ", n_samples[i], ")"))
  lines(density(data_norm_chisq[[i]]), col = "blue")
}

for (i in seq_along(n_samples)) {
  hist(data_chisq_ab[[i]], main = paste0("Chi-square with a and b (n = ", n_samples[i], ")"))
  lines(density(data_chisq_ab[[i]]), col = "blue")
}

for (i in seq_along(n_samples)) {
  hist(data_mixed[[i]], main = paste0("Mixed distribution (n = ", n_samples[i], ")"))
  lines(density(data_mixed[[i]]), col = "blue")
}

Hasil observasi di atas menunjukkan bahwa semakin besar jumlah sampel, kecenderungan sebaran gugus data akan mendekati sebaran normal. Jika dilihat dari garis kurva, pada bangkitan sebaran normal, dengan menggunakan n yang tidak besar, kurva sudah mendekati sebaran normal. Ditambah jika n yang digunakan semakin besar, tentu kurva akan menunjukkan pola yang mendekati kesimetrisan

Problem 2

Bangkitan data dengan distribusi multivariat normal dan data outliers sebaran tertentu

# Fungsi untuk mengganti data dengan outlier secara acak
generate_data <- function(n) {
  # Set parameters
  mean <- c(10, 12) #Vektor rataan(10,12)
  sd <- c(2, 3) #Vektor simp baku(2,3)
  correlation <- 0.85 #Korelasi
  outlier_prop <- 0.2 #Outlier sebanyak 20%
  outlier_mean <- c(25, 25) #Mean outlier = 25
  outlier_sd <- c(2, 2) #Std outlier = 2
  
  # Generate normal data
  data <- MASS::mvrnorm(n, mu = mean, Sigma = matrix(c(sd[1]^2, sd[1]*sd[2]*correlation, sd[1]*sd[2]*correlation, sd[2]^2), nrow = 2))

  # Generate outliers
  num_outliers <- round(n*outlier_prop)
  outlier_index <- sample(1:n, num_outliers)
  data[outlier_index,] <- MASS::mvrnorm(num_outliers, mu = outlier_mean, Sigma = matrix(c(outlier_sd[1]^2, 0, 0, outlier_sd[2]^2), nrow = 2))
  
  # Randomly replace values with outliers
  replace_index <- sample(1:n, num_outliers)
  data[replace_index,] <- MASS::mvrnorm(num_outliers, mu = outlier_mean, Sigma = matrix(c(outlier_sd[1]^2, 0, 0, outlier_sd[2]^2), nrow = 2))
  
  return(data)
}

Dari bangkitan di atas, nilai korelasi yang ditetapkan yakni sebesar 0.85. Perlu diperhatikan bahwa outlier sebesar 20% bisa jadi digolongkan bukan sebagai outlier. Akan tetapi, outlier sebesar 20% hanya berupa asumsi saja, jangan dihiraukan.

Fungsi untuk membangkitkan data acak dan data outlier

# Fungsi untuk mengganti data dengan outlier secara acak
replace_outliers_randomly <- function(data) {
  # Hitung nilai batas untuk masing-masing variabel
  q1 <- apply(data, 2, quantile, probs=0.25, na.rm=TRUE)
  q3 <- apply(data, 2, quantile, probs=0.75, na.rm=TRUE)
  iqr <- q3 - q1
  upper_bound <- q3 + 1.5 * iqr
  lower_bound <- q1 - 1.5 * iqr
  
  # Ganti nilai outlier dengan data acak yang memiliki rentang nilai yang sama
  for (i in seq_along(upper_bound)) {
    outlier_index <- which(data[,i] > upper_bound[i] | data[,i] < lower_bound[i])
    n_outliers <- length(outlier_index)
    if (n_outliers > 0) {
      random_data <- runif(n_outliers, min=lower_bound[i], max=upper_bound[i])
      data[outlier_index,i] <- random_data
    }
  }
  return(data)
}

# Fungsi untuk menghitung korelasi
calculate_correlations <- function(data) {
  # Hitung korelasi Pearson dan Spearman
  cor_pearson <- cor(data, method="pearson", use="pairwise.complete.obs")
  cor_spearman <- cor(data, method="spearman", use="pairwise.complete.obs")
  
  # Masukkan hasil perhitungan ke dalam sebuah list
  result <- list(pearson=cor_pearson[1,2], spearman=cor_spearman[1,2])
  return(result)
}

Gabungan data acak dengan data outlier

set.seed(123)
n_samples <- c(10, 20, 30, 40, 50, 60, 70, 80, 90, 100)
n_trials <- 1000
correlations <- matrix(0, nrow=n_trials, ncol=2)

results <- data.frame()

for (i in seq_along(n_samples)) {
  for (j in 1:n_trials) {
    data <- generate_data(n_samples[i])
    data <- replace_outliers_randomly(data)
    cor <- calculate_correlations(data)
    correlations[j, 1] <- cor$pearson
    correlations[j, 2] <- cor$spearman
  }
  mean_correlations <- apply(correlations, 2, mean)
  
  # Tambahkan hasil perhitungan korelasi ke dataframe
  results <- rbind(results, data.frame(sample_size=n_samples[i], 
                                       mean_correlation_pearson=mean_correlations[1], 
                                       mean_correlation_spearman=mean_correlations[2]))
}

results

##    sample_size mean_correlation_pearson mean_correlation_spearman
## 1           10                0.9471366                 0.8783879
## 2           20                0.9511094                 0.8959835
## 3           30                0.9514735                 0.9015751
## 4           40                0.9512274                 0.9042769
## 5           50                0.9514346                 0.9059175
## 6           60                0.9508133                 0.9068146
## 7           70                0.9510756                 0.9082908
## 8           80                0.9506812                 0.9076597
## 9           90                0.9502162                 0.9084890
## 10         100                0.9502068                 0.9075938

Dari dataframe di atas, dapat disimpulkan bahwa korelasi pearson menunjukkan nilai yang lebih besar, yakni >0.94. Begitupun dengan progresi atau perubahan nilai korelasi jika menggunakan jumlah sampel yang berbeda. Akan tetapi, kesimpulan yang digunakan berdasarkan pada selisih (absolut) terhadap nilai yang sudah ditetapkan sebelumnya, yakni 0.85. Oleh karena itu, dapat dikatakan bahwa korelasi spearman adalah metode yang lebih baik dibandingkan dengan korelasi pearson

Pembangkitan Bilangan Acak Multivariat Normal dan Acak Campuran (Non-simetris)

2023-05-26