# Pengaplikasian pada Software R

par(mfrow=c(3,1))
library(probs)
## Warning: package 'probs' was built under R version 4.4.3
## 
## Attaching package: 'probs'
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, union
# a). Hampiran Normal terhadap Distribusi Geometrik
set.seed(42)
# Membentuk populasi distribusi geometrik
populasi <- rgeom(20, 0.19)

# Ukuran sampel n = 3
n1 <- 3
contoh_geo1 <- urnsamples(populasi,
                          size = n1,
                          replace = FALSE,
                          ordered = FALSE)

mean_geo1 <- matrix(apply(contoh_geo1, 1, mean))

# Ukuran sampel n = 6
n2 <- 6
contoh_geo2 <- urnsamples(populasi,
                          size = n2,
                          replace = FALSE,
                          ordered = FALSE)

mean_geo2 <- matrix(apply(contoh_geo2, 1, mean))

# Histogram hasil rataan
hist(mean_geo1,
     main = "Hampiran Normal Terhadap Geometrik (n = 3)",
     xlab = "xbar",
     col = "orange")

hist(mean_geo2,
     main = "Hampiran Normal Terhadap Geometrik (n = 6)",
     xlab = "xbar",
     col = "lightgreen")

# b). Hampiran Normal Terhadap Eksponensial
set.seed(42)
# Membentuk populasi dari distribusi eksponensial
populasi = rexp(20)

# Menentukan ukuran sampel pertama
n1 = 3
contoh_exp1 = urnsamples(populasi,
                         size = n1,
                         replace = FALSE,
                         ordered = FALSE)
mean_exp1 = matrix(apply(contoh_exp1, 1, mean))

# Menentukan ukuran sampel kedua
n2 = 6
contoh_exp2 = urnsamples(populasi,
                         size = n2,
                         replace = FALSE,
                         ordered = FALSE)
mean_exp2 = matrix(apply(contoh_exp2, 1, mean))

# Histogram distribusi rataan untuk n = 3
hist(mean_exp1,
     main = "Hampiran Normal Terhadap Eksponensial (n = 3)",
     xlab = "xbar",
     col = "yellow")

# Histogram distribusi rataan untuk n = 6
hist(mean_exp2,
     main = "Hampiran Normal Terhadap Eksponensial (n = 6)",
     xlab = "xbar",
     col = "gold")

# c). Hampiran Normal Terhadap Seragam
set.seed(42)
# Membentuk populasi distribusi seragam
populasi = runif(20)

# Ukuran sampel pertama
n1 = 3
contoh_unif1 = urnsamples(populasi,
                          size = n1,
                          replace = FALSE,
                          ordered = FALSE)
mean_unif1 = matrix(apply(contoh_unif1, 1, mean))

# Ukuran sampel kedua
n2 = 6
contoh_unif2 = urnsamples(populasi,
                          size = n2,
                          replace = FALSE,
                          ordered = FALSE)
mean_unif2 = matrix(apply(contoh_unif2, 1, mean))

# Histogram rataan sampel n = 3
hist(mean_unif1,
     main = "Hampiran Normal Terhadap Seragam (n = 3)",
     xlab = "xbar",
     col = "blue")

# Histogram rataan sampel n = 6
hist(mean_unif2,
     main = "Hampiran Normal Terhadap Seragam (n = 6)",
     xlab = "xbar",
     col = "lightgreen")

#Kesimpulan : Berdasarkan simulasi menggunakan distribusi geometrik, eksponensial, dan seragam, diperoleh bahwa semakin besar ukuran sampel (n), maka distribusi rataan sampel akan semakin mendekati distribusi normal. Hal ini terlihat pada histogram rataan sampel yang semakin membentuk pola menyerupai kurva normal ketika ukuran sampel diperbesar dari n = 3 menjadi n = 6. Hasil tersebut sesuai dengan konsep Central Limit Theorem (CLT), yaitu distribusi rataan sampel cenderung berdistribusi normal meskipun data awal berasal dari distribusi yang berbeda.
#Sebaran Percontohan Sebaran Normal

set.seed(42)
populasi     = rnorm(20,5,sqrt(12))

n1           = 3
contoh_norm1 = urnsamples(populasi, size = 3, replace = F, ordered = F)
head(contoh_norm1)
##         X1       X2        X3
## 1 9.749139 3.043828  6.257914
## 2 9.749139 3.043828  7.192300
## 3 9.749139 3.043828  6.400427
## 4 9.749139 3.043828  4.632374
## 5 9.749139 3.043828 10.236066
## 6 9.749139 3.043828  4.672091
# Menghitung rata-rata dan varians rataan sampel untuk n = 3
mean_norm1   = matrix(apply(contoh_norm1, 1, mean))
mean_xbar1   = mean(mean_norm1)
var_xbar1    = var(mean_norm1)

