R을 이용한 통계 프로그래밍 기초 제4장 난수 발생과 모의실험

프로그램 소스

# R을 이용한 통계 프로그래밍 기초

# 4 난수 발생과 모의실험

# 4.1 연속분포

# 4.1.1 정규분포

rnorm(1, 100, 16) # 평균이 100이고 표준편차가 16인 정규분포에서 1개 난수 발생
rnorm(5, mean = 280, sd = 10) # 평균이 280이고 표준편차가 10인 정규분포에서 5개 난수 발생

pnorm(265.5393, mean = 280, sd = 10) # 평균이 280이고 표준편차가 10인 분포에서 X <= 295.5393의 확률

qnorm(0.05, mean = 0, sd = 1) # 평균이 0이고, 표준편차가 1인 분포에서 누적확률이 0.05인 x
qnorm(0.95, mean = 0, sd = 1) # 평균이 0이고, 표준편차가 1인 분포에서 누적확률이 0.95인 x
qnorm(0.975, mean = 0, sd = 1) # 평균이 0이고, 표준편차가 1인 분포에서 누적확률이 0.975인 x

dnorm(0, mean = 0, sd = 1) # 평균이 0이고, 표준편차가 1인 분포에서 x = 0의 확률

x = rnorm(100) # 평균이 0이고 표준편차가 1인 분포에서 난수 100개 뽑기
hist(x, probability = TRUE, col = gray(.9), main = "normal mu = 0, sigma = 1")
curve(dnorm(x), add = T)

# 4.1.2 t 분포

dt(2, df = 10) # 자유도 10인 t 분포에서 x = 3에서의 확률
pt(2, df = 10) # 자유도 10인 t 분포에서 x <= 3인 확률
qt(0.99, df = 10) # 자유도 10인 t 분포에서 누적확률분포가 99%인 x
rt(10, df = 10) # 자유도 10인 t 분포에서 난수 10개 발생

# 4.1.3 균일분포

dunif(0.5, min = 0, max = 2) # 최소값 0, 최대값 2인 균일분포에서 x = 0.5의 확률
punif(0.5, min = 0, max = 2) # 최소값 0, 최대값 2인 균일분포에서 x <= 0.5의 누적확률분포
punif(0.5, min = 0, max = 2, lower.tail = TRUE) # 최소값 0, 최대값 2인 균일분포에서 x <= 0.5의 누적확률분포, 만일 lower.tail = FALSE이면 x > 0.5의 누적확률분포
qunif(0.5, min = 0, max = 2) # 최소값 0, 최대값 2인 균일분포에서 누적확률이 0.5인 x값
runif(10, min = 0, max = 2) # 최소값 0, 최대값 2인 균일분포에서 10개의 난수 발생

# 4.1.4 지수분포

dexp(3, rate = 2) # 평균이 1/2인 지수분포에서 x=3의 확률
pexp(3, rate = 2) # 평균이 1/2인 지수분포에서 x<=3의 누적확률
qexp(0.5, rate = 2) # 평균이 1/2인 지수분포에서 누적확률이 0.5인 x
rexp(10, rate = 2) # 평균이 1/2인 지수분포에서 10개의 난수 발생

op = par(mfrow=c(2,2)) # 행 2, 열 2로 분할

x = rt(1000, df = 5) # 자유도 5인 t 분포에서 1000개의 난수 발생
y = dt(x, df = 5) # 그에 해당하는 확률밀도
plot(x, y, sub = "t-dist") # 그래프

x = rnorm(1000, 0, 1) # 평균이 0이고 표준편차 1인 정규분포에서 1000개 난수 발생
y = dnorm(x, 0, 1) # 그에 해당하는 확률밀도
plot(x, y, sub = "표준정규분포") # 그래프

x = rexp(1000, rate = 1) # 평균이 1인 지수분포에서 1000개 난수 발생
y = dexp(x, rate = 1)
plot(x, y, sub = "지수분포")

x = rpois(1000, lambda = 3) # 람다가 3인 포아송분포에서 1000개 난수 발생
y = dpois(x, lambda = 3)
plot(x, y, sub = "포아송분포")

par(op)

# 4.2 이산분포

# 4.2.1 베르누이 분포와 이항분포

dbinom(0:5, size = 5, p = 0.5) # 시행횟수 5, 성공확률 0.5인 이항분포에서 각 x 의 확률
pbinom(0:5, size = 5, p = 0.5) # 시행횟수 5, 성공확률 0.5인 이항분포에서 각 x 의 누적확률
qbinom(0.5, size = 5, p = 0.5) # 시행횟수 5, 성공확률 0.5인 이항분포에서 누적확률 0.5인 x
rbinom(3, size = 5, p = 0.5) # 시행횟수 5, 성공확률 0.5인 이항분포에서 3개의 난수 발생

