MSR Kelp Kuliah No 3

Syntaks Soal No 3

# Packages
library(tidyverse)

## Warning: package 'tidyverse' was built under R version 4.3.2

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.3     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.4     ✔ tibble    3.2.1
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(ggpubr)
library(knitr)
library(kableExtra)

## Warning: package 'kableExtra' was built under R version 4.3.3

## 
## Attaching package: 'kableExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     group_rows

# Mengatur parameter awal
n_sim <- 1000
n_obs <- c(5, 15, 30, 50, 100, 500, 1000)
mu <- 10
stadev <- 2
true_variance <- stadev^2

# Mendefinisikan penduga
W1 <- function(x) { 
  sum((x - mean(x))^2) / (length(x) - 1)
}
W2 <- function(x) { 
  (median(abs(x - median(x))))^2
}

# Ukuran performa
bias <- function(estimate, parameter) {
  mean(estimate) - parameter
}
relative_bias <- function(estimate, parameter) {
  (mean(estimate) - parameter) / parameter
}
empirical_variance <- function(estimate) { var(estimate) }

# Pembangkitan data
set.seed(1082)
gen_data <- map(n_obs, function(i) {
  x <- replicate(n = n_sim,
                 expr = rnorm(i, mean = mu, sd = stadev)
  ) %>% as.data.frame()
  return(x)
})

# Fungsi untuk menambahkan pencilan
add_outliers <- function(data) {
  data_with_outliers <- data
  for (i in seq_along(data)) {
    median_data <- median(data[[i]])
    iqr_data <- IQR(data[[i]])
    medium <- median_data + 2 * iqr_data
    ekstrim <- median_data + 5 * iqr_data
    data_with_outliers[[i]][1] <- medium  # Menambahkan pencilan medium
    data_with_outliers[[i]][2] <- ekstrim  # Menambahkan pencilan ekstrim
  }
  return(data_with_outliers)
}

# Menambahkan pencilan pada setiap dataset
gen_data_with_outliers <- map(gen_data, add_outliers)

# Menghitung W1 dan W2 untuk setiap ukuran sampel setelah menambahkan pencilan
W1_results <- map(gen_data_with_outliers, function(data) {
  apply(data, 2, W1)
})

W2_results <- map(gen_data_with_outliers, function(data) {
  apply(data, 2, W2)
})

# Menghitung ukuran performa untuk W1 dan W2
W1_performance <- map(W1_results, function(estimate) {
  tibble(
    bias = bias(estimate, true_variance),
    relative_bias = relative_bias(estimate, true_variance),
    empirical_variance = empirical_variance(estimate)
  )
})

W2_performance <- map(W2_results, function(estimate) {
  tibble(
    bias = bias(estimate, true_variance),
    relative_bias = relative_bias(estimate, true_variance),
    empirical_variance = empirical_variance(estimate)
  )
})

# Menggabungkan hasil performa dalam satu data frame
performance_results_W1 <- map2_df(n_obs, W1_performance, ~ mutate(.y, n_obs = .x, estimator = "W1"))
performance_results_W2 <- map2_df(n_obs, W2_performance, ~ mutate(.y, n_obs = .x, estimator = "W2"))
performance_results <- bind_rows(performance_results_W1, performance_results_W2)

# Menampilkan hasil dalam bentuk tabel
kable(performance_results, caption = "Ukuran Performa Estimator W1 dan W2 dalam Pendugaan Ragam") %>% 
  kable_styling(full_width = FALSE)

Ukuran Performa Estimator W1 dan W2 dalam Pendugaan Ragam

bias	relative_bias	empirical_variance	n_obs	estimator
24.7740290	6.1935072	1027.5418658	5	W1
11.4645706	2.8661426	59.5891842	15	W1
6.1336539	1.5334135	10.4004846	30	W1
3.8190349	0.9547587	3.5302453	50	W1
2.0622689	0.5155672	0.8819524	100	W1
0.4008482	0.1002120	0.0808392	500	W1
0.1908640	0.0477160	0.0365599	1000	W1
3.8613382	0.9653346	82.0984459	5	W2
-1.2955102	-0.3238776	2.6835843	15	W2
-1.8454408	-0.4613602	0.7179872	30	W2
-2.0258780	-0.5064695	0.3907648	50	W2
-2.0846609	-0.5211652	0.1866237	100	W2
-2.1603259	-0.5400815	0.0369479	500	W2
-2.1787309	-0.5446827	0.0184469	1000	W2