1. Ketakbiasan Penduga

Jenis Penduga : Nilai tengah -> rataan, median Ragam -> ragam contoh (n-1), ragam contoh(n-1)

Normal(4,16)

set.seed(117)
n = 10  # ukuran sampel
k = 1000 # jumlah ulangan
N = 20 # ukuran populasi
populasi = rnorm(N, 4, sqrt(16)) # rnorm(n, mean, sd)
hist(populasi)

Penduga rataan

miu = mean(populasi) 

# Mengambil sampel dan menghitung rata-rata sampel
sampel = matrix(NA, k, n)
for (i in 1:k) sampel[i,] = sample(populasi, n)
xbar = apply(sampel, 1, mean)
E.xbar = mean(xbar) # estimasi harapan dari xbar

# Output hasil
cat("Mean Populasi:", miu, "\n")

## Mean Populasi: 3.389474

cat("E(x̄) dari Sampel:", E.xbar, "\n")

## E(x̄) dari Sampel: 3.379834

Terlihat bahwa nilai E(Xbar) mendekati nilai miu yaitu nilai rata rat populasi sehingga pada kasus ini menyatakan E(xbar) tidak bias terhadapa nilai miu.

Penduga Ragam

sigma2 = var(populasi) * (length(populasi) - 1) / length(populasi) # varians populasi

# Mengambil sampel 
sampel = matrix(NA, k, n)
for (i in 1:k) sampel[i,] = sample(populasi, n)

# Menghitung ragam (n)
s2_biased = apply(sampel, 1, function(x) var(x) * (length(x) - 1) / length(x)) 

E.s2_biased = mean(s2_biased) 

# Mengambil sampel dan menghitung varians dengan koreksi (tak bias)
s2_unbiased = apply(sampel, 1, var) # Ragam dengan koreksi (n-1)
E.s2_unbiased = mean(s2_unbiased) # Ekspektasi varians dengan koreksi

# Output hasil
cat("Varians Populasi:", sigma2, "\n")

## Varians Populasi: 17.16209

cat("E(s^2) tanpa koreksi (bias):", E.s2_biased, "\n")

## E(s^2) tanpa koreksi (bias): 16.19388

cat("E(s^2) dengan koreksi (tak bias):", E.s2_unbiased, "\n")

## E(s^2) dengan koreksi (tak bias): 17.9932

Karena populasinya hanya 20 angka, maka varians populasi (sigma²) yang rentan terhadap fluktuasi (noise). maka terlihat hasilnya tidak terlalu jauh

Seragam(5,10)

set.seed(121)
n = 10  # ukuran sampel
k = 1000 # jumlah ulangan
N = 20 # ukuran populasi
populasi = runif(N, 5, 10) # runif(n,min,max)
hist(populasi)

Penduga Rataan

miu = mean(populasi) 

# Mengambil sampel dan menghitung rata-rata sampel
sampel = matrix(NA, k, n)
for (i in 1:k) sampel[i,] = sample(populasi, n)
xbar = apply(sampel, 1, mean)
E.xbar = mean(xbar) # estimasi harapan dari xbar

# Output hasil
cat("Mean Populasi:", miu, "\n")

## Mean Populasi: 7.691463

cat("E(x̄) dari Sampel:", E.xbar, "\n")

## E(x̄) dari Sampel: 7.69513

Terlihat bahwa nilai E(Xbar) mendekati nilai miu yaitu nilai rata rat populasi sehingga pada kasus ini menyatakan E(xbar) tidak bias terhadapa nilai miu.

E(X) untuk sebaran seragam adalah (b+a)/2 yaitu (15)/2 = 7.5

Penduga Ragam

sigma2 = var(populasi) * (length(populasi) - 1) / length(populasi) # varians populasi

# Mengambil sampel 
sampel = matrix(NA, k, n)
for (i in 1:k) sampel[i,] = sample(populasi, n)

# Menghitung ragam (n)
s2_biased = apply(sampel, 1, function(x) var(x) * (length(x) - 1) / length(x)) 

E.s2_biased = mean(s2_biased) 

