library(probs)
## 
## Attaching package: 'probs'
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, union
par(mfrow=c(3,1))
set.seed(123)
populasi = rgeom(20,0.1)
n1 = 2
contoh_geo1 = urnsamples(populasi, size = 2, replace = F, ordered = F)
mean_geo1 = matrix(apply(contoh_geo1, 1, mean))
hist(mean_geo1,main = paste("Hampiran Normal Terhadap Geometrik (n = 2)"),xlab = "xbar")

n2 = 5
contoh_geo2 = urnsamples(populasi, size = 5, replace = F, ordered = F)
mean_geo2 = matrix(apply(contoh_geo2, 1, mean))
hist(mean_geo2,main = paste("Hampiran Normal Terhadap Geometrik (n = 5)"),xlab = "xbar")

n3 = 10
contoh_geo3 = urnsamples(populasi, size = 10, replace = F, ordered = F)
mean_geo3 = matrix(apply(contoh_geo3, 1, mean))
hist(mean_geo3,main = paste("Hampiran Normal Terhadap Geometrik (n = 10)"),xlab = "xbar")

library(probs)
set.seed(123)
par(mfrow=c(3,1))
populasi = rexp(20)
n1 = 2
contoh_exp1 = urnsamples(populasi, size = 2, replace = F, ordered = F)
mean_exp1 = matrix(apply(contoh_exp1, 1, mean))
hist(mean_exp1,main = paste("Hampiran Normal Terhadap Eksponensial (n = 2)"),xlab = "xbar")

n2 = 5
contoh_exp2 = urnsamples(populasi, size = 5, replace = F, ordered = F)
mean_exp2 = matrix(apply(contoh_exp2, 1, mean))
hist(mean_exp2,main = paste("Hampiran Normal Terhadap Eksponensial (n = 5)"),xlab = "xbar")

n3 = 10
contoh_exp3 = urnsamples(populasi, size = 10, replace = F, ordered = F)
mean_exp3 = matrix(apply(contoh_exp3, 1, mean))
hist(mean_exp3,main = paste("Hampiran Normal Terhadap Eksponensial (n = 10)"),xlab = "xbar")

library(probs)
set.seed(123)
par(mfrow=c(3,1))
populasi = runif(20)
n1 = 2
contoh_unif1 = urnsamples(populasi, size = 2, replace = F, ordered = F)
mean_unif1   = matrix(apply(contoh_unif1, 1, mean))
hist(mean_unif1,main = paste("Hampiran Normal Terhadap Seragam (n = 2)"),xlab = "xbar")

n2 = 5
contoh_unif2 = urnsamples(populasi, size = 5, replace = F, ordered = F)
mean_unif2   = matrix(apply(contoh_unif2, 1, mean))
hist(mean_unif2,main = paste("Hampiran Normal Terhadap Seragam (n = 5)"),xlab = "xbar")

n3 = 10
contoh_unif3 = urnsamples(populasi, size = 10, replace = F, ordered = F)
mean_unif3   = matrix(apply(contoh_unif3, 1, mean))
hist(mean_unif3,main = paste("Hampiran Normal Terhadap Seragam (n = 10)"),xlab = "xbar")

library(probs)
set.seed(123)
par(mfrow=c(3,1))
populasi = rnorm(20,5,sqrt(12)) # Membangkitkan bil. acak ~ Normal (miu = 5, sigma2 =12) 
n1 = 3
contoh_norm1 = urnsamples(populasi, size = 3, replace = F, ordered = F)
mean_norm1 = matrix(apply(contoh_norm1, 1, mean))
mean_xbar1 = mean(mean_norm1)
var_xbar1 = var(mean_norm1)
hist(mean_norm1,main = paste("(n = 3)"),xlab = "xbar")

n2 = 4
contoh_norm2 = urnsamples(populasi, size = 4, replace = F, ordered = F)
mean_norm2 = matrix(apply(contoh_norm2, 1, mean))
mean_xbar2 = mean(mean_norm2)
var_xbar2 = var(mean_norm2)
hist(mean_norm2,main = paste("(n = 4)"),xlab = "xbar")

n3 = 15
contoh_norm3 = urnsamples(populasi, size = 15, replace = F, ordered = F)
mean_norm3 = matrix(apply(contoh_norm3, 1, mean))
mean_xbar3 = mean(mean_norm3)
var_xbar3 = var(mean_norm3)
hist(mean_norm3,main = paste("(n = 15)"),xlab = "xbar")

