Komputasi Statistika

~ Probability ~

Kontak	: \(\downarrow\)
Email	dsciencelabs@outlook.com
Instagram	https://www.instagram.com/dsciencelabs/
RPubs	https://rpubs.com/dsciencelabs/

Distribusi Binomial

Kasus 1

Berapa peluang munculnya 2 kepala dalam 10 pelemparan koin dimana peluang munculnya kepala adalah 0,3?

dbinom(x = 2, size = 10, prob = 0.3)          # exact

## [1] 0.2334744

mean(rbinom(n=10000,size=10,prob=0.3)==2)     # simulated

## [1] 0.233

Detail

library(dplyr)
library(ggplot2)
library(scales)
data.frame(heads = 0:10, prob = dbinom(x = 0:10, size = 10, prob = 0.3)) %>%
  mutate(Heads = ifelse(heads == 2, "2", "other")) %>%
  ggplot(aes(x = factor(heads), y = prob, fill = Heads)) +
  theme_minimal()+
  geom_col() +
  geom_text(
    aes(label = round(prob,2), y = prob + 0.01),
    position = position_dodge(0.9),
    size = 3,
    vjust = 0) +
  labs(title = "Probability of X = 2 successes",
       subtitle = "b(10, .3)",
       x = "Successes (x)",
       y = "probability")

Kasus 2

Berapa probabilitas <=5 kepala dalam 10 pelemparan koin di mana probabilitas kepala adalah 0,3?

dbinom(0, size=10, prob=0.3) + 
  dbinom(1, size=10, prob=0.3) + 
  dbinom(2, size=10, prob=0.3) + 
  dbinom(3, size=10, prob=0.3) + 
  dbinom(4, size=10, prob=0.3) +
  dbinom(5, size=10, prob=0.3)                # manual calculation for binomial distribution

## [1] 0.952651

pbinom(5, size=10, prob=0.3)                  # alternative way for binomial distribution

## [1] 0.952651

pbinom(q=5,size=10,p=0.3,lower.tail=TRUE)     # alternative way for binomial distribution

## [1] 0.952651

mean(rbinom(n=10000,size=10,prob=0.3)<= 5)    # simulated

## [1] 0.953

Detail

library(dplyr)
library(ggplot2)
data.frame(heads = 0:10, 
           pmf = dbinom(x = 0:10, size = 10, prob = 0.3),
           cdf = pbinom(q = 0:10, size = 10, prob = 0.3, 
                        lower.tail = TRUE)) %>%
  mutate(Heads = ifelse(heads <= 5, "<=5", "other")) %>%
  ggplot(aes(x = factor(heads), y = cdf, fill = Heads)) +
  geom_col() +
  theme_minimal()+
  geom_text(
    aes(label = round(cdf,2), y = cdf + 0.01),
    position = position_dodge(0.9),
    size = 3,
    vjust = 0) +
  labs(title = "Probability of X <= 5 successes",
       subtitle = "b(10, .3)",
       x = "Successes (x)",
       y = "probability")

Kasus 3

Berapa peluang >=5 kepala dalam 10 pelemparan koin dimana peluang munculnya kepala adalah 0,3?

library(dplyr)
library(ggplot2)
data.frame(heads = 0:10, 
           pmf = dbinom(x = 0:10, size = 10, prob = 0.3),
           cdf = pbinom(q = -1:9, size = 10, prob = 0.3, 
                        lower.tail = FALSE)) %>%
  mutate(Heads = ifelse(heads >= 5, ">=5", "other")) %>%
  ggplot(aes(x = factor(heads), y = cdf, fill = Heads)) +
  geom_col() +
  theme_minimal()+
  geom_text(
    aes(label = round(cdf,2), y = cdf + 0.01),
    position = position_dodge(0.9),
    size = 3,
    vjust = 0) +
  labs(title = "Probability of X >= 5 successes",
       subtitle = "b(10, .3)",
       x = "Successes (x)",
       y = "probability")

Kasus 4

Berapa jumlah kepala yang diharapkan dalam 25 pelemparan koin di mana peluang munculnya kepala adalah 0,3?

25 * 0.3                                      # exact expected number of heads in 25 coin flips