# Mengambil sampel dan menghitung varians dengan koreksi (tak bias)
s2_unbiased = apply(sampel, 1, var) # Ragam dengan koreksi (n-1)
E.s2_unbiased = mean(s2_unbiased) # Ekspektasi varians dengan koreksi

# Output hasil
cat("Varians Populasi:", sigma2, "\n")

## Varians Populasi: 1.501897

cat("E(s^2) tanpa koreksi (bias):", E.s2_biased, "\n")

## E(s^2) tanpa koreksi (bias): 1.419181

cat("E(s^2) dengan koreksi (tak bias):", E.s2_unbiased, "\n")

## E(s^2) dengan koreksi (tak bias): 1.576868

Var(X) untuk sebaran seragam adalah (b-a)^2/12 = 25/12 = 2.0833

Tidak seperti distribusi normal, distribusi uniform tidak punya “bentuk simetris tajam” di tengah.

Untuk distribusi yang tidak normal, performa estimator unbiased bisa lebih fluktuatif, terutama jika n kecil (sampel = 10).

Kesimpulan based on gpt : 1. Populasi kecil → varians populasi acak & tidak stabil.

Distribusi uniform → varians data lebih merata & tidak terpusat.
Estimator unbiased memang bisa lebih variatif dalam kasus kecil & non-normal.

Eksponensial(1)

set.seed(121)
n = 10  # ukuran sampel
k = 1000 # jumlah ulangan
N = 20 # ukuran populasi
populasi = rexp(N,1)
hist(populasi)

Terlihat bahwa nilai E(Xbar) mendekati nilai miu yaitu nilai rata rat populasi sehingga pada kasus ini menyatakan E(xbar) tidak bias terhadapa nilai miu.

Penduga Rataan

miu = mean(populasi) 

# Mengambil sampel dan menghitung rata-rata sampel
sampel = matrix(NA, k, n)
for (i in 1:k) sampel[i,] = sample(populasi, n)
xbar = apply(sampel, 1, mean)
E.xbar = mean(xbar) # estimasi harapan dari xbar

# Output hasil
cat("Mean Populasi:", miu, "\n")

## Mean Populasi: 0.8305081

cat("E(x̄) dari Sampel:", E.xbar, "\n")

## E(x̄) dari Sampel: 0.8298502

E(X) untuk sebaran eksponensial adalah 1/lambda = 1

Penduga Ragam

sigma2 = var(populasi) * (length(populasi) - 1) / length(populasi) # varians populasi

# Mengambil sampel 
sampel = matrix(NA, k, n)
for (i in 1:k) sampel[i,] = sample(populasi, n)

# Menghitung ragam (n)
s2_biased = apply(sampel, 1, function(x) var(x) * (length(x) - 1) / length(x)) 

E.s2_biased = mean(s2_biased) 

# Mengambil sampel dan menghitung varians dengan koreksi (tak bias)
s2_unbiased = apply(sampel, 1, var) # Ragam dengan koreksi (n-1)
E.s2_unbiased = mean(s2_unbiased) # Ekspektasi varians dengan koreksi

# Output hasil
cat("Varians Populasi:", sigma2, "\n")

## Varians Populasi: 0.4810118

cat("E(s^2) tanpa koreksi (bias):", E.s2_biased, "\n")

## E(s^2) tanpa koreksi (bias): 0.4607114

cat("E(s^2) dengan koreksi (tak bias):", E.s2_unbiased, "\n")

## E(s^2) dengan koreksi (tak bias): 0.5119016

E(X) untuk sebaran eksponensial adalah 1/lambda^2 = 1

Estimator tak bias s²_unbiased seharusnya menuju ke varian distribusi sebenarnya (1), bukan ke varian dari 20 angka acak itu.

2. Penduga selang kepercayaan

Normal

n=10

n    = 10
k    = 100 # ulangan
alpha  = 0.05
mu   = 50
std  = 10

set.seed(111)
sampel = matrix(rnorm(n*k,mu,std),k)
xbar = apply(sampel,1,mean)
s    = apply(sampel,1,sd)
SE   = s/sqrt(n)
z    = qnorm(1-alpha/2)
SK   = (xbar-z*SE < mu & mu < xbar+z*SE)