# Menentukan ukuran sampel kedua
n3           = 6

# Mengambil seluruh kemungkinan sampel ukuran 6
contoh_norm3 = urnsamples(populasi,
                          size = 6,
                          replace = FALSE,
                          ordered = FALSE)

# Menghitung rata-rata dan varians rataan sampel
mean_norm3   = matrix(apply(contoh_norm3, 1, mean))
mean_xbar3   = mean(mean_norm3)
var_xbar3    = var(mean_norm3)

# Histogram distribusi rataan sampel untuk n = 3
hist(mean_norm1,
     main = "(n = 3)",
     xlab = "xbar",
     col = "lightblue")

# Histogram distribusi rataan sampel untuk n = 6
hist(mean_norm3,
     main = "(n = 6)",
     xlab = "xbar",
     col = "lightpink")

# Membuat tabel hasil mean dan varians
hasil = data.frame(
  "." = c("mean", "varian"),
  "Populasi" = c(5, 12),
  "n=3" = c(mean_xbar1, var_xbar1),
  "n=6" = c(mean_xbar3, var_xbar3)
)

# Menampilkan hasil
hasil
##        . Populasi      n.3      n.6
## 1   mean        5 5.664830 5.664830
## 2 varian       12 5.863321 2.412253
#Kesimpulan :
#- Berdasarkan output yang diperoleh, distribusi rataan sampel yang diambil dari populasi normal dengan mean = 5 dan varians = 12 menunjukkan bahwa semakin besar ukuran sampel, maka mean dari rataan sampel akan semakin mendekati mean populasi (μ), sedangkan varians rataan sampel akan semakin kecil dan mendekati σ²/n.

#- Berdasarkan histogram distribusi rataan sampel, terlihat bahwa ketika ukuran sampel meningkat dari n = 3 menjadi n = 6, bentuk histogram semakin menyerupai distribusi normal. Hal ini menunjukkan bahwa distribusi rataan sampel menjadi lebih stabil dan lebih terpusat di sekitar mean populasi.
# Penerapan ketakbiasan penduga (mean)

#POPULASI TERHINGGA
library(probs)
# 1. Sebaran Normal
set.seed(42)

populasi1      = rnorm(20)
mean_pop1      = mean(populasi1)

sampel_normal1 = urnsamples(populasi1,
                            size = 5,
                            replace = FALSE,
                            ordered = FALSE)

mean_normal1   = matrix(apply(sampel_normal1, 1, mean))
median_normal1 = matrix(apply(sampel_normal1, 1, median))

harapan_mean_norm1   = mean(mean_normal1)
harapan_median_norm1 = mean(median_normal1)

# 2. Sebaran Eksponensial
set.seed(42)

populasi2    = rexp(20)
mean_pop2    = mean(populasi2)

sampel_exp1  = urnsamples(populasi2,
                          size = 5,
                          replace = FALSE,
                          ordered = FALSE)

mean_exp1    = matrix(apply(sampel_exp1, 1, mean))
median_exp1  = matrix(apply(sampel_exp1, 1, median))

harapan_mean_exp1   = mean(mean_exp1)
harapan_median_exp1 = mean(median_exp1)

# 3. Sebaran Seragam
set.seed(42)

populasi3    = runif(20)
mean_pop3    = mean(populasi3)

sampel_unif1 = urnsamples(populasi3,
                          size = 5,
                          replace = FALSE,
                          ordered = FALSE)

mean_unif1   = matrix(apply(sampel_unif1, 1, mean))
median_unif1 = matrix(apply(sampel_unif1, 1, median))

harapan_mean_unif1   = mean(mean_unif1)
harapan_median_unif1 = mean(median_unif1)

hasil = data.frame(
  "Hasil" = c("mean_populasi",
              "harapan_mean_contoh",
              "harapan_median_contoh"),

  "Sebaran Normal" =
    c(mean_pop1,
      harapan_mean_norm1,
      harapan_median_norm1),

  "Sebaran Eksponensial" =
    c(mean_pop2,
      harapan_mean_exp1,
      harapan_median_exp1),

  "Sebaran Seragam" =
    c(mean_pop3,
      harapan_mean_unif1,
      harapan_median_unif1)
)