## [1] 7.5

mean(rbinom(n = 10000, size = 25, prob = .3)) # exact expected number of heads in 25 coin flips

## [1] 7.5079

25 * 0.3 * (1 - 0.3)                          # variance

## [1] 5.25

var(rbinom(n = 10000, size = 25, prob = .3))  # variance

## [1] 5.196304

Detail

library(dplyr)
library(ggplot2)
data.frame(heads = 0:25, 
           pmf = dbinom(x = 0:25, size = 25, prob = 0.3)) %>%
  mutate(Heads = ifelse(heads == 7, "7", "other"))%>%
  ggplot(aes(x = factor(heads), y = pmf, fill = Heads)) +
  geom_col() +
  theme_minimal()+
  geom_text(
    aes(label = round(pmf,2), y = pmf + 0.01),
    position = position_dodge(0.9),
    size = 3,
    vjust = 0) +
  labs(title = "Probability of X = x successes.",
       subtitle = "b(25, .3)",
       x = "Successes (x)",
       y = "probability")

Latihan 1

Misalkan ada dua puluh pertanyaan pilihan ganda dalam kuis kelas Statistik. Setiap pertanyaan memiliki lima kemungkinan jawaban, dan hanya satu yang benar. Temukan probabilitas memiliki empat atau kurang jawaban yang benar jika seorang siswa mencoba menjawab setiap pertanyaan secara acak. Jawab:
20 pertanyaan dan hanya satu yang benar dari lima kemungkinan jawaban yang benar secara acak = 1/5=0.2
Menemukan probabilitas memiliki empat atau kurang jawaban yang benar.

dbinom(4, size = 20, prob = 0.2)

## [1] 0.2181994

Untuk menemukan probabilitas <=4 yang benar dengan acak, menggunakan fungsi dbinom dengan \(x=0,...,4.\)

dbinom(0, size=20, prob=0.2) + 
+ dbinom(1, size=20, prob=0.2) + 
+ dbinom(2, size=20, prob=0.2) + 
+ dbinom(3, size=20, prob=0.2) + 
+ dbinom(4, size=20, prob=0.2)

## [1] 0.6296483

pbinom(4, size=20, prob=0.2)

## [1] 0.6296483

pbinom(q=4, size=20, p=0.2, lower.tail=TRUE)

## [1] 0.6296483

mean(rbinom(n=1000, size=20, prob=0.2)<=4)

## [1] 0.617

Detail visualisasi

library(dplyr)
library(ggplot2)
data.frame(heads = 0:10, 
           pmf = dbinom(x = 0:10, size = 20, prob = 0.2),
           cdf = pbinom(q = 0:10, size = 20, prob = 0.2, 
                        lower.tail = TRUE)) %>%
  mutate(Heads = ifelse(heads <= 4, "<=4", "other")) %>%
  ggplot(aes(x = factor(heads), y = cdf, fill = Heads)) +
  geom_col() +
  theme_minimal()+
  geom_text(
    aes(label = round(cdf,2), y = cdf + 0.01),
    position = position_dodge(0.9),
    size = 3,
    vjust = 0) +
  labs(title = "Probability of X <= 4 successes",
       subtitle = "b(20, .2)",
       x = "Successes (x)",
       y = "probability")

Distribusi Poisson

Distribusi Poisson adalah distribusi probabilitas diskrit . Ini memberikan probabilitas suatu peristiwa terjadi beberapa kali ( k ) dalam interval waktu atau ruang tertentu.
Distribusi Poisson hanya memiliki satu parameter , (lambda), yang merupakan jumlah rata -rata kejadian. Grafik di bawah ini menunjukkan contoh distribusi Poisson dengan nilai yang berbeda.

Kasus 5

Berapa peluang menghasilkan 2 sampai 4 penjualan dalam seminggu jika tingkat penjualan rata-rata adalah 3 per minggu?

library(ggplot2)
library(dplyr)
# Using cumulative probability
ppois(q = 4, lambda = 3, lower.tail = TRUE) - 
  ppois(q = 1, lambda = 3, lower.tail = TRUE)

## [1] 0.616115

# Using exact probability
dpois(x = 2, lambda = 3) +
  dpois(x = 3, lambda = 3) +
  dpois(x = 4, lambda = 4)