propz = sum(SK)/k # proporsi banyaknya SK yang memuat mu
propz

## [1] 0.93

matplot(rbind(xbar-z*SE, xbar+z*SE), rbind(1:k,1:k), col=ifelse(SK,"blue","red"), type = "l", lty = 1, 
        main='Selang Kepercayaan', xlab='SK', ylab='Banyak Ulangan')
abline(v=mu)

cat("Proporsi SK yang memuat mu (σ diketahui, pakai Z):", propz, "\n")

## Proporsi SK yang memuat mu (σ diketahui, pakai Z): 0.93

# --- 2. Selang kepercayaan saat σ TIDAK diketahui (pakai s dan distribusi t)
s <- apply(sampel, 1, sd)
SE_unknown <- s / sqrt(n)
t <- qt(1 - alpha/2, df = n - 1)
SK_t <- (xbar - t * SE_unknown < mu) & (mu < xbar + t * SE_unknown)
prop_t <- mean(SK_t)
matplot(rbind(xbar - t * SE_unknown,xbar + t * SE_unknown), rbind(1:k,1:k), col=ifelse(SK,"blue","red"), type = "l", lty = 1, 
        main='Selang Kepercayaan', xlab='SK', ylab='Banyak Ulangan')
abline(v=mu)

cat("Proporsi SK yang memuat mu (σ tidak diketahui, pakai t):", round(prop_t, 4), "\n")

## Proporsi SK yang memuat mu (σ tidak diketahui, pakai t): 0.94

n=30

n    = 30
k    = 100 # ulangan
alpha  = 0.05
mu   = 0
std  = 1

set.seed(111)
sampel = matrix(rnorm(n*k,mu,std),k)
xbar = apply(sampel,1,mean)
s    = apply(sampel,1,sd)
SE   = s/sqrt(n)
z    = qnorm(1-alpha/2)
SK   = (xbar-z*SE < mu & mu < xbar+z*SE)

sum(SK)/k # proporsi banyaknya SK yang memuat mu

## [1] 0.94

## [1] 0.94
matplot(rbind(xbar-z*SE, xbar+z*SE), rbind(1:k,1:k), col=ifelse(SK,"blue","red"), type = "l", lty = 1, 
        main='Selang Kepercayaan', xlab='SK', ylab='Banyak Ulangan')
abline(v=mu)

n = 100

n    = 100
k    = 100 # ulangan
alpha  = 0.05
mu   = 0
std  = 1

set.seed(111)
sampel = matrix(rnorm(n*k,mu,std),k)
xbar = apply(sampel,1,mean)
s    = apply(sampel,1,sd)
SE   = s/sqrt(n)
z    = qnorm(1-alpha/2)
SK   = (xbar-z*SE < mu & mu < xbar+z*SE)

sum(SK)/k # proporsi banyaknya SK yang memuat mu

## [1] 0.95

## [1] 0.94
matplot(rbind(xbar-z*SE, xbar+z*SE), rbind(1:k,1:k), col=ifelse(SK,"blue","red"), type = "l", lty = 1, 
        main='Selang Kepercayaan', xlab='SK', ylab='Banyak Ulangan')
abline(v=mu)

Seragam(0,1)

Uniform <- function(n, reps = 1000) {
  means <- replicate(reps, mean(runif(n, 0, 1)))
  return(means)
}

Membuat fungsi untuk membangkitkan data dengan distribusi uniform(0,1) yang dilakukan dengan pengulangan 1000x dan diambil rataannya

set.seed(123)
sample_sizes <- c(10, 30, 100)
par(mfrow = c(1, 3), mar = c(4, 4, 2, 1))
for (n in sample_sizes) {
  means <- Uniform(n)
  hist(
    means,
    breaks = 30,
    probability = TRUE,
    col = "lightgreen",
    main = paste("Uniform (n=", n, ")", sep = ""),
    xlab = "Rata-rata Sampel",
    ylab = "Kerapatan"
  )
}