# 4.2.2 포아송 분포

dpois(3, lambda = 3) # 람다가 3인 포아송분포에서 x = 3의 확률
ppois(3, lambda = 3) # 람다가 3인 포아송분포에서 x <= 3의 누적확률
qpois(0.6, lambda = 3) # 람다가 3인 포아송분포에서 누적확률이 0.6이 되는 x
rpois(10, lambda = 3) # 람다가 3인 포아송분포에서 10개의 난수 발생

plot(table(rpois(100, 3)) # 람다가 3인 포아송분포에서 100개의 난수 발생하여 빈도분석
, type = "h" # 그래프 형태는 수직선
, col = "red" # 색깔은 빨강
, lwd = 10, # 선두께는 10
, main = "rpois(100, lambda = 3)") # 제목

# 4.3 정규성 검정

# 4.3.1 히스토그램을 이용한 정규성 검토

op = par(mfrow=c(2,3)) # 2 by 3개의 그래프를 그릴 것

x = seq(6.5, 13.5, by = 0.01) # 정규분포를 그리기 위한 기초자료

for (i in 1:6) { # 6번 반복
y <- rnorm(1000, mean = 10, sd = 1) # 평균이 10이고, 표준편차가 1인 정규분포에서 1000개 난수 발생
hist(y, prob = TRUE, xlim = c(6.5, 13.5), ylim = c(0, 0.5)) # 히스토그램, 확률, x축 한계, y축 한계
lines(x, dnorm(x, mean = 10, sd = 1)) # 평균이 10이고 표준편차가 1인 정규분포 그래프 추가
}

par(op) # 2 by 3 그래프 마지막

# 4.3.2 정규확률그림

x = rnorm(20) # 평균이 0이고 표준편차가 1인 정규분포에서 20 난수발생
hist(x) # 히스토그램
qqnorm(x) # Q-Q 플롯
qqline(x) # Q-Q 라인

# Shapiro 검정 - 정규성 검정

x = rnorm(10) # 10개의 난수 발생
shapiro.test(x) # Shapiro 검정

# 여러 분포에서 발생한 난수를 이용하여 정규성 검정

op = par(mfrow=c(2,2)) # 2 by 2 분할 그래프

x = rnorm(100) # 정규분포에서 난수 발생
qqnorm(x, sub = "Normal") # Q-Q Plot
qqline(x) # Q-Q Line

x = runif(100, 0, 1) # 균일분포에서 난수 발생
qqnorm(x, sub = "Uniform") # Q-Q Plot
qqline(x) # Q-Q Line

x = rexp(100, 1) # 지수분포에서 난수 발생
qqnorm(x, sub = "Exponential") # Q-Q Plot
qqline(x) # Q-Q Line

x = rpois(100, 3) # 포아송분포에서 난수 발생
qqnorm(x, sub = "Poisson") # Q-Q Plot
qqline(x) # Q-Q Line

par(op) # 2 by 2 분할 그래프 끝

# 4.4 모의 실험

# 4.4.1 이항분포의 정규분포 근사

m = 100
p = 0.5

try.n = c(5, 10, 15, 30) # 시행 횟수

op = par(mfrow=c(2,2)) # 2 by 2 분할 그래프 시작

for(i in 1:4){ # 네 번 반복
n = try.n[i] # 시행 횟수 조절
res = rbinom(m, n, p) # 이항분포에서 m개의 난수발생
hist(res, prob = TRUE) # 히스토그램, 확률
curve(dnorm(x, n*p, sqrt(n*p*(1-p))), add = TRUE) # 곡선추가
}

par(op) # 2 by 2 분할 그래프 끝

# 4.4.2 중심극한정리(central limit theorem)와 모의실험(simulation)

m = 50
mx = rep(0, m) # 0이 50개 있는 벡터 생성
n.value = c(5, 10, 15, 30, 50) # 횟수

# 균일분포 난수들의 표본평균들의 정규분포화

plot(0, 0
, type = "n" # 그래프 형태 nothing, 즉 빈 그래프 생성
, xlim = c(0, 1) # x 축 값의 한계
, ylim = c(0, 10) # y 축 값의 한계
, ylab = "density" # y 축 라벨
, xlab = "mx" # x 축 라벨
, main = "uniform mean to normal") # 주 제목