## [1] 0.6434504

# variance = lambda = 3
options(scipen = 999, digits = 2) # sig digits
events <- 0:10
density <- dpois(x = events, lambda = 3)
prob <- ppois(q = events, lambda = 3, lower.tail = TRUE)
df <- data.frame(events, density, prob)
ggplot(df, aes(x = factor(events), y = density)) +
  theme_minimal()+
  geom_col() +
  geom_text(
    aes(label = round(density,2), y = density + 0.01),
    position = position_dodge(0.9),
    size = 3,
    vjust = 0) +
  labs(title = "PMF and CDF of Poisson Distribution",
       subtitle = "P(3).",
       x = "Events (x)",
       y = "Density") +
  geom_line(data = df, aes(x = events, y = prob))

Kasus 6

library(ggplot2)
library(dplyr)
ppois(q=150,lambda=.300*500,lower.tail=TRUE)  # probability of x <= 150

## [1] 0.52

dpois(x=150,lambda=.300*500)                  # probability of x = 150

## [1] 0.033

ppois(q=150,lambda=.300*500,lower.tail=FALSE) # probability of x > 150

## [1] 0.48

options(scipen = 999, digits = 2)             # sig digits
hits <- 0:100 * 3
density <- dpois(x = hits, lambda = .300 * 500)
prob <- ppois(q = hits, lambda = .300 * 500, lower.tail = TRUE)
df <- data.frame(hits, density, prob)
ggplot(df, aes(x = hits, y = density)) +
  geom_col() +
  theme_minimal()+
  labs(title = "Poisson(150)",
       subtitle = "PMF and CDF of Poisson(3) distribution.",
       x = "Hits (x)",
       y = "Density") +
  geom_line(data = df, aes(x = hits, y = prob))

Kasus 7

Apa distribusi keberhasilan dari sampel? \(n=50\) bila peluang berhasil adalah \(p=.03?\)

library(ggplot2)
library(dplyr)
library(tidyr)
options(scipen = 999, digits = 2) # sig digits
n = 0:10
df <- data.frame(events = 0:10, 
                 Poisson = dpois(x = n, lambda = .03 * 50),
                 Binomial = dbinom(x = n, size = 50, p = .03))
df_tidy <- gather(df, key = "Distribution", value = "density", -c(events))
ggplot(df_tidy, aes(x = factor(events), y = density, fill = Distribution)) +
  geom_col(position = "dodge") +
  theme_minimal()+
  labs(title = "Poisson(15) and Binomial(50, .03)",
       subtitle = "Poisson approximates binomial when n >= 20 and p <= .05.",
       x = "Events (x)",
       y = "Density")

Latihan 2

Jika rata-rata dua puluh mobil melintasi jembatan per menit, visualisasikan dan temukan probabilitas memiliki tiga belas atau lebih mobil melintasi jembatan dalam satu menit tertentu.
Jawab:

ppois(q=13, lambda = 20, lower=FALSE) # upper tail

## [1] 0.93

Sehingga jika terdapat 20 mobil yang melintasi jembatan rata rata per menit, kemungkinan memiliki 13 atau lebih mobil yang melintasi jembatan pada menit tertentu adalah 96%.

Latihan 3

Misalkan probabilitas suatu obat menghasilkan efek samping tertentu adalah p = = 0,1% dan n = 1.000 pasien dalam uji klinis menerima obat tersebut. Berapa probabilitas 0 orang mengalami efek samping dengan menggunakan teknik visualisasi?

dpois(x= 0 , lambda= 1000 * .001)

## [1] 0.37

library(ggplot2)
library(dplyr)
options(scipen = 999, digits  =2)

as <- 0:10
density <- dpois(x = as, lambda = 1000 * 0.001)
prob <- ppois(q =as, lambda = 1000 * 0.001, lower.tail = TRUE)
df <- data.frame(as, density, prob)
ggplot(df, aes(x = as, y = density)) +
  geom_col() +
  geom_text(
    aes(label = round(density,2), y = density + 0.01),
    position = position_dodge(0.9),
    size = 3,
    vjust = 0
  )+
  
  labs(title = "Poisson(1)",
       subtitle = "PMF dan CDF of Poisson(1) distribution.",
       x = "Events (x)",
       y = "Density")+
  geom_line (data=df, aes(x=n, y= prob))