n = 10

set.seed(123)
n <- 10
k <- 100       
B <- 1000
mu <- 0.5 # nilai tengah dari sebaran uni(0,1)
alpha <- 0.05

CI_lower <- numeric(k)
CI_upper <- numeric(k)
contains_mu <- logical(k)

for (i in 1:k) {
  sampel <- runif(n, 0,1)
  boots <- replicate(B, mean(sample(sampel, n, replace = TRUE)))
  CI <- quantile(boots, probs = c(alpha/2, 1 - alpha/2))
  CI_lower[i] <- CI[1]
  CI_upper[i] <- CI[2]
  contains_mu[i] <- CI[1] < mu && mu < CI[2]
}

# Plot CI menggunakan matplot
matplot(rbind(CI_lower, CI_upper), 
        rbind(1:k, 1:k), 
        type = "l", lty = 1, 
        col = ifelse(contains_mu, "blue", "red"),
        xlab = "Interval", ylab = "Simulasi ke-", 
        main = "Visualisasi Selang Kepercayaan Bootstrap")
abline(v = mu, col = "black", lty = 2)  # garis vertikal di mu

# Proporsi CI yang memuat mu
prop <- mean(contains_mu)
cat("Proporsi Bootstrap CI yang memuat mu =", round(prop, 4), "\n")

## Proporsi Bootstrap CI yang memuat mu = 0.9

n = 30

set.seed(123)
n <- 30
k <- 100       
B <- 1000
mu <- 0.5
alpha <- 0.05

CI_lower <- numeric(k)
CI_upper <- numeric(k)
contains_mu <- logical(k)

for (i in 1:k) {
  sampel <- runif(n, 0,1)
  boots <- replicate(B, mean(sample(sampel, n, replace = TRUE)))
  CI <- quantile(boots, probs = c(alpha/2, 1 - alpha/2))
  CI_lower[i] <- CI[1]
  CI_upper[i] <- CI[2]
  contains_mu[i] <- CI[1] < mu && mu < CI[2]
}

# Plot CI menggunakan matplot
matplot(rbind(CI_lower, CI_upper), 
        rbind(1:k, 1:k), 
        type = "l", lty = 1, 
        col = ifelse(contains_mu, "blue", "red"),
        xlab = "Interval", ylab = "Simulasi ke-", 
        main = "Visualisasi Selang Kepercayaan Bootstrap")
abline(v = mu, col = "black", lty = 2)  # garis vertikal di mu

# Proporsi CI yang memuat mu
prop <- mean(contains_mu)
cat("Proporsi Bootstrap CI yang memuat mu =", round(prop, 4), "\n")

## Proporsi Bootstrap CI yang memuat mu = 0.95

n = 100

set.seed(123)
n <- 100
k <- 100       
B <- 1000
mu <- 0.5
alpha <- 0.05

CI_lower <- numeric(k)
CI_upper <- numeric(k)
contains_mu <- logical(k)

for (i in 1:k) {
  sampel <- runif(n, 0,1)
  boots <- replicate(B, mean(sample(sampel, n, replace = TRUE)))
  CI <- quantile(boots, probs = c(alpha/2, 1 - alpha/2))
  CI_lower[i] <- CI[1]
  CI_upper[i] <- CI[2]
  contains_mu[i] <- CI[1] < mu && mu < CI[2]
}

# Plot CI menggunakan matplot
matplot(rbind(CI_lower, CI_upper), 
        rbind(1:k, 1:k), 
        type = "l", lty = 1, 
        col = ifelse(contains_mu, "blue", "red"),
        xlab = "Interval", ylab = "Simulasi ke-", 
        main = "Visualisasi Selang Kepercayaan Bootstrap")
abline(v = mu, col = "black", lty = 2)  # garis vertikal di mu

# Proporsi CI yang memuat mu
prop <- mean(contains_mu)
cat("Proporsi Bootstrap CI yang memuat mu =", round(prop, 4), "\n")

## Proporsi Bootstrap CI yang memuat mu = 0.95

Eksponensial(1)

