# Packages
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.1 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ ggplot2 3.4.2 ✔ tibble 3.2.1
## ✔ lubridate 1.9.2 ✔ tidyr 1.3.0
## ✔ purrr 1.0.1
## ── 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)
##
## 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))))}
# Ukuran performa
bias <- function(estimate, parameter) {
mean(estimate) - parameter
}
relative_bias <- function(estimate, parameter) {
(mean(estimate) - parameter) / parameter
}
empirical_variance <- function(estimate) { var(estimate) }
Data Bangkitan dengan Pencilan Medium
# Pembangkitan data
set.seed(1082)
gen_data <- map(n_obs-1, 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 %>% as.list()
for (i in seq_along(data_with_outliers)) {
median_data <- median(data_with_outliers[[i]])
iqr_data <- IQR(data_with_outliers[[i]])
outlier_value <- median_data + 2 * iqr_data
data_with_outliers[[i]] <- c(data_with_outliers[[i]], outlier_value) # Menambahkan pencilan
}
data_with_outliers <- as.data.frame(data_with_outliers)
return(data_with_outliers)
}
# Menambahkan pencilan pada setiap dataset
gen_data_with_outliersmed <- map(gen_data, add_outliers)
# Menghitung W1 dan W2 untuk
W1_results <- map(gen_data_with_outliersmed, function(data) {
apply(data, 2, W1)})
W2_results <- map(gen_data_with_outliersmed, 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 <- cbind(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
|
bias
|
relative_bias
|
empirical_variance
|
n_obs
|
estimator
|
|
2.9476339
|
0.7369085
|
40.2054279
|
5
|
W1
|
-2.349078
|
-0.5872695
|
0.9360484
|
5
|
W2
|
|
1.4306983
|
0.3576746
|
4.8322440
|
15
|
W1
|
-2.593501
|
-0.6483751
|
0.1857778
|
15
|
W2
|
|
0.7661446
|
0.1915361
|
1.6656363
|
30
|
W1
|
-2.625190
|
-0.6562976
|
0.0771155
|
30
|
W2
|
|
0.4716430
|
0.1179107
|
0.9004991
|
50
|
W1
|
-2.644421
|
-0.6611053
|
0.0509759
|
50
|
W2
|
|
0.2992956
|
0.0748239
|
0.3945772
|
100
|
W1
|
-2.640270
|
-0.6600676
|
0.0238734
|
100
|
W2
|
|
0.0460549
|
0.0115137
|
0.0648254
|
500
|
W1
|
-2.648537
|
-0.6621343
|
0.0047764
|
500
|
W2
|
|
0.0142375
|
0.0035594
|
0.0329000
|
1000
|
W1
|
-2.652585
|
-0.6631462
|
0.0023516
|
1000
|
W2
|
Data Bangkitan dengan Pencilan Ekstrim
# Pembangkitan data
set.seed(1082)
gen_data <- map(n_obs-1, 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 %>% as.list()
for (i in seq_along(data_with_outliers)) {
median_data <- median(data_with_outliers[[i]])
iqr_data <- IQR(data_with_outliers[[i]])
outlier_value <- median_data + 5 * iqr_data
data_with_outliers[[i]] <- c(data_with_outliers[[i]], outlier_value) # Menambahkan pencilan
}
data_with_outliers <- as.data.frame(data_with_outliers)
return(data_with_outliers)
}
# Menambahkan pencilan pada setiap dataset
gen_data_with_outlierseks <- map(gen_data, add_outliers)
# Menghitung W1 dan W2
W1_results <- map(gen_data_with_outlierseks, function(data) {
apply(data, 2, W1)})
W2_results <- map(gen_data_with_outlierseks, 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 <- cbind(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
|
bias
|
relative_bias
|
empirical_variance
|
n_obs
|
estimator
|
|
22.8934361
|
5.7233590
|
782.5180483
|
5
|
W1
|
-2.349078
|
-0.5872695
|
0.9360484
|
5
|
W2
|
|
10.4882509
|
2.6220627
|
52.6814024
|
15
|
W1
|
-2.593501
|
-0.6483751
|
0.1857778
|
15
|
W2
|
|
5.5213740
|
1.3803435
|
9.7628038
|
30
|
W1
|
-2.625190
|
-0.6562976
|
0.0771155
|
30
|
W2
|
|
3.4270395
|
0.8567599
|
3.3153056
|
50
|
W1
|
-2.644421
|
-0.6611053
|
0.0509759
|
50
|
W2
|
|
1.8134835
|
0.4533709
|
0.8027089
|
100
|
W1
|
-2.640270
|
-0.6600676
|
0.0238734
|
100
|
W2
|
|
0.3515677
|
0.0878919
|
0.0751511
|
500
|
W1
|
-2.648537
|
-0.6621343
|
0.0047764
|
500
|
W2
|
|
0.1663511
|
0.0415878
|
0.0354857
|
1000
|
W1
|
-2.652585
|
-0.6631462
|
0.0023516
|
1000
|
W2
|