Distribusi Seragam Berkelanjutan

Distribusi seragam adalah distribusi probabilitas yang berkelanjutan dan berkaitan dengan peristiwa yang sama -sama mungkin terjadi.

Kasus 8

Pilih seratus angka acak antara minus satu dan lima.
Kami menerapkan fungsi generasi runifdari distribusi seragam untuk menghasilkan sepuluh angka acak antara minus satu dan lima.

rand.unif <- runif(100, min=-3, max=5)        # ten random numbers between minus one and five
hist(rand.unif, col = "cornflowerblue",               # plot the results as a histogram
     freq = FALSE, 
     xlab = 'x', 
     density = 20)

Selanjutnya, kami memplot keduanya, histogram densitas dari atas serta distribusi probabilitas seragam untuk interval [-3,5], dengan menerapkan dunif()fungsi.

a <- -3
b <- 5
hist(rand.unif, 
     freq = FALSE,
     col = "azure4",
     xlab = 'x',  
     ylim = c(0, 0.2),
     xlim = c(-4,6),
     density = 20,
     main = "Uniform distribution for the interval [-3,5]")
curve(dunif(x, min = a, max = b), 
      from = -4, to = 6, 
      n = 100000, 
      col = "darkblue", 
      lwd = 2, 
      add = TRUE)

Apa itu nilai \(x\) sesuai dengan nilai yang membagi distribusi seragam yang diberikan menjadi dua bagian yang sama, atau ditulis lebih formal

qunif(0.5, min=a, max=b)                      # yield value uniform dist into two equal parts

## [1] 1

punif(1, min= a, max = b, lower.tail = FALSE) # measurements yield values `>=1`

## [1] 0.5

Distribusi Eksponensial

Jika X memiliki distribusi eksponensial dengan mean \(\mu\) maka parameter peluruhannya adalah \(m=\frac{1}{\mu}\), kita tulis \(X\) ~ \(Exp(m)\) where \(x>=0\) dan \(m>0\). Fungsi Kepadatan dari X adalah \(f(x)=\frac{1}{\mu}e^\frac{-x}{\mu}\). fungsi distribusi kumulatif X adalah \(P(X<=x)=1-e^{-mx}\)
Mari kita lihat bagaimana distribusi eksponensial terlihat. Sampel acak diambil dari distribusi eksponensial yang dihasilkan rexp(n)dengan menggunakan ukuran sampel (n).

set.seed(1)                                   # to generate the same random number results
rexp(10)                                      # random number with 10 sample size

##  [1] 0.76 1.18 0.15 0.14 0.44 2.89 1.23 0.54 0.96 0.15

par(mfrow=c(2,2))                             # partition the graph by 2x2 rows and column 
hist(rexp(10),col = "azure1")                 # histogram plot with 10 sample size 
hist(rexp(100),col = "azure2")                # histogram plot with 100 sample size
hist(rexp(1000),col = "azure3")               # histogram plot with 1000 sample size
hist(rexp(10000),col = "azure4")              # histogram plot with 10000 sample size

Kasus 9

Mari, Kami ingin menghasilkan distribusi rata-rata dari 40 eksponensial acak dengan 10.000 simulasi dengan lambda (parameter laju) = 0,2.

library(ggplot2)
set.seed(1)                                   # to generate the same random number results
mns=NULL
for (i in 1 : 10000) mns = c(mns, mean(rexp(40,0.2)))
Exp <- data.frame(mns,size=40)
ggplot(Exp,aes(x=mns,fill=size))+
  geom_histogram(aes(y=..density..),binwidth=.25,col="black") + 
  ylim(c(0,0.6))+
  stat_function(fun=dnorm,arg=list(mean=5,sd=sd(mns)))+
  theme_minimal()+
  geom_vline(aes(xintercept=mean(mns),colour="red"))+
  geom_text(aes(x=mean(mns),
                label="\n sample mean",
                y=0.2),colour="black",
            angle=90, text = element_text(size=11))+
  xlab("Averages of the distribution") + ylab("Frequency")+
  ggtitle("Distribution of the averages of \n 40 random exponentials (10000 simulations)")