Expo <- function(n, reps = 1000) {
  means <- replicate(reps, mean(rexp(n, rate = 1)))
  return(means)
}

set.seed(123)
par(mfrow = c(1, 3), mar = c(4, 4, 2, 1))
for (n in sample_sizes) {
  means <- Expo(n)
  hist(
    means,
    breaks = 30,
    probability = TRUE,
    col = "lightblue",
    main = paste("Eksponensial (n=", n, ")", sep = ""),
    xlab = "Rata-rata Sampel",
    ylab = "Kerapatan"
  )
}

n = 10

set.seed(123)
n <- 10
k <- 100       # jumlah simulasi
B <- 1000      # jumlah bootstrap
mu <- 0.5      # mean dari Uniform(0,1)
alpha <- 0.05

CI_lower <- numeric(k)
CI_upper <- numeric(k)
contains_mu <- logical(k)

for (i in 1:k) {
  sampel <- runif(n, min = 0, max = 1)
  boots <- replicate(B, mean(sample(sampel, n, replace = TRUE)))
  CI <- quantile(boots, probs = c(alpha/2, 1 - alpha/2))
  CI_lower[i] <- CI[1]
  CI_upper[i] <- CI[2]
  contains_mu[i] <- CI[1] < mu && mu < CI[2]
}

# Visualisasi selang kepercayaan
matplot(rbind(CI_lower, CI_upper),
        rbind(1:k, 1:k),
        type = "l", lty = 1,
        col = ifelse(contains_mu, "blue", "red"),
        xlab = "Interval", ylab = "Simulasi ke-",
        main = "Selang Kepercayaan Bootstrap - Uniform(0,1)")
abline(v = mu, col = "black", lty = 2)

# Proporsi yang memuat mu
prop <- mean(contains_mu)
cat("Proporsi Bootstrap CI yang memuat mu =", round(prop, 4), "\n")

## Proporsi Bootstrap CI yang memuat mu = 0.9

n = 30

set.seed(123)
n <- 30
k <- 100       
B <- 1000
mu <- 1
alpha <- 0.05

CI_lower <- numeric(k)
CI_upper <- numeric(k)
contains_mu <- logical(k)

for (i in 1:k) {
  sampel <- rexp(n, rate = 1)
  boots <- replicate(B, mean(sample(sampel, n, replace = TRUE)))
  CI <- quantile(boots, probs = c(alpha/2, 1 - alpha/2))
  CI_lower[i] <- CI[1]
  CI_upper[i] <- CI[2]
  contains_mu[i] <- CI[1] < mu && mu < CI[2]
}

# Plot CI menggunakan matplot
matplot(rbind(CI_lower, CI_upper), 
        rbind(1:k, 1:k), 
        type = "l", lty = 1, 
        col = ifelse(contains_mu, "blue", "red"),
        xlab = "Interval", ylab = "Simulasi ke-", 
        main = "Visualisasi Selang Kepercayaan Bootstrap")
abline(v = mu, col = "black", lty = 2)  # garis vertikal di mu

# Proporsi CI yang memuat mu
prop <- mean(contains_mu)
cat("Proporsi Bootstrap CI yang memuat mu =", round(prop, 4), "\n")

## Proporsi Bootstrap CI yang memuat mu = 0.91

n = 100

set.seed(123)
n <- 100
k <- 100       
B <- 1000
mu <- 1
alpha <- 0.05

CI_lower <- numeric(k)
CI_upper <- numeric(k)
contains_mu <- logical(k)

for (i in 1:k) {
  sampel <- rexp(n, rate = 1)
  boots <- replicate(B, mean(sample(sampel, n, replace = TRUE)))
  CI <- quantile(boots, probs = c(alpha/2, 1 - alpha/2))
  CI_lower[i] <- CI[1]
  CI_upper[i] <- CI[2]
  contains_mu[i] <- CI[1] < mu && mu < CI[2]
}

# Plot CI menggunakan matplot
matplot(rbind(CI_lower, CI_upper), 
        rbind(1:k, 1:k), 
        type = "l", lty = 1, 
        col = ifelse(contains_mu, "blue", "red"),
        xlab = "Interval", ylab = "Simulasi ke-", 
        main = "Selang Kepercayaan Bootstrap Exp(1)")