hasil
##                   Hasil Sebaran.Normal Sebaran.Eksponensial Sebaran.Seragam
## 1         mean_populasi      0.1919200            0.6965990       0.6131464
## 2   harapan_mean_contoh      0.1919200            0.6965990       0.6131464
## 3 harapan_median_contoh      0.2692767            0.5680749       0.6339751
#Kesimpulan :
#Berdasarkan hasil simulasi pada distribusi normal, eksponensial, dan seragam, diperoleh bahwa nilai harapan rataan sampel (x̄) mendekati bahkan sama dengan mean populasi (μ), sedangkan nilai harapan median sampel masih berbeda terhadap μ. Hal ini menunjukkan bahwa rataan sampel (x̄) merupakan penduga tak bias bagi mean populasi.

#Pada populasi terhingga, seluruh kemungkinan sampel dapat dibentuk sehingga nilai harapan penduga dapat sama persis dengan parameter populasi. Sementara pada populasi tak hingga, sampel yang diperoleh hanya sebagian dari seluruh kemungkinan yang ada sehingga hasil pendugaan tidak selalu sama persis, namun tetap mendekati parameter populasi.
# Penerapan ketakbiasan penduga (Ragam)

# POPULASI TERHINGGA
library(probs)

# Sebaran Normal
set.seed(42)

n        = 5
populasi = rnorm(20)

# Ragam populasi
sigma2   = var(populasi)*(20-1)/20

# Membentuk seluruh kemungkinan sampel
sampel = urnsamples(populasi,
                    size = n,
                    replace = FALSE,
                    ordered = FALSE)

# Ragam contoh dengan pembagi (n-1)
s2.n1   = matrix(apply(sampel, 1, var))
E.s2.n1 = mean(s2.n1)

# Ragam contoh dengan pembagi (n)
s2.n    = s2.n1*(n-1)/n
E.s2.n  = mean(s2.n)

# Sebaran Eksponensial
set.seed(42)

populasi2  = rexp(20)

# Ragam populasi
sigma2.exp = var(populasi2)*(20-1)/20

# Membentuk seluruh kemungkinan sampel
sampel_exp = urnsamples(populasi2,
                        size = n,
                        replace = FALSE,
                        ordered = FALSE)

# Ragam contoh dengan pembagi (n-1)
s2.n1.exp   = matrix(apply(sampel_exp, 1, var))
E.s2.n1.exp = mean(s2.n1.exp)

# Ragam contoh dengan pembagi (n)
s2.n.exp   = s2.n1.exp*(n-1)/n
E.s2.n.exp = mean(s2.n.exp)

hasil = data.frame(
  "." = c("ragam populasi",
          "nilai harapan ragam contoh (n-1)",
          "nilai harapan ragam contoh (n)"),

  "Sebaran Normal" =
    c(sigma2,
      E.s2.n1,
      E.s2.n),

  "Sebaran Eksponensial" =
    c(sigma2.exp,
      E.s2.n1.exp,
      E.s2.n.exp)
)

hasil
##                                  . Sebaran.Normal Sebaran.Eksponensial
## 1                   ragam populasi       1.636844            0.3504577
## 2 nilai harapan ragam contoh (n-1)       1.722993            0.3689028
## 3   nilai harapan ragam contoh (n)       1.378395            0.2951223
#Kesimpulan : Berdasarkan hasil simulasi pada populasi terhingga dengan distribusi normal dan eksponensial, diperoleh bahwa nilai harapan ragam sampel dengan penyebut (n−1) lebih mendekati ragam populasi dibandingkan ragam sampel dengan penyebut n. Hal ini menunjukkan bahwa ragam sampel dengan penyebut (n−1) merupakan penduga tak bias bagi ragam populasi (σ²). Karena populasi bersifat terhingga, seluruh kombinasi sampel dapat dibentuk tanpa pengulangan sehingga hasil pendugaan menjadi lebih stabil dan mendekati parameter populasi. Namun, nilai harapan penduga ragam tidak selalu sama persis dengan ragam populasi karena adanya faktor koreksi dalam perhitungan ragam populasi.
# Selang Kepercayaan Sebaran Normal

# Ukuran sampel pertama
n1     = 15
k      = 100
alpha  = 0.05
mu     = 50
std    = 10

# Mengatur seed random
set.seed(42)

# Membentuk data sampel
sampel.norm1 = matrix(rnorm(n1*k, mu, std), k)

# Menghitung rata-rata sampel
xbar.norm1   = apply(sampel.norm1, 1, mean)

# Menghitung standar deviasi sampel
s.norm1      = apply(sampel.norm1, 1, sd)

# Menghitung standard error
SE.norm1     = s.norm1/sqrt(n1)

# Menentukan nilai z
z.norm1      = qnorm(1-alpha/2)

