#   Buktikan dengan pendekatan simulasi jika σ^2 diketahui maka rata-rata menyebar N(μ,σ^2/n)!
set.seed(123)
M = 10000       # Jumlah ulangan simulasi
n = 30          # Ukuran sampel
mu = 50         # Rata-rata populasi
sigma = 10      # Simpangan baku populasi (diketahui)

# Vektor untuk menyimpan hasil nilai Z
Z_values = numeric(M)

for (i in 1:M) {
  # Bangkitkan sampel acak
  sampel = rnorm(n, mean = mu, sd = sigma)
  
  # Hitung rata-rata sampel
  x_bar = mean(sampel)
  
  # Hitung nilai Z (karena sigma diketahui)
  Z_values[i] = (x_bar - mu) / (sigma / sqrt(n))
}

# Membuat Plot
hist(Z_values, prob = TRUE, breaks = 50, col = "lightblue", 
     main = "Pembuktian Sebaran Normal (Sigma Diketahui)", 
     xlab = "Nilai Z")
# Overlay kurva normal baku
curve(dnorm(x, mean = 0, sd = 1), col = "red", lwd = 2, add = TRUE)

# Buktikan dengan pendekatan simulai jika σ^2 tidak diketahui maka rata-rata menyebar t-student dengan derajat bebas n-1!
set.seed(123)
M = 10000
n = 5 
mu = 50
t_values = numeric(M)

for (i in 1:M) {
  sampel = rnorm(n, mean = mu, sd = 10) 
  x_bar = mean(sampel)
  s = sd(sampel) 
  t_values[i] = (x_bar - mu) / (s / sqrt(n))
}

hist(t_values, prob = TRUE, breaks = 50, xlim = c(-5,5), main = "Sebaran t-student")
curve(dt(x, df = n-1), col = "blue", lwd = 2, add = TRUE)

# Buktikan dengan pendekatan simulasi, Berdasarkan dalil limit pusat, walau σ2 tidak diketahui asalkan ukuran sampel besar (n>30) maka sebaran dari rata-rata dapat juga diaproksimasi dengan sebaran N(μ,s2/n)!
set.seed(123)
M = 10000
n = 500 # Sampel Besar
mu = 50
t_clt = numeric(M)

for (i in 1:M) {
  sampel = rnorm(n, mean = mu, sd = 10)
  x_bar = mean(sampel)
  s = sd(sampel)  
  t_clt[i] = (x_bar - mu) / (s / sqrt(n))
}

hist(t_clt, prob = TRUE, breaks = 50, xlim = c(-4,4), main = "Aproksimasi Normal (CLT)")
curve(dnorm(x, mean = 0, sd = 1), col = "red", lwd = 2, add = TRUE)

# Buktikan dengan pendekatan simulasi bahwa nilai harapan rata-rata contoh adalah μ!
set.seed(123)
M = 10000
n = 30
mu_pop = 50
x_bar_values = numeric(M)

for (i in 1:M) {
  sampel = rnorm(n, mean = mu_pop, sd = 10)
  x_bar_values[i] = mean(sampel)
}

Nilai_Harapan = mean(x_bar_values)
print(Nilai_Harapan)
## [1] 50.01865
# Tunjukkan dengan pendekatan  simulasi, apakah s merupakan penduga tak bias dari σ!
set.seed(123)
M = 10000
n = 10
sigma_kuadrat = 100 

var_tak_bias = numeric(M)
var_bias = numeric(M)

for (i in 1:M) {
  sampel = rnorm(n, mean = 50, sd = sqrt(sigma_kuadrat))
  var_tak_bias[i] = var(sampel) 
  var_bias[i] = var(sampel) * (n - 1) / n
}

print(mean(var_tak_bias)) # Akan mendekati 100
## [1] 100.1031
print(mean(var_bias))     # Akan bias ke bawah
## [1] 90.09275
# Tunjukan dengan menggunakan pendekatan simulasi, jika sampel diambil dari populasi normal,  maka sebaran dari rata-rata jika σ^2 tidak diketahui adalah N(μ,S^2/n)!
set.seed(123)
M = 10000
n = 8 
mu_norm = 0 
t_norm_pop = numeric(M)

for (i in 1:M) {
  # Bangkitkan dari distribusi Normal Baku
  sampel = rnorm(n, mean = mu_norm, sd = 1) 
  x_bar = mean(sampel)
  s = sd(sampel)
  t_norm_pop[i] = (x_bar - mu_norm) / (s / sqrt(n))
}

hist(t_norm_pop, prob = TRUE, breaks = 50, xlim = c(-4,4), main = "Dari Populasi Normal")
curve(dt(x, df = n-1), col = "red", lwd = 2, add = TRUE)