abline(v = mu, col = "black", lty = 2)  # garis vertikal di mu

# Proporsi CI yang memuat mu
prop <- mean(contains_mu)
cat("Proporsi Bootstrap CI yang memuat mu =", round(prop, 4), "\n")

## Proporsi Bootstrap CI yang memuat mu = 0.98

3. Batas atas dan batas bawah selang kepercayaan

Simetris

# Parameter awal
mu <- 50
sigma <- 10
n <- 100
se <- sigma / sqrt(n)
alpha <- 0.05
z_crit <- qnorm(1 - alpha/2)

# Jumlah simulasi
N <- 100

# Vektor untuk menyimpan hasil
lower_sym <- numeric(N)
upper_sym <- numeric(N)
contains_mu <- logical(N)

set.seed(123)

for (i in 1:N) {
  x <- rnorm(n, mean = mu, sd = sigma)
  x_bar <- mean(x)
  
  lower <- x_bar - z_crit * se
  upper <- x_bar + z_crit * se
  
  lower_sym[i] <- lower
  upper_sym[i] <- upper
  contains_mu[i] <- (mu >= lower & mu <= upper)
}

# Hitung proporsi selang yang mengandung mu
coverage_rate <- mean(contains_mu)
cat("Proporsi selang yang mengandung mu:", coverage_rate, "\n")

## Proporsi selang yang mengandung mu: 0.97

cat("Selang simetris  : [", lower, ",", upper, "]\n")

## Selang simetris  : [ 47.69487 , 51.6148 ]

cat("Jarak ke bawah   :", x_bar - lower, "\n")

## Jarak ke bawah   : 1.959964

cat("Jarak ke atas    :", upper - x_bar, "\n\n")

## Jarak ke atas    : 1.959964

Tidak simetris

# Parameter awal
mu <- 50
sigma <- 10
n <- 100
se <- sigma / sqrt(n)
alpha_lower <- 0.01
alpha_upper <- 0.04

# Hitung z-score untuk batas bawah dan atas
z_lower <- qnorm(alpha_lower)
z_upper <- qnorm(1 - alpha_upper)

# Jumlah simulasi
N <- 100

# Vektor untuk menyimpan hasil
lower_bounds <- numeric(N)
upper_bounds <- numeric(N)
contains_mu <- logical(N)

set.seed(123)

for (i in 1:N) {
  x <- rnorm(n, mean = mu, sd = sigma)
  x_bar <- mean(x)
  
  lower <- x_bar + z_lower * se
  upper <- x_bar + z_upper * se
  
  lower_bounds[i] <- lower
  upper_bounds[i] <- upper
  contains_mu[i] <- (mu >= lower & mu <= upper)
}

# Hitung proporsi selang yang mengandung mu
coverage_rate <- mean(contains_mu)
cat("Proporsi selang yang mengandung mu (tidak simetris):", coverage_rate, "\n")

## Proporsi selang yang mengandung mu (tidak simetris): 0.96

cat("Selang simetris  : [", lower, ",", upper, "]\n")

## Selang simetris  : [ 47.32849 , 51.40552 ]

cat("Jarak ke bawah   :", x_bar - lower, "\n")

## Jarak ke bawah   : 2.326348

cat("Jarak ke atas    :", upper - x_bar, "\n\n")

## Jarak ke atas    : 1.750686

Tugas MSR (1)

Muhammad Aqil Siradj

2025-05-04

1. Ketakbiasan Penduga

Normal(4,16)

Penduga rataan

Penduga Ragam

Seragam(5,10)

Penduga Rataan

Penduga Ragam

Eksponensial(1)

Penduga Rataan

Penduga Ragam

2. Penduga selang kepercayaan

Normal

n=10

n=30

n = 100

Seragam(0,1)

n = 10

n = 30

n = 100

Eksponensial(1)

n = 10

n = 30

n = 100

3. Batas atas dan batas bawah selang kepercayaan

Simetris

Tidak simetris