# Membentuk selang kepercayaan
SK.norm1     = (xbar.norm1-z.norm1*SE.norm1 < mu &
                mu < xbar.norm1+z.norm1*SE.norm1)

# Proporsi SK yang memuat mean populasi
x.norm1      = sum(SK.norm1)/k

# Ukuran sampel kedua
n2     = 40

set.seed(42)

sampel.norm2 = matrix(rnorm(n2*k, mu, std), k)

xbar.norm2   = apply(sampel.norm2, 1, mean)

s.norm2      = apply(sampel.norm2, 1, sd)

SE.norm2     = s.norm2/sqrt(n2)

z.norm2      = qnorm(1-alpha/2)

SK.norm2     = (xbar.norm2-z.norm2*SE.norm2 < mu &
                mu < xbar.norm2+z.norm2*SE.norm2)

x.norm2      = sum(SK.norm2)/k

# Ukuran sampel ketiga
n3     = 80

set.seed(42)

sampel.norm3 = matrix(rnorm(n3*k, mu, std), k)

xbar.norm3   = apply(sampel.norm3, 1, mean)

s.norm3      = apply(sampel.norm3, 1, sd)

SE.norm3     = s.norm3/sqrt(n3)

z.norm3      = qnorm(1-alpha/2)

SK.norm3     = (xbar.norm3-z.norm3*SE.norm3 < mu &
                mu < xbar.norm3+z.norm3*SE.norm3)

x.norm3      = sum(SK.norm3)/k

hasil = data.frame(
  "n" = c(15, 40, 80),
  "Ketepatan SK Sebaran Normal" =
    c(x.norm1, x.norm2, x.norm3)
)

hasil
##    n Ketepatan.SK.Sebaran.Normal
## 1 15                        0.92
## 2 40                        0.95
## 3 80                        0.93
# Membuat plot selang kepercayaan
matplot(
  rbind(xbar.norm2-z.norm2*SE.norm2,
        xbar.norm2+z.norm2*SE.norm2),

  rbind(1:k,1:k),

  col = ifelse(SK.norm2, "skyblue", "red"),

  type = "l",
  lty = 1,

  main = "Selang Kepercayaan 95% (n = 40)",

  xlab = "SK",
  ylab = "Banyak Ulangan"
)

# Garis vertikal mean populasi
abline(v = mu)

# Gambar di atas merupakan hasil simulasi selang kepercayaan 95% untuk rata-rata populasi (μ = 50) dengan ukuran sampel n = 40. Simulasi dilakukan sebanyak k = 100 kali sehingga terbentuk 100 selang kepercayaan.

# Garis vertikal pada x = 50 menunjukkan nilai mean populasi sebenarnya. Setiap garis horizontal merepresentasikan satu selang kepercayaan dari hasil pengambilan sampel. Selang kepercayaan yang memotong garis vertikal berarti berhasil memuat nilai μ, sedangkan yang tidak memotong berarti gagal memuat μ.

# Berdasarkan hasil simulasi, sebagian besar selang kepercayaan berhasil mencakup nilai mean populasi. Hal ini menunjukkan bahwa selang kepercayaan 95% memiliki tingkat ketepatan yang mendekati nilai (1 − α) = 0.95. Semakin besar ukuran sampel, maka estimasi yang diperoleh akan semakin stabil dan proporsi selang kepercayaan yang memuat parameter populasi cenderung semakin mendekati nilai teoritisnya.
# Interval Kepercayaan

# Memanggil package car
library(car)
## Warning: package 'car' was built under R version 4.4.3
## Loading required package: carData
## Warning: package 'carData' was built under R version 4.4.3
# Memanggil dataset Prestige
data("Prestige")

# Menghitung rata-rata income
m <- mean(Prestige$income)
m
## [1] 6797.902
# Menghitung standard error
p <- dim(Prestige)[1]

se <- sd(Prestige$income)/sqrt(p)
se
## [1] 420.4089
# Menentukan nilai kritis t
tval <- qt(0.975, df = p-1)

# Menghitung interval kepercayaan 95%
cat(
  paste(
    "KI 95%: [",
    round(m - tval*se, 2),
    ",",
    round(m + tval*se, 2),
    "]"
  )
)
## KI 95%: [ 5963.92 , 7631.88 ]
#Berdasarkan hasil perhitungan interval kepercayaan 95%, diperoleh bahwa rata-rata pendapatan populasi diperkirakan berada pada rentang interval yang dihasilkan. Artinya, dengan tingkat kepercayaan 95%, interval tersebut diyakini mampu mencakup nilai rata-rata pendapatan populasi yang sebenarnya.