for(k in 1:length(n.value)){ # 5회 반복

n = n.value[k]

for(i in 1:m){ # 50회 반복
mx[i] = mean(runif(n, 0, 1)) # 균일분포에서 n개의 난수발생, 평균계산, 저장
}

lines(density(mx) # 라인추가
, lty = k # 라인의 형태
, col = k) # 색깔의 종류

}

# 지수분포 난수들의 표본평균들의 정규분포화

plot(0, 0
, type = "n" # 그래프 형태 nothing, 즉 빈 그래프 생성
, xlim = c(0, 3) # x 축 값의 한계
, ylim = c(0, 3.7) # y 축 값의 한계
, ylab = "density" # y 축 라벨
, xlab = "mx" # x 축 라벨
, main = "exponential to normal") # 주 제목

for(k in 1:length(n.value)){ # 5회 반복

n = n.value[k]

for(i in 1:m){ # 50회 반복
mx[i] = mean(rexp(n, rate = 1)) # 지수분포에서 n개의 난수발생, 평균계산, 저장
}

lines(density(mx) # 라인추가
, lty = k # 라인의 형태
, col = k) # 색깔의 종류

}

프로그램 결과

> # R을 이용한 통계 프로그래밍 기초
>
> # 4 난수 발생과 모의실험
>
> # 4.1 연속분포
>
> # 4.1.1 정규분포
>
> rnorm(1, 100, 16) # 평균이 100이고 표준편차가 16인 정규분포에서 1개 난수 발생
[1] 115.8602
> rnorm(5, mean = 280, sd = 10) # 평균이 280이고 표준편차가 10인 정규분포에서 5개 난수 발생
[1] 274.6592 292.4493 278.5038 268.5236 281.6879
>
> pnorm(265.5393, mean = 280, sd = 10) # 평균이 280이고 표준편차가 10인 분포에서 X <= 295.5393의 확률
[1] 0.07407878
>
> qnorm(0.05, mean = 0, sd = 1) # 평균이 0이고, 표준편차가 1인 분포에서 누적확률이 0.05인 x
[1] -1.644854
> qnorm(0.95, mean = 0, sd = 1) # 평균이 0이고, 표준편차가 1인 분포에서 누적확률이 0.95인 x
[1] 1.644854
> qnorm(0.975, mean = 0, sd = 1) # 평균이 0이고, 표준편차가 1인 분포에서 누적확률이 0.975인 x
[1] 1.959964
>
> dnorm(0, mean = 0, sd = 1) # 평균이 0이고, 표준편차가 1인 분포에서 x = 0의 확률
[1] 0.3989423
>
> x = rnorm(100) # 평균이 0이고 표준편차가 1인 분포에서 난수 100개 뽑기
> hist(x, probability = TRUE, col = gray(.9), main = "normal mu = 0, sigma = 1")
> curve(dnorm(x), add = T)
>
>
> # 4.1.2 t 분포
>
> dt(2, df = 10) # 자유도 10인 t 분포에서 x = 3에서의 확률
[1] 0.06114577
> pt(2, df = 10) # 자유도 10인 t 분포에서 x <= 3인 확률
[1] 0.963306
> qt(0.99, df = 10) # 자유도 10인 t 분포에서 누적확률분포가 99%인 x
[1] 2.763769
> rt(10, df = 10) # 자유도 10인 t 분포에서 난수 10개 발생
[1] -1.4084871 -0.8089772 -0.5935148 -2.6261923 -0.9957817 0.2826226
[7] 0.3196900 0.7430285 -0.7519106 0.8259090
>
> # 4.1.3 균일분포
>
> dunif(0.5, min = 0, max = 2) # 최소값 0, 최대값 2인 균일분포에서 x = 0.5의 확률
[1] 0.5
> punif(0.5, min = 0, max = 2) # 최소값 0, 최대값 2인 균일분포에서 x <= 0.5의 누적확률분포
[1] 0.25
> punif(0.5, min = 0, max = 2, lower.tail = TRUE) # 최소값 0, 최대값 2인 균일분포에서 x <= 0.5의 누적확률분포, 만일 lower.tail = FALSE이면 x > 0.5의 누적확률분포
[1] 0.25
> qunif(0.5, min = 0, max = 2) # 최소값 0, 최대값 2인 균일분포에서 누적확률이 0.5인 x값
[1] 1
> runif(10, min = 0, max = 2) # 최소값 0, 최대값 2인 균일분포에서 10개의 난수 발생
[1] 1.4216475 0.9142103 0.5447351 0.9028850 0.6442939 1.9250584 1.9905732
[8] 0.8617541 1.8252440 0.4394071
>
> # 4.1.4 지수분포
>
> dexp(3, rate = 2) # 평균이 1/2인 지수분포에서 x=3의 확률
[1] 0.004957504
> pexp(3, rate = 2) # 평균이 1/2인 지수분포에서 x<=3의 누적확률
[1] 0.9975212
> qexp(0.5, rate = 2) # 평균이 1/2인 지수분포에서 누적확률이 0.5인 x
[1] 0.3465736
> rexp(10, rate = 2) # 평균이 1/2인 지수분포에서 10개의 난수 발생
[1] 0.07753687 0.01760280 0.47249418 0.17080055 0.55388220 0.08268597
[7] 0.13543911 0.31719901 0.02750894 0.01618359
>
>
> op = par(mfrow=c(2,2)) # 행 2, 열 2로 분할
>
> x = rt(1000, df = 5) # 자유도 5인 t 분포에서 1000개의 난수 발생
> y = dt(x, df = 5) # 그에 해당하는 확률밀도
> plot(x, y, sub = "t-dist") # 그래프
>
> x = rnorm(1000, 0, 1) # 평균이 0이고 표준편차 1인 정규분포에서 1000개 난수 발생
> y = dnorm(x, 0, 1) # 그에 해당하는 확률밀도
> plot(x, y, sub = "표준정규분포") # 그래프
>
> x = rexp(1000, rate = 1) # 평균이 1인 지수분포에서 1000개 난수 발생
> y = dexp(x, rate = 1)
> plot(x, y, sub = "지수분포")
>
> x = rpois(1000, lambda = 3) # 람다가 3인 포아송분포에서 1000개 난수 발생
> y = dpois(x, lambda = 3)
> plot(x, y, sub = "포아송분포")
>
> par(op)
>
> # 4.2 이산분포
>
> # 4.2.1 베르누이 분포와 이항분포
>
> dbinom(0:5, size = 5, p = 0.5) # 시행횟수 5, 성공확률 0.5인 이항분포에서 각 x 의 확률
[1] 0.03125 0.15625 0.31250 0.31250 0.15625 0.03125
> pbinom(0:5, size = 5, p = 0.5) # 시행횟수 5, 성공확률 0.5인 이항분포에서 각 x 의 누적확률
[1] 0.03125 0.18750 0.50000 0.81250 0.96875 1.00000
> qbinom(0.5, size = 5, p = 0.5) # 시행횟수 5, 성공확률 0.5인 이항분포에서 누적확률 0.5인 x
[1] 2
> rbinom(3, size = 5, p = 0.5) # 시행횟수 5, 성공확률 0.5인 이항분포에서 3개의 난수 발생
[1] 1 3 2
>
> # 4.2.2 포아송 분포
>
> dpois(3, lambda = 3) # 람다가 3인 포아송분포에서 x = 3의 확률
[1] 0.2240418
> ppois(3, lambda = 3) # 람다가 3인 포아송분포에서 x <= 3의 누적확률
[1] 0.6472319
> qpois(0.6, lambda = 3) # 람다가 3인 포아송분포에서 누적확률이 0.6이 되는 x
[1] 3
> rpois(10, lambda = 3) # 람다가 3인 포아송분포에서 10개의 난수 발생
[1] 3 3 6 3 2 2 6 3 6 5
>
> plot(table(rpois(100, 3)) # 람다가 3인 포아송분포에서 100개의 난수 발생하여 빈도분석
+ , type = "h" # 그래프 형태는 수직선
+ , col = "red" # 색깔은 빨강
+ , lwd = 10, # 선두께는 10
+ , main = "rpois(100, lambda = 3)") # 제목
>
>
> # 4.3 정규성 검정
>
> # 4.3.1 히스토그램을 이용한 정규성 검토
>
> op = par(mfrow=c(2,3)) # 2 by 3개의 그래프를 그릴 것
>
> x = seq(6.5, 13.5, by = 0.01) # 정규분포를 그리기 위한 기초자료
>
> for (i in 1:6) { # 6번 반복
+ y <- rnorm(1000, mean = 10, sd = 1) # 평균이 10이고, 표준편차가 1인 정규분포에서 1000개 난수 발생
+ hist(y, prob = TRUE, xlim = c(6.5, 13.5), ylim = c(0, 0.5)) # 히스토그램, 확률, x축 한계, y축 한계
+ lines(x, dnorm(x, mean = 10, sd = 1)) # 평균이 10이고 표준편차가 1인 정규분포 그래프 추가
+ }
>
> par(op) # 2 by 3 그래프 마지막
>
> # 4.3.2 정규확률그림
>
> x = rnorm(20) # 평균이 0이고 표준편차가 1인 정규분포에서 20 난수발생
> hist(x) # 히스토그램
> qqnorm(x) # Q-Q 플롯
> qqline(x) # Q-Q 라인
>
> # Shapiro 검정 - 정규성 검정
>
> x = rnorm(10) # 10개의 난수 발생
> shapiro.test(x) # Shapiro 검정