Latihan 4

Misalkan waktu checkout rata-rata kasir supermarket adalah tiga menit. Temukan probabilitas checkout pelanggan diselesaikan oleh kasir dalam waktu kurang dari dua menit.

pexp(2, rate=1/3)

## [1] 0.49

Distribusi Normal

Distribusi normal adalah istilah yang tepat untuk kurva lonceng probabilitas. Dalam distribusi normal, rata-ratanya adalah nol dan simpangan bakunya adalah 1. Memiliki kemiringan nol dan kurtosis 3. Distribusi normal adalah simetris, tetapi tidak semua distribusi simetris adalah normal.
Mari kita ambil sampel dari distribusi normal dengan \(\bar{Y}=100\) dan \(\sigma=5\):

mydata <- rnorm(n=10000, mean=100, sd=5)
mean(mydata)

## [1] 100

sd(mydata)

## [1] 5

# Make a histogram of your sample
hist(mydata, main ="", 
     col = "blue", 
     xlim = c(80, 120), 
     freq = FALSE,
     xlab = "")
# ...and add a density curve
curve(dnorm(x, 
            mean=mean(mydata), 
            sd=sd(mydata)), add=TRUE, 
            col="black", lwd=3)

Kasus 10

Asumsikan bahwa nilai ujian dari ujian masuk perguruan tinggi sesuai dengan distribusi normal. Selanjutnya, nilai rata-rata tes adalah 75, dan standar deviasi adalah 14,2. Berapa persentase siswa yang mendapat nilai 86 atau lebih dalam ujian?
Kami menerapkan fungsi pnorm distribusi normal dengan mean 75 dan standar deviasi 14.2. Karena kami mencari persentase siswa yang mendapat nilai lebih tinggi dari 86, kami tertarik pada bagian atas dari distribusi normal.

pnorm(86, mean=75, sd=14.2, lower.tail=FALSE)

## [1] 0.22

Persentase siswa yang mendapat nilai 84 atau lebih dalam ujian masuk perguruan tinggi adalah 21,9%

Distribusi Chi-Kuadrat

Distribusi chi-kuadrat (Χ 2 ) adalah keluarga distribusi probabilitas kontinu. Mereka banyak digunakan dalam uji hipotesis , termasuk uji kecocokan chi-kuadrat dan uji independensi chi-kuadrat.
Bentuk distribusi chi-kuadrat ditentukan oleh parameter k , yang mewakili derajat kebebasan. Ketika derajat kebebasan meningkat, distribusi chi-kuadrat berubah dari kurva ke bawah menjadi bentuk punuk. Saat derajat kebebasan meningkat lebih jauh, punuk berubah dari sangat miring ke kanan menjadi kira-kira normal.

Kasus 11

Temukan persentil ke-95 dari distribusi Chi-Squared dengan 6 derajat kebebasan.
Kami menerapkan fungsi kuantil qchisqdari distribusi Chi-Squared terhadap nilai desimal 0,95.

qchisq(.95, df=6)                                      # 6 degrees of freedom

## [1] 13

Persentil ke-95 dari distribusi Chi-Squared dengan 6 derajat kebebasan adalah 12.59159.

Distribusi T-Student

Distribusi t-student adalah distribusi probabilitas kontinu.

Kasus 12

qt(c(.025, .975), df=5)                                # 5 degrees of freedom

## [1] -2.6  2.6

Distribusi 1.8 F

Jika V1 dan V2 adalah dua variabel acak independen yang memiliki distribusi Chi-Squared dengan m1 dan m2 derajat kebebasan masing-masing, maka kuantitas berikut mengikuti an F distribusi dengan m1 derajat kebebasan pembilang dan m2 derajat kebebasan penyebut yaitu (m1, m2) derajat kebebasan.

Kasus 13

Temukan persentil ke-95 dari distribusi F dengan (5, 2) derajat kebebasan. Kami menerapkan fungsi kuantil qf dari F distribusi terhadap nilai desimal 0,95.

qf(.95, df1=5, df2=2)

## [1] 19