hasil = data.frame("."=c("mean","varian"),"Populasi"=c(5,12),"n=3"=c(mean_xbar1,var_xbar1),"n=4"=c(mean_xbar2,var_xbar2),"n=15"=c(mean_xbar3,var_xbar3))
hasil
##        . Populasi      n.3      n.4      n.15
## 1   mean        5 5.490599 5.490599 5.4905992
## 2 varian       12 3.219489 2.271056 0.1892278
#POPULASI TERHINGGA
#1. Sebaran Normal 
library(probs)
set.seed(123)
n = 10
populasi1 = rnorm(20)     
mean_pop1 = mean(populasi1)
sampel_normal1 = urnsamples(populasi1, size = 10, replace = F, ordered = F)
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
library(probs)
set.seed(123)
n = 10
populasi2 = rexp(20)     
mean_pop2 = mean(populasi2)    
sampel_exp1 = urnsamples(populasi2, size = 10, replace = F, ordered = F)  
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. Uniform
library(probs)
set.seed(123)
n = 10              
populasi3 = runif(20)      
mean_pop3 = mean(populasi3)     
sampel_unif1 = urnsamples(populasi3, size = 10, replace = F, ordered = F)
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.1416238            0.8111726       0.5508084
## 2   harapan_mean_contoh      0.1416238            0.8111726       0.5508084
## 3 harapan_median_contoh      0.1174878            0.4931612       0.5504018
# POPULASI TERHINGGA
#Sebaran Normal
library(probs)
set.seed(123)
n = 10
populasi = rnorm(20) 
sigma2 = var(populasi)*(20-1)/20 #fungsi var pada R adalah varian contoh (penyebut n-1) sehingga perlu dikali (n-1)/n
sampel = urnsamples(populasi, size = 10, replace = F, ordered = F)
## Pembagi (n-1)
s2.n1 = matrix(apply(sampel, 1, var))       
E.s2.n1 = mean(s2.n1)
## Pembagi (n)
s2.n = s2.n1*(10-1)/10       
E.s2.n = mean(s2.n)
#Sebaran Eksponensial
library(probs)
set.seed(123)
n = 10
populasi2 = rexp(20) 
sigma2.exp = var(populasi2)*(20-1)/20
sampel_exp = urnsamples(populasi2, size = 10, replace = F, ordered = F)
## Pembagi (n-1)
s2.n1.exp = matrix(apply(sampel_exp, 1, var))
E.s2.n1.exp = mean(s2.n1.exp)
## Pembagi (n)
s2.n.exp = s2.n1.exp*(10-1)/10
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      0.8987739            0.9439256
## 2 nilai harapan ragam contoh (n-1)      0.9460778            0.9936059
## 3   nilai harapan ragam contoh (n)      0.8514700            0.8942453
set.seed(123)
n1 = 10     
k = 100 #ulangan      
alpha  = 0.05
mu = 50     
std = 10
sampel.norm1 = matrix(rnorm(n1*k,mu,std),k)
xbar.norm1 = apply(sampel.norm1,1,mean)
s.norm1 = apply(sampel.norm1,1,sd)      
SE.norm1 = s.norm1/sqrt(n1)   
z.norm1 = qnorm(1-alpha/2)     
SK.norm1 = (xbar.norm1-z.norm1*SE.norm1 < mu & mu < xbar.norm1+z.norm1*SE.norm1)    
x.norm1 = sum(SK.norm1)/k #proporsi banyaknya SK yang memuat mu     
set.seed(123)
n2 = 30     
k = 100 #ulangan     
alpha = 0.05
mu = 50
std = 10
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 #proporsi banyaknya SK yang memuat mu
set.seed(123)
n3 = 100     
k = 100 #ulangan     
alpha  = 0.05
mu = 50
std = 10
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 = s.norm3/sqrt(n3)      
SK.norm3 = (xbar.norm3-z.norm3*SE.norm3 < mu & mu < xbar.norm3+z.norm3*SE.norm3)     
x.norm3 = sum(SK.norm3)/k #proporsi banyaknya SK yang memuat mu
hasil = data.frame("n" =c(10,30,100),"Ketepatan SK Sebaran Normal"=c(x.norm1, x.norm2, x.norm3))
hasil
##     n Ketepatan.SK.Sebaran.Normal
## 1  10                        0.93
## 2  30                        0.93
## 3 100                        0.75
matplot(rbind (xbar.norm2-z.norm2*SE.norm2, xbar.norm2+z.norm2*SE.norm2), rbind(1:k,1:k), col=ifelse(SK.norm2,"blue","red"), type = "l", lty = 1,main='Selang Kepercayaan 95% (n=100)', xlab='SK', ylab='banyak ulangan')
abline(v=mu)

# Interval Kepercayaan
library(car)
## Loading required package: carData
# Menghitung rata-rata
data("Prestige")
m <- mean(Prestige$income)
m
## [1] 6797.902
# Menghitung standar error
p <- dim(Prestige)[1]
se <- sd(Prestige$income)/sqrt(p)
se
## [1] 420.4089
# Menghitung nilai kritis t
tval <- qt(0.975, df=p-1)
# Menghitung interval kepercayaan
cat(paste("KI: [", round(m-tval*se, 2),",",round(m+tval*se,2),"]"))
## KI: [ 5963.92 , 7631.88 ]