> ##### 3. Hyper-Geometric Distributions : 초기하분포 #####
> N <- 100; m <- 10; n <- 5
> x <- 0:min(m, n)
> hg1 <- choose(m, x)*choose(N-m, n-x)/choose(N, n) # (1) Mass Function
> hg1
[1] 5.837524e-01 3.393909e-01 7.021881e-02 6.383528e-03 2.510376e-04 3.347168e-06
> hg2 <- dhyper(x, m, N-m, n) # (2) dhyper() 함수
> hg2
[1] 5.837524e-01 3.393909e-01 7.021881e-02 6.383528e-03 2.510376e-04 3.347168e-06
> format(hg1, scientific=F, digit=3) # (3) format() 함수: 소수 형태로 표시
[1] "0.58375237" "0.33939091" "0.07021881" "0.00638353" "0.00025104" "0.00000335"
> format(hg2, scientific=F, digit=3)
[1] "0.58375237" "0.33939091" "0.07021881" "0.00638353" "0.00025104" "0.00000335"
> cumsum(hg1) # (4) cumsum() 함수 이용
[1] 0.5837524 0.9231433 0.9933621 0.9997456 0.9999967 1.0000000
> phyper(x, m, N-m, n) # (5) hyper() 함수
[1] 0.5837524 0.9231433 0.9933621 0.9997456 0.9999967 1.0000000
> p <- m/N
> n * p # E(X)
[1] 0.5
> n * p * (1-p) * (N-n)/(N-1) # V(X)
[1] 0.4318182
> sapply(10^(0:5), # 난수 발생
+ FUN=function(x) {set.seed(1);
+ rnd <- rhyper(x, m, N-m, n);
+ data.frame(E=mean(rnd), V=var(rnd))})
[,1] [,2] [,3] [,4] [,5] [,6]
E 0 0.6 0.48 0.508 0.4986 0.50093
V NA 0.4888889 0.3329293 0.4684044 0.4366417 0.4332635
> ##### report. 초기하분포 -> 이항분포 #####
>
> n <- 5; p <- 0.4 # 상수값 고정
> N <- 20; # N의 초기값 설정
> m <- N*p; # m을 변수 값으로 설정
> x <- 0:min(m,n) # x의 범위 지정, 초기하분포는 비복원추출이므로 0부터 m이나 n 중에
> # 작은 값까지가 범위
> # 초기하분포와 이항분포의 그래프 생성
> # 원활한 비교를 위해 초기하분포 그래프 각각에 이항분포 삽입
> mypar <- par(mfrow=c(2,3)) # 화면분할 6등분
> # 초기하분포 변수인 hg.1 생성
> N <- 20; m <- N*p; hg.1 <- dhyper(x, m, N-m, n)
> # 초기하분포 그래프 생성: x가 0부터 5까지이므로 별도 지정하며
> # 초기하분포의 최대값이 약 0.47므로 y축 범위를 0부터 0.5까지로 고정
> # type(그래프유형), cex(점크기), lwd(선굵기), pch(점모양) 지정
> plot(hg.1, type="o", xlab="X", ylab="P(X=x)", x=0:5, ylim=c(0,0.5), cex=2, col="Red",
+ lwd=2, pch=1, main="초기하분포의 이항 근사; H(N,m,n), N=20, n=5, p=0.4")
> # 이항분포 결과에 해당하는 점 생성
> points(x, dbinom(x, n, p), pch=16)
> # 막대(type='h')형 그래프, lty(선종류)는 1~6까지를 모두 활용
> points(x, dbinom(x, n, p), type='h', lty=1, lwd=2)
> # legend() : 범례 함수, x=2가 좌측이고 y=0.5가 상단인 범례 상자 생성
> # col(색상), pch(점모양), lwd(선굵기)는 각각의 분포와 같게끔 설정
> # legend=c("",expression(""))으로 범례의 이름 설정
> legend(2, 0.5, col=c("Red", "Black"), lty=1, lwd=2, pch=c(1, 16),
+ legend=c("X~H(N=20,m=8,n=5)", "X~B(n=5,p=0.4)"))
> # N의 크기를 무한으로 향하게 늘려 m의 값을 조정해주면서 위 과정을 반복
> N <- 50; m <- N*p; hg.2 <- dhyper(x, m, N-m, n)
> plot(hg.2, type="o", xlab="X", ylab="P(X=x)", x=0:5, ylim=c(0,0.5), cex=2, col="Green",
+ lwd=2, pch=3, main="초기하분포의 이항 근사; H(N,m,n), N=50, n=5, p=0.4")
> points(x, dbinom(x, n, p), pch=16)
> points(x, dbinom(x, n, p), type='h', lty=2, lwd=2)
> legend(2, 0.5, col=c("Green", "Black"), lty=1:2, lwd=2, pch=c(3, 16),
+ legend=c("X~H(N=50,m=20,n=5)", "X~B(n=5,p=0.4)"))
> N <- 100; m <- N*p; hg.3 <- dhyper(x, m, N-m, n)
> plot(hg.3, type="o", xlab="X", ylab="P(X=x)", x=0:5, ylim=c(0,0.5), cex=2, col="Blue",
+ lwd=2, pch=5, main="초기하분포의 이항 근사; H(N,m,n), N=100, n=5, p=0.4")
> points(x, dbinom(x, n, p), pch=16)
> points(x, dbinom(x, n, p), type='h', lty=3, lwd=2)
> legend(2, 0.5, col=c("Blue", "Black"), lty=c(1,3), lwd=2, pch=c(5, 16),
+ legend=c("X~H(N=100,m=40,n=5)", "X~B(n=5,p=0.4)"))
> # col(색깔) 지정방식을 [col=""]에서 [col=고유번호]로 변경
> N <- 500; m <- N*p; hg.4 <- dhyper(x, m, N-m, n)
> plot(hg.4, type="o", xlab="X", ylab="P(X=x)", x=0:5, ylim=c(0,0.5), cex=2, col=6,
+ lwd=2, pch=15, main="초기하분포의 이항 근사; H(N,m,n), N=500, n=5, p=0.4")
> points(x, dbinom(x, n, p), pch=16)
> points(x, dbinom(x, n, p), type='h', lty=4, lwd=2)
> legend(2, 0.5, col=c(6, 1), lty=c(1,4), lwd=2, pch=c(15, 16),
+ legend=c("X~H(N=500,m=200,n=5)", "X~B(n=5,p=0.4)"))
> N <- 1000; m <- N*p; hg.5 <- dhyper(x, m, N-m, n)
> plot(hg.5, type="o", xlab="X", ylab="P(X=x)", x=0:5, ylim=c(0,0.5), cex=2, col=7,
+ lwd=2, pch=19, main="초기하분포의 이항 근사; H(N,m,n), N=1000, n=5, p=0.4")
> points(x, dbinom(x, n, p), pch=16)
> points(x, dbinom(x, n, p), type='h', lty=5, lwd=2)
> legend(2, 0.5, col=c(7, 1), lty=c(1,5), lwd=2, pch=c(19, 16),
+ legend=c("X~H(N=1000,m=400,n=5)", "X~B(n=5,p=0.4)"))
> N <- 10000; m <- N*p; hg.6 <- dhyper(x, m, N-m, n)
> plot(hg.6, type="o", xlab="X", ylab="P(X=x)", x=0:5, ylim=c(0,0.5), cex=2, col=8,
+ lwd=2, pch=23, main="초기하분포의 이항 근사; H(N,m,n), N=10000, n=5, p=0.4")
> points(x, dbinom(x, n, p), pch=16)
> points(x, dbinom(x, n, p), type='h', lty=6, lwd=2)
> legend(2, 0.5, col=c(8, 1), lty=c(1,6), lwd=2, pch=c(23, 16),
+ legend=c("X~H(N=10000,m=4000,n=5)", "X~B(n=5,p=0.4)"))
> par(mypar)
> # 초기하분포와 이항분포의 차이(gap) 계산
> binom <- dbinom(x=0:5, 5, p=0.4) # 이항분포 변수값 생성
> gap1 <- hg.1-binom # N=20일 때 초기하분포와 이항분포의 차이 계산
> gap2 <- hg.2-binom # N=50일 때 초기하분포와 이항분포의 차이 계산
> gap3 <- hg.3-binom # N=100일 때 초기하분포와 이항분포의 차이 계산
> gap4 <- hg.4-binom # N=500일 때 초기하분포와 이항분포의 차이 계산
> gap5 <- hg.5-binom # N=1,000일 때 초기하분포와 이항분포의 차이 계산
> gap6 <- hg.6-binom # N=10,000일 때 초기하분포와 이항분포의 차이 계산
>
> # 계산한 차이를 데이터 프레임으로 표시, x가 0부터 시작이므로 x도 데이터 프레임에 포함
> gap <- data.frame(x=0:min(m,n), gap1, gap2, gap3, gap4, gap5, gap6)
> # 소수점으로 표기하기 위해 format()을 활용, scientific=F로 표시를 소수 형태로 변경
> format(gap, scientific=F, digit=3)
x gap1 gap2 gap3 gap4 gap5 gap6
1 0 -0.02668 -0.010501 -0.005218 -0.00103818 -0.00051875 -0.0000518435
2 1 -0.00378 -0.000511 -0.000121 -0.00000465 -0.00000116 -0.0000000115
3 2 0.05172 0.018481 0.008929 0.00173922 0.00086679 0.0000864279
4 3 0.00799 0.003652 0.001878 0.00038243 0.00019161 0.0000191962
5 4 -0.02262 -0.008199 -0.003967 -0.00077301 -0.00038525 -0.0000384125
6 5 -0.00663 -0.002923 -0.001500 -0.00030581 -0.00015325 -0.0000153565
> # 결론 : N의 값이 커질수록 초기하분포와 이항분포의 차이는 점점 줄어든다.
> # hg.1에서 쓴 변수값을 써야 하므로 N, m 등의 변수 재설정
> N <- 20; m <- N*p; hg.1 <- dhyper(x, m, N-m, n)
> # hg.1 plot 생성
> plot(hg.1, type="o", xlab="X", ylab="P(X=x)", x=0:5, ylim=c(0,0.5), cex=2, col="Red",
+ lwd=2, pch=1, main="초기하분포의 이항 근사; H(N,m,n), n=5, p=0.4")
> N <- 50; m <- N*p; hg.2 <- dhyper(x, m, N-m, n)
> # hg.2부터는 point() 함수를 이용하여 그래프 생성
> points(hg.2, type="o", x=0:5, ylim=c(0,0.5), cex=2, col="Green", lwd=2, pch=3)
> N <- 100; m <- N*p; hg.3 <- dhyper(x, m, N-m, n)
> points(hg.3, type="o", x=0:5, ylim=c(0,0.5), cex=2, col="Blue", lwd=2, pch=5)
> N <- 500; m <- N*p; hg.4 <- dhyper(x, m, N-m, n)
> points(hg.4, type="o", x=0:5, ylim=c(0,0.5), cex=2, col=6, lwd=2, pch=15)
> N <- 1000; m <- N*p; hg.5 <- dhyper(x, m, N-m, n)
> points(hg.5, type="o", x=0:5, ylim=c(0,0.5), cex=2, col=7, lwd=2, pch=19)
> N <- 10000; m <- N*p; hg.6 <- dhyper(x, m, N-m, n)
> points(hg.6, type="o", x=0:5, ylim=c(0,0.5), cex=2, col=8, lwd=2, pch=23)
> # 합쳐진 초기하분포 그래프에 이항분포(type='h') 그래프 겹치기
> points(x, dbinom(x, n, p), pch=16)
> points(x, dbinom(x, n, p), type='h', lty=6, lwd=2)
> # 범례 생성 : 초기하분포 6개와 이항분포 하나의 그래프 모두 포함
> legend(3, 0.5, col=c("Red", "Green", "Blue", 6, 7, 8, 1), lwd=2, pch=c(1,3,5,15,19,23,16),
+ legend=c("X~H(N=20,m=8,n=5)", "X~H(N=50,m=20,n=5)",
+ "X~H(N=100,m=40,n=5)", "X~H(N=500,m=200,n=5)",
+ "X~H(N=1000,m=400,n=5)", "X~H(N=10000,m=4000,n=5)",
+ "X~B(n=5,p=0.4)"))
> cov <- var(gap) # 참고. 공분산행렬
> format(cov, scientific=F, digit=1)
x gap1 gap2
x "3.500000000000000000" "0.000000000000000000" "0.000000000000000031"
gap1 "0.000000000000000000" "0.000804004122784053" "0.000294366944192417"
gap2 "0.000000000000000031" "0.000294366944192417" "0.000108233228176942"
gap3 "0.000000000000000032" "0.000143222334736306" "0.000052726865003160"
gap4 "0.000000000000000005" "0.000028045532954817" "0.000010334865353541"
gap5 "0.000000000000000009" "0.000013986353020685" "0.000005154623311419"
gap6 "0.000000000000000022" "0.000001395376833478" "0.000000514315876031"
gap3 gap4 gap5
x "0.000000000000000032" "0.000000000000000005" "0.000000000000000009"
gap1 "0.000143222334736306" "0.000028045532954817" "0.000013986353020685"
gap2 "0.000052726865003160" "0.000010334865353541" "0.000005154623311419"
gap3 "0.000025696109748295" "0.000005038086417808" "0.000002512887490211"
gap4 "0.000005038086417808" "0.000000988006082314" "0.000000492809140497"
gap5 "0.000002512887490211" "0.000000492809140497" "0.000000245809872250"
gap6 "0.000000250737773870" "0.000000049174047984" "0.000000024527755522"
gap6
x "0.000000000000000022"
gap1 "0.000001395376833478"
gap2 "0.000000514315876031"
gap3 "0.000000250737773870"
gap4 "0.000000049174047984"
gap5 "0.000000024527755522"
gap6 "0.000000002447470378"
> var(gap1)
[1] 0.0008040041
> var(gap2)
[1] 0.0001082332
> var(gap3)
[1] 2.569611e-05
> var(gap4)
[1] 9.880061e-07
> var(gap5)
[1] 2.458099e-07
> var(gap6)
[1] 2.44747e-09
> ##### Hypergeometric --> Binomial --> Poisson (초기하 -> 이항 -> 포아송) #####
> N <- 100; m <- 30; n <- 5
> x <- 0:n
> eps <- 5
> # round() 함수를 사용하여 소수점 몇 째 자리만큼 표시
> round(p.hyper <- dhyper(x, m, N-m, n), eps) # 초기하분포의 pdf
[1] 0.16076 0.36536 0.31628 0.13023 0.02548 0.00189
> round(p.binom <- dbinom(x, n, m/N), eps) # 이항분포의 pdf
[1] 0.16807 0.36015 0.30870 0.13230 0.02835 0.00243
> round(p.pois <- dpois(x, n*m/N), eps) # 포아송분포의 pdf
[1] 0.22313 0.33470 0.25102 0.12551 0.04707 0.01412
> N.times <- n.times <- 10^(0:5)
> p.hyper <- sapply(N.times, FUN=function(N.times)
+ dhyper(x, m*N.times, (N-m)*N.times, n))
> dimnames(p.hyper) = list(paste("X=", x, sep=""),
+ paste("N*", N.times, sep=""))
> round(p.binom, eps) - round(p.hyper, eps)
N*1 N*10 N*100 N*1000 N*10000 N*1e+05
X=0 0.00731 0.00072 7e-05 1e-05 0 0
X=1 -0.00521 -0.00052 -5e-05 -1e-05 0 0
X=2 -0.00758 -0.00074 -7e-05 -1e-05 0 0
X=3 0.00207 0.00019 2e-05 0e+00 0 0
X=4 0.00287 0.00028 3e-05 0e+00 0 0
X=5 0.00054 0.00006 1e-05 0e+00 0 0
> p.binom <- sapply(n.times, FUN=function(n.times)
+ dbinom(x, n*n.times, m/N/n.times))
> dimnames(p.binom) = list(paste0("X=", x),
+ paste("N*", n.times, sep=""))
> round(p.pois, eps) - round(p.binom, eps)
N*1 N*10 N*100 N*1000 N*10000 N*1e+05
X=0 0.05506 0.00506 0.00050 5e-05 0e+00 0
X=1 -0.02545 -0.00251 -0.00025 -2e-05 0e+00 0
X=2 -0.05768 -0.00450 -0.00044 -5e-05 -1e-05 0
X=3 -0.00679 -0.00093 -0.00009 -1e-05 0e+00 0
X=4 0.01872 0.00112 0.00011 1e-05 0e+00 0
X=5 0.01169 0.00105 0.00010 1e-05 0e+00 0
'데이터 [Data] > R' 카테고리의 다른 글
| R plot: 이산형 분포의 확률밀도함수 (0) | 2021.06.08 |
|---|---|
| 워드클라우드 자체실습: wordcloud2() (0) | 2021.06.07 |
| R Distributions: 포아송분포 (0) | 2021.06.06 |
| R Distributions: 이항분포의 누적분포함수 (0) | 2021.06.05 |
| R Distributions: 이항분포 (0) | 2021.06.04 |
댓글