Shapiro-Wilk normality test

data: x
W = 0.9265, p-value = 0.4144

>
>
> # 여러 분포에서 발생한 난수를 이용하여 정규성 검정
>
> op = par(mfrow=c(2,2)) # 2 by 2 분할 그래프
>
> x = rnorm(100) # 정규분포에서 난수 발생
> qqnorm(x, sub = "Normal") # Q-Q Plot
> qqline(x) # Q-Q Line
>
> x = runif(100, 0, 1) # 균일분포에서 난수 발생
> qqnorm(x, sub = "Uniform") # Q-Q Plot
> qqline(x) # Q-Q Line
>
> x = rexp(100, 1) # 지수분포에서 난수 발생
> qqnorm(x, sub = "Exponential") # Q-Q Plot
> qqline(x) # Q-Q Line
>
> x = rpois(100, 3) # 포아송분포에서 난수 발생
> qqnorm(x, sub = "Poisson") # Q-Q Plot
> qqline(x) # Q-Q Line
>
> par(op) # 2 by 2 분할 그래프 끝
>
> # 4.4 모의 실험
>
> # 4.4.1 이항분포의 정규분포 근사
>
> m = 100
> p = 0.5
>
> try.n = c(5, 10, 15, 30) # 시행 횟수
>
> op = par(mfrow=c(2,2)) # 2 by 2 분할 그래프 시작
>
> for(i in 1:4){ # 네 번 반복
+ n = try.n[i] # 시행 횟수 조절
+ res = rbinom(m, n, p) # 이항분포에서 m개의 난수발생
+ hist(res, prob = TRUE) # 히스토그램, 확률
+ curve(dnorm(x, n*p, sqrt(n*p*(1-p))), add = TRUE) # 곡선추가
+ }
>
> par(op) # 2 by 2 분할 그래프 끝
>
> # 4.4.2 중심극한정리(central limit theorem)와 모의실험(simulation)
>
> m = 50
> mx = rep(0, m) # 0이 50개 있는 벡터 생성
> n.value = c(5, 10, 15, 30, 50) # 횟수
>
> # 균일분포 난수들의 표본평균들의 정규분포화
>
> plot(0, 0
+ , type = "n" # 그래프 형태 nothing, 즉 빈 그래프 생성
+ , xlim = c(0, 1) # x 축 값의 한계
+ , ylim = c(0, 10) # y 축 값의 한계
+ , ylab = "density" # y 축 라벨
+ , xlab = "mx" # x 축 라벨
+ , main = "uniform mean to normal") # 주 제목
>
> for(k in 1:length(n.value)){ # 5회 반복
+
+ n = n.value[k]
+
+ for(i in 1:m){ # 50회 반복
+ mx[i] = mean(runif(n, 0, 1)) # 균일분포에서 n개의 난수발생, 평균계산, 저장
+ }
+
+ lines(density(mx) # 라인추가
+ , lty = k # 라인의 형태
+ , col = k) # 색깔의 종류
+
+ }
>
>
> # 지수분포 난수들의 표본평균들의 정규분포화
>
> plot(0, 0
+ , type = "n" # 그래프 형태 nothing, 즉 빈 그래프 생성
+ , xlim = c(0, 3) # x 축 값의 한계
+ , ylim = c(0, 3.7) # y 축 값의 한계
+ , ylab = "density" # y 축 라벨
+ , xlab = "mx" # x 축 라벨
+ , main = "exponential to normal") # 주 제목
>
> for(k in 1:length(n.value)){ # 5회 반복
+
+ n = n.value[k]
+
+ for(i in 1:m){ # 50회 반복
+ mx[i] = mean(rexp(n, rate = 1)) # 지수분포에서 n개의 난수발생, 평균계산, 저장
+ }
+
+ lines(density(mx) # 라인추가
+ , lty = k # 라인의 형태
+ , col = k) # 색깔의 종류
+
+ }
>