정리
\begin{align*} \text{Geometric Distribution: } \;\;\; \text{X} & \thicksim Geo(p) \\ p(X = k) & = q^{k-1} \cdot p \\ E\left[ X \right] & = \frac{1}{p} \\ V\left[ X \right] & = \frac{q}{p^2} \\ \\ \text{Binomial Distribution: } \;\;\; \text{X} & \thicksim B(n, p) \\ p(X = r) & = \binom{n}{r} \cdot p^{r} \cdot q^{n-r} \\ E\left[ X \right] & = {n}{p} \\ V\left[ X \right] & = {n}{p}{q} \\ \\ \text{Poisson Distribution: } \;\;\; \text{X} & \thicksim P( \lambda ) \\ P(X=r) & = e^{- \lambda} \cdot \dfrac{\lambda^{r}} {r!} \\ E\left[ X \right] & = \lambda \\ V\left[ X \right] & = \lambda \\ \end{align*}

The probability of Chad making a clear run down the slope is 0.2, and he's going to keep on trying until he succeeds. After he’s made his first successful run down the slopes, he’s going to stop snowboarding, and head back to the lodge triumphantly

P(X = 1) = P(success in the first trial) = 0.2
P(X = 2) = P(success in the second trial union failure in the first trial) = 0.8 * 0.2 = 0.16
1회 혹은 2회에서 성공할 확률
P(X <= 2) = P(X = 1) + P(X = 2) = 0.2 + 0.16 = 0.36

$P(X = r) = 0.8^{r-1} × 0.2$
$P(X = r) = q^{r-1} × p $

This formula is called the geometric distribution.

• You run a series of independent trials. * (각 시행이 독립적임 = 이번 시행이 이전 시행과 상관없이 일어남)
• There can be either a success or failure for each trial, and the probability of success is the same for each trial. (성공/실패만 일어나는 상황에서, 성공하는 확률은 p로 실행햇수동안 일정)
• The main thing you’re interested in is how many trials are needed in order to get the first successful outcome. (성공하면 중단하고 성공할 때까지의 확률을 분포로 봄)

$ P(X=r) = {p \cdot q^{r-1}} $
$ P(X=r) = {p \cdot (1-p)^{r-1}} $

p = 0.20
n = 29
## geometric . . . . 
## note that it starts with 0 rather than 1
## since the function uses p * q^(r), 
## rather than p * q^(r-1)
dgeom(x = 0:n, prob = p)
hist(dgeom(x = 0:n, prob = p))

> p = 0.20
> n = 29
> # exact
> dgeom(0:n, prob = p)
 [1] 0.2000000000 0.1600000000 0.1280000000 0.1024000000 0.0819200000 0.0655360000 0.0524288000
 [8] 0.0419430400 0.0335544320 0.0268435456 0.0214748365 0.0171798692 0.0137438953 0.0109951163
[15] 0.0087960930 0.0070368744 0.0056294995 0.0045035996 0.0036028797 0.0028823038 0.0023058430
[22] 0.0018446744 0.0014757395 0.0011805916 0.0009444733 0.0007555786 0.0006044629 0.0004835703
[29] 0.0003868563 0.0003094850 
> 
> hist(dgeom(x = 0:n, prob = p))

r번 시도한 이후, 그 이후 어디서든지 간에 성공을 얻을 확률
$$ P(X > r) = q^{r} $$

예, 20번 시도 후에 어디선가 성공할 확률은?

Solution.

• 21번째 성공 + 22번째 + 23번째 + . . . .
• 위는 구할 수 없음
• 따라서
• 전체 확률이 1이고 20번째까지 성공한 확률을 1에서 빼면 원하는 확률이 됨
• 1 - (1번째 성공 + 2번째 성공 + . . . + 20번째 성공)
• 그런데 이것은
• 20번까지는 실패하는 확률 = $q^{r} $ 과 같다

p <- .2
q <- 1-p
n <- 19
s <- dgeom(x = 0:n, prob = p)
# 20번째까지 성공할 확률을 모두 더한 확률
sum(s)
# 따라서 아래는 20번 이후 어디서든지 간에서 성공할 확률
1-sum(s)
## 혹은 (교재가 이야기하는) 20번까지 실패하는 확률
q^20

> p <- .2
> q <- 1-p
> n <- 19
> s <- dgeom(x = 0:n, prob = p)
> # 20번째까지 성공할 확률
> sum(s)
[1] 0.9884708
> # 따라서 아래는 20번 이후 어디서든지 간에서 성공할 확률
> 1-sum(s)
[1] 0.01152922
> ## 혹은 (교재가 이야기하는) 20번까지 실패하는 확률
> q^20
[1] 0.01152922
> 

그렇다면
r 번 이전에 성공이 있을 확률은? = r 번까지의 실패할 확률의 보수
$$ P(X \le r) = 1 - q^{r} $$

혹은 1번째 성공 + 2번째 성공 + . . . + r 번째 성공으로 구해도 된다

# r = 20 이라고 하면 
p <- .2
q <- 1-p
n <- 19
s <- dgeom(x = 0:n, prob = p)
sum(s)

Note that
$$P(X > r) + P(X \le r) = 1 $$

X가 성공할 확률 p를 가진 Geometric distribution을 따른다 :: $X \sim \text{Geo}(p)$

Reminding . . . Expected value in discrete probability
$E(X) = \sum x*P(X=x)$

p <- .2
q <- 1-p
trial <- c(1:8)
px <- q^(trial-1)*p
px
## npx <- trial*(q^(trial-1))*p
## 위는 아래와 같음
npx <- trial*px
npx
## plex <- cumsum(trial*(q^(trial-1))*p)
## 위는 아래와 같음
plex <- cumsum(npx)
plex
sumgeod <- data.frame(trial,px,npx,plex)
round(sumgeod,3)

> p <- .2
> q <- 1-p
> trial <- c(1,2,3,4,5,6,7,8)
> px <- q^(trial-1)*p
> px
[1] 0.20000000 0.16000000 0.12800000 0.10240000 0.08192000 0.06553600 0.05242880 0.04194304
> npx <- trial*(q^(trial-1))*p
> npx
[1] 0.2000000 0.3200000 0.3840000 0.4096000 0.4096000 0.3932160 0.3670016 0.3355443
> plex <- cumsum(trial*(q^(trial-1))*p)
> plex
[1] 0.200000 0.520000 0.904000 1.313600 1.723200 2.116416 2.483418 2.818962
> sumgeod <- data.frame(trial,px,npx,plex)
> round(sumgeod,3)
  trial    px   npx  plex
1     1 0.200 0.200 0.200
2     2 0.160 0.320 0.520
3     3 0.128 0.384 0.904
4     4 0.102 0.410 1.314
5     5 0.082 0.410 1.723
6     6 0.066 0.393 2.116
7     7 0.052 0.367 2.483
8     8 0.042 0.336 2.819
> 

> p <- .2
> q <- 1-p
> trial <- c(1:100)
> px <- q^(trial-1)*p
> round(px, 3)
  [1] 0.200 0.160 0.128 0.102 0.082 0.066 0.052 0.042 0.034 0.027 0.021 0.017
 [13] 0.014 0.011 0.009 0.007 0.006 0.005 0.004 0.003 0.002 0.002 0.001 0.001
 [25] 0.001 0.001 0.001 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
 [37] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
 [49] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
 [61] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
 [73] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
 [85] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
 [97] 0.000 0.000 0.000 0.000
> npx <- trial*(q^(trial-1))*p
> round(npx, 3)
  [1] 0.200 0.320 0.384 0.410 0.410 0.393 0.367 0.336 0.302 0.268 0.236 0.206
 [13] 0.179 0.154 0.132 0.113 0.096 0.081 0.068 0.058 0.048 0.041 0.034 0.028
 [25] 0.024 0.020 0.016 0.014 0.011 0.009 0.008 0.006 0.005 0.004 0.004 0.003
 [37] 0.002 0.002 0.002 0.001 0.001 0.001 0.001 0.001 0.000 0.000 0.000 0.000
 [49] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
 [61] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
 [73] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
 [85] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
 [97] 0.000 0.000 0.000 0.000
> plex <- cumsum(trial*(q^(trial-1))*p)
> round(plex, 3)
  [1] 0.200 0.520 0.904 1.314 1.723 2.116 2.483 2.819 3.121 3.389 3.626 3.832
 [13] 4.010 4.164 4.296 4.409 4.505 4.586 4.654 4.712 4.760 4.801 4.835 4.863
 [25] 4.887 4.906 4.923 4.936 4.947 4.957 4.964 4.971 4.976 4.980 4.984 4.987
 [37] 4.989 4.991 4.993 4.994 4.995 4.996 4.997 4.997 4.998 4.998 4.999 4.999
 [49] 4.999 4.999 4.999 4.999 5.000 5.000 5.000 5.000 5.000 5.000 5.000 5.000
 [61] 5.000 5.000 5.000 5.000 5.000 5.000 5.000 5.000 5.000 5.000 5.000 5.000
 [73] 5.000 5.000 5.000 5.000 5.000 5.000 5.000 5.000 5.000 5.000 5.000 5.000
 [85] 5.000 5.000 5.000 5.000 5.000 5.000 5.000 5.000 5.000 5.000 5.000 5.000
 [97] 5.000 5.000 5.000 5.000
> sumgeod <- data.frame(trial,px,npx,plex)
> round(sumgeod,3)
    trial    px   npx  plex
1       1 0.200 0.200 0.200
2       2 0.160 0.320 0.520
3       3 0.128 0.384 0.904
4       4 0.102 0.410 1.314
5       5 0.082 0.410 1.723
6       6 0.066 0.393 2.116
7       7 0.052 0.367 2.483
8       8 0.042 0.336 2.819
9       9 0.034 0.302 3.121
10     10 0.027 0.268 3.389
11     11 0.021 0.236 3.626
12     12 0.017 0.206 3.832
13     13 0.014 0.179 4.010
14     14 0.011 0.154 4.164
15     15 0.009 0.132 4.296
16     16 0.007 0.113 4.409
17     17 0.006 0.096 4.505
18     18 0.005 0.081 4.586
19     19 0.004 0.068 4.654
20     20 0.003 0.058 4.712
21     21 0.002 0.048 4.760
22     22 0.002 0.041 4.801
23     23 0.001 0.034 4.835
24     24 0.001 0.028 4.863
25     25 0.001 0.024 4.887
26     26 0.001 0.020 4.906
27     27 0.001 0.016 4.923
28     28 0.000 0.014 4.936
29     29 0.000 0.011 4.947
30     30 0.000 0.009 4.957
31     31 0.000 0.008 4.964
32     32 0.000 0.006 4.971
33     33 0.000 0.005 4.976
34     34 0.000 0.004 4.980
35     35 0.000 0.004 4.984
36     36 0.000 0.003 4.987
37     37 0.000 0.002 4.989
38     38 0.000 0.002 4.991
39     39 0.000 0.002 4.993
40     40 0.000 0.001 4.994
41     41 0.000 0.001 4.995
42     42 0.000 0.001 4.996
43     43 0.000 0.001 4.997
44     44 0.000 0.001 4.997
45     45 0.000 0.000 4.998
46     46 0.000 0.000 4.998
47     47 0.000 0.000 4.999
48     48 0.000 0.000 4.999
49     49 0.000 0.000 4.999
50     50 0.000 0.000 4.999
51     51 0.000 0.000 4.999
52     52 0.000 0.000 4.999
53     53 0.000 0.000 5.000
54     54 0.000 0.000 5.000
55     55 0.000 0.000 5.000
56     56 0.000 0.000 5.000
57     57 0.000 0.000 5.000
58     58 0.000 0.000 5.000
59     59 0.000 0.000 5.000
60     60 0.000 0.000 5.000
61     61 0.000 0.000 5.000
62     62 0.000 0.000 5.000
63     63 0.000 0.000 5.000
64     64 0.000 0.000 5.000
65     65 0.000 0.000 5.000
66     66 0.000 0.000 5.000
67     67 0.000 0.000 5.000
68     68 0.000 0.000 5.000
69     69 0.000 0.000 5.000
70     70 0.000 0.000 5.000
71     71 0.000 0.000 5.000
72     72 0.000 0.000 5.000
73     73 0.000 0.000 5.000
74     74 0.000 0.000 5.000
75     75 0.000 0.000 5.000
76     76 0.000 0.000 5.000
77     77 0.000 0.000 5.000
78     78 0.000 0.000 5.000
79     79 0.000 0.000 5.000
80     80 0.000 0.000 5.000
81     81 0.000 0.000 5.000
82     82 0.000 0.000 5.000
83     83 0.000 0.000 5.000
84     84 0.000 0.000 5.000
85     85 0.000 0.000 5.000
86     86 0.000 0.000 5.000
87     87 0.000 0.000 5.000
88     88 0.000 0.000 5.000
89     89 0.000 0.000 5.000
90     90 0.000 0.000 5.000
91     91 0.000 0.000 5.000
92     92 0.000 0.000 5.000
93     93 0.000 0.000 5.000
94     94 0.000 0.000 5.000
95     95 0.000 0.000 5.000
96     96 0.000 0.000 5.000
97     97 0.000 0.000 5.000
98     98 0.000 0.000 5.000
99     99 0.000 0.000 5.000
100   100 0.000 0.000 5.000
> plot(npx, type="l")
> plot(plex, type="l")

위에서 $X \sim \text{Geo}(p)$ 일때, 기대값은 $E(X) = \displaystyle \frac{1}{p}$ 임을 알 수 있다. 아래는 그 증명이다.

기하분포의 기대값(expected value)

기하분포의 분산 증명
아래는 R에서의 시뮬레이션
\begin{eqnarray*} Var(X) & = & E((X-E(X))^{2}) \\ & = & E(X^{2} - 2XE(X) + E(X)^{2}) \\ & = & E(X^{2}) - E(2X)(X) + E(X)^{2} \\ & = & E(X^{2}) - 2E(X)E(X) + E(X)^{2} \\ & = & E(X^{2}) - 2E(X)^{2} + E(X)^{2} \\ & = & E(X^{2}) - E(X)^{2} \end{eqnarray*}

Var(X) 계산의 앞부분 $E[X^{2}]$ 부분은
$E(X^{2}) = \sum {x^{2} * P(X=x)}$
왜냐하면, $E(X) = \sum {x * P(X=x)}$ 이므로

이를 구한 후에 여기에서 $E[X]^{2}$ 를 빼준 값이 Var(X)
참고로 $E[X]^{2} = 25$

\begin{eqnarray*} E(X)^{2} & = & \left(\frac{1}{p}\right)^{2} \\ & = & (5^{2}) \\ & = & 25 \end{eqnarray*}

p <- .2
q <- 1-p
# 100 trials
trial <- c(1:100)   
# px = probability of geometric distribution
px <- q^(trial-1)*p 
px

# for variance, Var(X) = E(X^2)-E(X)^2
# E(X^2) = sum(x^2*p(X=x))
# 위에서 x 는 trial 에 해당하는 숫자
# after summing up the above
# subtract E(X)^2, that is, squared value
# of expected value of X

# for the X^2 part in E(X^2)
# that is, x^2*p(x) part
nspx <- (trial^2)*px
nspx

# summing up the nspx
# sum of x^2*p(x) 
plex <- cumsum(nspx)
plex

sumgeod <- data.frame(trial,px,nspx,plex)
round(sumgeod,3)

plot(nspx, type="l")
plot(plex, type="l")

> p <- .2
> q <- 1-p
> trial <- c(1:100)
> px <- q^(trial-1)*p
> px
  [1] 2.000000e-01 1.600000e-01 1.280000e-01 1.024000e-01 8.192000e-02 6.553600e-02 5.242880e-02
  [8] 4.194304e-02 3.355443e-02 2.684355e-02 2.147484e-02 1.717987e-02 1.374390e-02 1.099512e-02
 [15] 8.796093e-03 7.036874e-03 5.629500e-03 4.503600e-03 3.602880e-03 2.882304e-03 2.305843e-03
 [22] 1.844674e-03 1.475740e-03 1.180592e-03 9.444733e-04 7.555786e-04 6.044629e-04 4.835703e-04
 [29] 3.868563e-04 3.094850e-04 2.475880e-04 1.980704e-04 1.584563e-04 1.267651e-04 1.014120e-04
 [36] 8.112964e-05 6.490371e-05 5.192297e-05 4.153837e-05 3.323070e-05 2.658456e-05 2.126765e-05
 [43] 1.701412e-05 1.361129e-05 1.088904e-05 8.711229e-06 6.968983e-06 5.575186e-06 4.460149e-06
 [50] 3.568119e-06 2.854495e-06 2.283596e-06 1.826877e-06 1.461502e-06 1.169201e-06 9.353610e-07
 [57] 7.482888e-07 5.986311e-07 4.789049e-07 3.831239e-07 3.064991e-07 2.451993e-07 1.961594e-07
 [64] 1.569275e-07 1.255420e-07 1.004336e-07 8.034690e-08 6.427752e-08 5.142202e-08 4.113761e-08
 [71] 3.291009e-08 2.632807e-08 2.106246e-08 1.684997e-08 1.347997e-08 1.078398e-08 8.627183e-09
 [78] 6.901746e-09 5.521397e-09 4.417118e-09 3.533694e-09 2.826955e-09 2.261564e-09 1.809251e-09
 [85] 1.447401e-09 1.157921e-09 9.263367e-10 7.410694e-10 5.928555e-10 4.742844e-10 3.794275e-10
 [92] 3.035420e-10 2.428336e-10 1.942669e-10 1.554135e-10 1.243308e-10 9.946465e-11 7.957172e-11
 [99] 6.365737e-11 5.092590e-11
> ## for variance, Var(X) = E(X^2)-E(X)^2
> ## E(X^2) = sum(x^2*p(X=x))
> ## after summing up the above
> ## substract E(X)^2, that is, squared value
> ## of expected value of X
> nspx <- (trial^2)*(q^(trial-1))*p
> nspx
  [1] 2.000000e-01 6.400000e-01 1.152000e+00 1.638400e+00 2.048000e+00 2.359296e+00 2.569011e+00
  [8] 2.684355e+00 2.717909e+00 2.684355e+00 2.598455e+00 2.473901e+00 2.322718e+00 2.155043e+00
 [15] 1.979121e+00 1.801440e+00 1.626925e+00 1.459166e+00 1.300640e+00 1.152922e+00 1.016877e+00
 [22] 8.928224e-01 7.806662e-01 6.800208e-01 5.902958e-01 5.107712e-01 4.406535e-01 3.791191e-01
 [29] 3.253461e-01 2.785365e-01 2.379321e-01 2.028241e-01 1.725589e-01 1.465404e-01 1.242298e-01
 [36] 1.051440e-01 8.885318e-02 7.497677e-02 6.317987e-02 5.316912e-02 4.468865e-02 3.751613e-02
 [43] 3.145910e-02 2.635147e-02 2.205030e-02 1.843296e-02 1.539448e-02 1.284523e-02 1.070882e-02
 [50] 8.920298e-03 7.424542e-03 6.174844e-03 5.131698e-03 4.261739e-03 3.536834e-03 2.933292e-03
 [57] 2.431190e-03 2.013795e-03 1.667068e-03 1.379246e-03 1.140483e-03 9.425461e-04 7.785568e-04
 [64] 6.427752e-04 5.304151e-04 4.374889e-04 3.606772e-04 2.972193e-04 2.448202e-04 2.015743e-04
 [71] 1.658998e-04 1.364847e-04 1.122418e-04 9.227042e-05 7.582485e-05 6.228826e-05 5.115057e-05
 [78] 4.199022e-05 3.445904e-05 2.826955e-05 2.318457e-05 1.900845e-05 1.557992e-05 1.276608e-05
 [85] 1.045747e-05 8.563983e-06 7.011443e-06 5.738841e-06 4.696008e-06 3.841704e-06 3.142039e-06
 [92] 2.569180e-06 2.100268e-06 1.716542e-06 1.402607e-06 1.145833e-06 9.358629e-07 7.642068e-07
 [99] 6.239059e-07 5.092590e-07
> ## summing up the nspx
> plex <- cumsum(nspx)
> plex
  [1]  0.200000  0.840000  1.992000  3.630400  5.678400  8.037696 10.606707 13.291062 16.008971
 [10] 18.693325 21.291781 23.765682 26.088400 28.243443 30.222564 32.024004 33.650929 35.110095
 [19] 36.410735 37.563656 38.580533 39.473355 40.254022 40.934042 41.524338 42.035109 42.475763
 [28] 42.854882 43.180228 43.458765 43.696697 43.899521 44.072080 44.218620 44.342850 44.447994
 [37] 44.536847 44.611824 44.675004 44.728173 44.772862 44.810378 44.841837 44.868188 44.890239
 [46] 44.908671 44.924066 44.936911 44.947620 44.956540 44.963965 44.970140 44.975271 44.979533
 [55] 44.983070 44.986003 44.988434 44.990448 44.992115 44.993495 44.994635 44.995578 44.996356
 [64] 44.996999 44.997529 44.997967 44.998327 44.998625 44.998870 44.999071 44.999237 44.999373
 [73] 44.999486 44.999578 44.999654 44.999716 44.999767 44.999809 44.999844 44.999872 44.999895
 [82] 44.999914 44.999930 44.999943 44.999953 44.999962 44.999969 44.999974 44.999979 44.999983
 [91] 44.999986 44.999989 44.999991 44.999992 44.999994 44.999995 44.999996 44.999997 44.999997
[100] 44.999998
> sumgeod <- data.frame(trial,px,nspx,plex)
> round(sumgeod,3)
    trial    px  nspx   plex
1       1 0.200 0.200  0.200
2       2 0.160 0.640  0.840
3       3 0.128 1.152  1.992
4       4 0.102 1.638  3.630
5       5 0.082 2.048  5.678
6       6 0.066 2.359  8.038
7       7 0.052 2.569 10.607
8       8 0.042 2.684 13.291
9       9 0.034 2.718 16.009
10     10 0.027 2.684 18.693
11     11 0.021 2.598 21.292
12     12 0.017 2.474 23.766
13     13 0.014 2.323 26.088
14     14 0.011 2.155 28.243
15     15 0.009 1.979 30.223
16     16 0.007 1.801 32.024
17     17 0.006 1.627 33.651
18     18 0.005 1.459 35.110
19     19 0.004 1.301 36.411
20     20 0.003 1.153 37.564
21     21 0.002 1.017 38.581
22     22 0.002 0.893 39.473
23     23 0.001 0.781 40.254
24     24 0.001 0.680 40.934
25     25 0.001 0.590 41.524
26     26 0.001 0.511 42.035
27     27 0.001 0.441 42.476
28     28 0.000 0.379 42.855
29     29 0.000 0.325 43.180
30     30 0.000 0.279 43.459
31     31 0.000 0.238 43.697
32     32 0.000 0.203 43.900
33     33 0.000 0.173 44.072
34     34 0.000 0.147 44.219
35     35 0.000 0.124 44.343
36     36 0.000 0.105 44.448
37     37 0.000 0.089 44.537
38     38 0.000 0.075 44.612
39     39 0.000 0.063 44.675
40     40 0.000 0.053 44.728
41     41 0.000 0.045 44.773
42     42 0.000 0.038 44.810
43     43 0.000 0.031 44.842
44     44 0.000 0.026 44.868
45     45 0.000 0.022 44.890
46     46 0.000 0.018 44.909
47     47 0.000 0.015 44.924
48     48 0.000 0.013 44.937
49     49 0.000 0.011 44.948
50     50 0.000 0.009 44.957
51     51 0.000 0.007 44.964
52     52 0.000 0.006 44.970
53     53 0.000 0.005 44.975
54     54 0.000 0.004 44.980
55     55 0.000 0.004 44.983
56     56 0.000 0.003 44.986
57     57 0.000 0.002 44.988
58     58 0.000 0.002 44.990
59     59 0.000 0.002 44.992
60     60 0.000 0.001 44.993
61     61 0.000 0.001 44.995
62     62 0.000 0.001 44.996
63     63 0.000 0.001 44.996
64     64 0.000 0.001 44.997
65     65 0.000 0.001 44.998
66     66 0.000 0.000 44.998
67     67 0.000 0.000 44.998
68     68 0.000 0.000 44.999
69     69 0.000 0.000 44.999
70     70 0.000 0.000 44.999
71     71 0.000 0.000 44.999
72     72 0.000 0.000 44.999
73     73 0.000 0.000 44.999
74     74 0.000 0.000 45.000
75     75 0.000 0.000 45.000
76     76 0.000 0.000 45.000
77     77 0.000 0.000 45.000
78     78 0.000 0.000 45.000
79     79 0.000 0.000 45.000
80     80 0.000 0.000 45.000
81     81 0.000 0.000 45.000
82     82 0.000 0.000 45.000
83     83 0.000 0.000 45.000
84     84 0.000 0.000 45.000
85     85 0.000 0.000 45.000
86     86 0.000 0.000 45.000
87     87 0.000 0.000 45.000
88     88 0.000 0.000 45.000
89     89 0.000 0.000 45.000
90     90 0.000 0.000 45.000
91     91 0.000 0.000 45.000
92     92 0.000 0.000 45.000
93     93 0.000 0.000 45.000
94     94 0.000 0.000 45.000
95     95 0.000 0.000 45.000
96     96 0.000 0.000 45.000
97     97 0.000 0.000 45.000
98     98 0.000 0.000 45.000
99     99 0.000 0.000 45.000
100   100 0.000 0.000 45.000
> plot(nspx, type="l")
> plot(plex, type="l")
> 

위에서 보듯이 plex column에 해당하는 것이 $E(X^2)$ 부분이고 이는 45가 된다 (45에 수렴한다). 또한 언급한 것처럼 $Var(X) = E(X^2) - E(X)^2$ 이고 $E(X)^2 = 25$ 이므로,
\begin{align*} Var(X) & = E(X^2) - E(X)^2 \\ & = 45 - 25 \\ & = 20 \end{align*}

일반적으로
\begin{align} Var(X) = \displaystyle \frac{q}{p^{2}} \end{align}

아래는 이를 R에서 계산해 본 것

q/(p^2)
[1] 20
> 

\begin{align} P(X = r) & = p * q^{r-1} \\ P(X > r) & = q^{r} \\ P(X \le r) & = 1 - q^{r} \\ E(X) & = \displaystyle \frac{1}{p} \\ Var(X) & = \displaystyle \frac{q}{p^{2}} \\ \end{align}
$(5)$, $(6)$에 대한 증명은 Mean and Variance of Geometric Distribution

$P(X = 2) = p * q^{2-1}$
$P(X \le 4) = 1 - q^{4}$
$P(X > 4) = q^{4}$
$E(X) = \displaystyle \frac{1}{p}$
$Var(X) = \displaystyle \frac{q}{p^{2}}$

1. 1번의 시행에서 특정 사건 A가 발생할 확률을 p라고 하면
2. n번의 (독립적인) 시행에서 사건 A가 발생할 때의 확률 분포를
3. 이항확률분포라고 한다.

아래를 보면

• 각 한문제를 맞힐 확률은 1/4, 틀릴 확률은 3/4
• 3문제를 풀면서 (3번의 시행) 각 문제를 맞힐 확률 분포를 말한다.

$$P(X = r) = {\huge\text{?} \cdot 0.25^{r} \cdot 0.75^{3-r}} $$
$$P(X = r) = {\huge_{3}C_{r}} \cdot 0.25^{r} \cdot 0.75^{3-r}$$

$_{n}C_{r}$은 n개의 사물에서 r개를 (순서없이) 고르는 방법의 수라고 할 때, 3개의 질문 중에서 한 개의 정답을 맞히는 방법은 $_{3}C_{1} = 3$ 세가지가 존재.

Probability for getting one question right
\begin{eqnarray*} P(X = r) & = & _{3}C_{1} \cdot 0.25^{1} \cdot 0.75^{3-1} \\ & = & \frac{3!}{1! \cdot (3-1)!} \cdot 0.25 \cdot 0.75^2 \\ & = & 3 \cdot 0.25 \cdot 0.5625 \\ & = & 3 \cdot 0.25 \cdot 0.5625 \\ & = & 0.421875 \end{eqnarray*}

$$P(X = r) = _{n}C_{r} \cdot 0.25^{r} \cdot 0.75^{n-r}$$
$$P(X = r) = _{n}C_{r} \cdot p^{r} \cdot q^{n-r}$$

1. You’re running a series of independent trials. (n번의 시행을 하게 된다)
2. There can be either a success or failure for each trial, and the probability of success is the same for each trial. (각 시행은 성공/실패로 구분되고 성공의 확률은 (반대로 실패의 확률도) 각 시행마다 동일하다)
3. There are a finite number of trials. Note that this is different from that of geometric distribution. (n번의 시행으로 한정된다. 무한대 시행이 아님)

X가 n번의 시행에서 성공적인 결과를 얻는 수를 나타낸다고 할 때, r번의 성공이 있을 확률을 구하려면 아래 공식을 이용한다.

\begin{eqnarray*} P(X = r) & = & _{n}C_{r} \cdot p^{r} \cdot q^{n-r} \;\;\; \text{Where,} \\ _{n}C_{r} & = & \frac {n!}{r!(n-r)!} \end{eqnarray*}

p = 각 시행에서 성공할 확률
n = 시행 숫자
r = r 개의 정답을 구할 확률

$$X \sim B(n,p)$$

$$X \sim B(1,p)$$

\begin{eqnarray*} E(X) & = & \sum{x * p(x)} \\ & = & (0*q) + (1*p) \\ & = & p \end{eqnarray*}

\begin{eqnarray*} Var(X) & = & E((X - E(X))^{2}) \\ & = & \sum_{x}(x-E(X))^2p(x) \ldots \ldots \ldots E(X) = p \\ & = & (0 - p)^{2}*q + (1 - p)^{2}*p \\ & = & (0^2 - 2p0 + p^2)*q + (1-2p+p^2)*p \\ & = & p^2*(1-p) + (1-2p+p^2)*p \\ & = & p^2 - p^3 + p - 2p^2 + p^3 \\ & = & p - p^2 \\ & = & p(1-p) \\ & = & pq \end{eqnarray*}

For generalization,

$$X \sim B(n,p)$$

\begin{eqnarray*} E(X) & = & E(X_{1}) + E(X_{2}) + ... + E(X_{n}) \\ & = & n * E(X_{i}) \\ & = & n * p \end{eqnarray*}

\begin{eqnarray*} Var(X) & = & Var(X_{1}) + Var(X_{2}) + ... + Var(X_{n}) \\ & = & n * Var(X_{i}) \\ & = & n * p * q \end{eqnarray*}

Ans 1.

p <- .25
q <- 1-p
r <- 2
n <-5
# combinations of 5,2
c <- choose(n,r) 
ans1 <- c*(p^r)*(q^(n-r))
ans1

> p <- .25
> q <- 1-p
> r <- 2
> n <-5
> # combinations of 5,2
> c <- choose(n,r)
> ans <- c*(p^r)*(q^(n-r))
> ans
[1] 0.2636719
> 
> 

Ans 2.

p <- .25
q <- 1-p
r <- 3
n <-5
# combinations of 5,3
c <- choose(n,r)
ans2 <- c*(p^r)*(q^(n-r))
ans2

> p <- .25
> q <- 1-p
> r <- 3
> n <-5
> # combinations of 5,3
> c <- choose(n,r)
> ans2 <- c*(p^r)*(q^(n-r))
> ans2
[1] 0.08789062
> 

Ans 3.

ans1 + ans2 

> ans1 + ans2 
[1] 0.3515625

Ans 4.

p <- .25
q <- 1-p
r <- 0
n <-5
# combinations of 5,3
c <- choose(n,r)
ans4 <- c*(p^r)*(q^(n-r))
ans4

> p <- .25
> q <- 1-p
> r <- 0
> n <-5
> # combinations of 5,3
> c <- choose(n,r)
> ans4 <- c*(p^r)*(q^(n-r))
> ans4
[1] 0.2373047
> 

Ans 5

p <- .25
q <- 1-p
n <- 5
exp.x <- n*p
exp.x

> p <- .25
> q <- 1-p
> n <- 5
> exp.x <- n*p
> exp.x
[1] 1.25

p <- .25
q <- 1-p
n <- 5
var.x <- n*p*q
var.x

> p <- .25
> q <- 1-p
> n <- 5
> var.x <- n*p*q
> var.x
[1] 0.9375
> 

Toss a fair coin once. What is the distribution of the number of heads?

• A single trial
• The trial can be one of two possible outcomes – success and failure
• P(success) = p
• P(failure) = 1-p

X = 0, 1 (failure and success)
$P(X=x) = p^{x}(1-p)^{1-x}$ or
$P(x) = p^{x}(1-p)^{1-x}$

참고.

When x = 0 (failure), $P(X = 0) = p^{0}(1-p)^{1-0} = (1-p)$ = Probability of failure
When x = 1 (success), $P(X = 1) = p^{1}(1-p)^{0} = p $ = Probability of success

$P(X=x) = p^{x}(1-p)^{1-x}$ or
$P(x) = p^{x}(1-p)^{1-x}$
X takes, x = 0, 1

\begin{eqnarray*} E(X) & = & \sum_{x}xP(x) \\ & = & 0*p^{0}(1-p)^{1-0} + 1*p^{1}(1-p)^{1-1} \\ & = & p \\ \\ Var(X) & = & E((X-\mu)^{2}) \\ & = & \sum_{x}(x-\mu)^2P(x) \\ \end{eqnarray*}
그런데
\begin{eqnarray*} E((X-\mu)^{2}) & = & E(X^2) - (E(X))^2 \\ \end{eqnarray*}

위에서
\begin{eqnarray*} E(X^{2}) & = & \sum x^2 p(x) \\ & = & 0^2*p^0(1-p)^{1-0} + 1^2*p^1(1-p)^{1-1} \\ & = & p \end{eqnarray*}

zero squared probability of zero occurring
one squared prob of one occurring

또한 $E(X) = p $ 임을 알고 있음
\begin{eqnarray*} Var(X) & = & E((X-\mu)^{2}) \\ & = & E(X^2) - (E(X))^2 \\ & = & p - p^2 \\ & = & p(1-p) \end{eqnarray*}

위는 First Head Statistics 에서 $X \sim (1, 0.25)$ 에서 E(X)와 Var(X)를 구한 후 (각각, p와 pq), X가 n가지가 있다고 확장하여 np와 npq를 구한 것과 같다. 즉, 교재는 Bernoulli distribution을 이야기(설명)하지 않고, 활용하여 binomial distribution의 기대값과 분산값을 구해낸 것이다.

$E(U_{i}) = p$ and $Var(U_{i}) = p(1-p)$ or $Var(U_{i}) = p \cdot q$

$$X = U_{1} + . . . . + U_{n}$$
\begin{eqnarray*} E(X) & = & E(U_{1} + . . . + U_{n}) \\ & = & E(U_{1}) + . . . + E(U_{n}) \\ & = & p + . . . + p \\ & = & np \end{eqnarray*}

\begin{eqnarray*} Var(X) & = & Var(U_{1} + . . . + U_{n}) \\ & = & Var(U_{1}) + . . . + Var(U_{n}) \\ & = & p(1-p) + . . . + p(1-p) \\ & = & np(1-p) \\ & = & npq \end{eqnarray*}

이항분포에서의 기댓값과 분산에 대한 수학적 증명, Mathematical proof of Binomial Distribution Expected value and Variance

$$X \sim Po(\lambda)$$

단위 시간, 단위 공간에 어떤 사건이 몇 번 발생할 것인가를 표현하는 이산 확률분포
모수(population parameter).

• 단위시간 또는 단위공간에서 평균발생횟수
• lambda (λ)로 표시
  • 한 시간 동안 은행에 다녀간 고객의 수
  • 한 시간 동안 사무실에 걸려온 전화의 수
  • 어떤 책의 한 페이지에 존재하는 오타의 수
  • 팝콘 기계가 일주일 동안 고장나는 횟수

조건

• 개별 사건이 주어진 구간에 임의로 그리고 독립적으로 발생
  • 일주일 동안
  • 1마일마다 등 시간이나 공간
• 해당 구간에서 사건이 발생하는 수의 평균값이나 비율을 알고 있음 (lambda($\lambda$))

$$ P(X=r) = e^{- \lambda} \dfrac{\lambda^{r}} {r!},\qquad k = 0, 1, 2, . . ., $$

For curiosity,
\begin{eqnarray*} \sum_{r=0}^{\infty} e^{- \lambda} \dfrac{\lambda^{r}} {r!} & = & e^{- \lambda} \sum_{r=0}^{\infty} \dfrac{\lambda^{r}} {r!} \\ & = & e^{- \lambda} \left(1 + \lambda + \dfrac{\lambda^{2}}{2!} + \dfrac{\lambda^{3}}{3!} + . . . \right) \\ & = & e^{- \lambda}e^{\lambda} \\ & = & 1 \end{eqnarray*}

왜 $e^{\lambda} = \left(1 + \lambda + \dfrac{\lambda^{2}}{2!} + \dfrac{\lambda^{3}}{3!} + . . . \right)$ 인지는 Taylor series 문서를 참조.
이것이 의미하는 것은 r이 0에서 무한대로 갈 때의 확률값의 분포를 말하므로 전체 분포가 1이 됨을 의미한다. 아래 “What does the Poisson distribution look like?” 참조

> e <- exp(1)
> e
[1] 2.718282

위의 그림은 lambda는 2, 즉 한달에 아주대학교 앞의 건널목 주변 찻길에서 교통사고가 날 횟수가 2회라고 할 때, X=3 이므로 3번 교통사고가 일어날 확률을 (P(X=3)) 묻는 문제이다.
\begin{eqnarray*} P(X = 3) & = & e^{-2} * \frac {2^{3}}{3!} \\ & = & 0.180 \end{eqnarray*}

\begin{eqnarray*} P(X=r) = e^{- \lambda} \dfrac{\lambda^{r}} {r!},\qquad r = 0, 1, 2, . . ., \end{eqnarray*}

마포 신한은행 지점에 시간당 은행에 방문하는 손님의 숫자: lambda = 30

> dpois(x=1:60, lambda=30)
 [1] 2.807287e-12 4.210930e-11 4.210930e-10 3.158198e-09 1.894919e-08
 [6] 9.474593e-08 4.060540e-07 1.522702e-06 5.075675e-06 1.522702e-05
[11] 4.152825e-05 1.038206e-04 2.395861e-04 5.133987e-04 1.026797e-03
[16] 1.925245e-03 3.397491e-03 5.662486e-03 8.940767e-03 1.341115e-02
[21] 1.915879e-02 2.612562e-02 3.407689e-02 4.259611e-02 5.111534e-02
[26] 5.897924e-02 6.553248e-02 7.021338e-02 7.263453e-02 7.263453e-02
[31] 7.029148e-02 6.589826e-02 5.990751e-02 5.285957e-02 4.530820e-02
[36] 3.775683e-02 3.061365e-02 2.416867e-02 1.859128e-02 1.394346e-02
[41] 1.020253e-02 7.287524e-03 5.084319e-03 3.466581e-03 2.311054e-03
[46] 1.507209e-03 9.620485e-04 6.012803e-04 3.681308e-04 2.208785e-04
[51] 1.299285e-04 7.495876e-05 4.242949e-05 2.357194e-05 1.285742e-05
[56] 6.887904e-06 3.625212e-06 1.875110e-06 9.534457e-07 4.767229e-07
> plot(dpois(x=1:60, lambda=30), type = "l")
> 

위에서 언급한

\begin{eqnarray*} \sum_{r=0}^{\infty} e^{- \lambda} \dfrac{\lambda^{r}} {r!} & = & e^{- \lambda} \sum_{r=0}^{\infty} \dfrac{\lambda^{r}} {r!} \\ & = & e^{- \lambda} \left(1 + \lambda + \dfrac{\lambda^{2}}{2!} + \dfrac{\lambda^{3}}{3!} + . . . \right) \\ & = & e^{- \lambda}e^{\lambda} \\ & = & 1 \end{eqnarray*}

에서 1 이란 이야기는 아래 그림의 그래프가 전체가 1이 됨을 의미함. 즉 위에서는 1부터 60까지 갔지만, 1부터 무한대로 하면 완전한 분포곡선이 되는데 이것이 1이라는 뜻 (가령 dpois(x=1:1000, lambda=30)과 같은 케이스).

Figure 1. lambda=30

lambda가 클 수록 좌우대칭의 종형분포를 이루고 ¹⁾, 작을 수록 오른 쪽으로 편향된 (skewed to the right) 혹은 양의방향으로 편향된(positively skewed) 분포를 ²⁾ 이룬다.

> dpois(x=1:60, lambda=.3)
 [1]  2.222455e-01  3.333682e-02  3.333682e-03  2.500261e-04  1.500157e-05
 [6]  7.500784e-07  3.214622e-08  1.205483e-09  4.018277e-11  1.205483e-12
[11]  3.287682e-14  8.219204e-16  1.896739e-17  4.064441e-19  8.128883e-21
[16]  1.524166e-22  2.689704e-24  4.482840e-26  7.078168e-28  1.061725e-29
[21]  1.516750e-31  2.068296e-33  2.697777e-35  3.372222e-37  4.046666e-39
[26]  4.669230e-41  5.188033e-43  5.558607e-45  5.750283e-47  5.750283e-49
[31]  5.564790e-51  5.216991e-53  4.742719e-55  4.184752e-57  3.586930e-59
[36]  2.989108e-61  2.423601e-63  1.913370e-65  1.471823e-67  1.103867e-69
[41]  8.077076e-72  5.769340e-74  4.025121e-76  2.744401e-78  1.829600e-80
[46]  1.193218e-82  7.616283e-85  4.760177e-87  2.914394e-89  1.748636e-91
[51]  1.028610e-93  5.934286e-96  3.359030e-98 1.866128e-100 1.017888e-102
[56] 5.452971e-105 2.869985e-107 1.484475e-109 7.548177e-112 3.774089e-114
> plot(dpois(x=1:60, lambda=.3), type = "l")
> 

Figure 2. lambda = 3

일반적으로 lambda가 1보다 작으면 geometric distribution 형태의 그래프를, 1보다 크면 정규분포 형태의 모양을 갖는다.

1. What’s the probability of the machine not malfunctioning next week?

$\lambda = 3.4$
$\text{malfunctioning} = 0$

\begin{eqnarray*} P(X=0) & = & e^{-3.4} * \frac{3.4^{0}} {0!} \\ & = & e^{-3.4} \\ & = & 0.03337327 \end{eqnarray*}

# R 에서 계산
> e^(-3.4)
[1] 0.03337327
> 
# 혹은 
> dpois(0, 3.4)
[1] 0.03337327
> 

포아송 분포를 따르는 확률에서 아무것도 일어나지 않을 때의 확률은 e^-lambda 가 된다. 예를 들면 119 전화가 한시간에 걸려오는 확률이 5번이라고 할 때, 지난 한 시간동안 한 건의 전화도 없을 확률은?
\begin{eqnarray*} P(X=0) & = & e^{-5} * \frac{5^{0}} {0!} \\ & = & e^{-5} \\ & = & 0.006737947 \end{eqnarray*}

> lamba <- 5
> e <- exp(1)
> px.0 <- e^(-lamba)
> 
> px.0
[1] 0.006737947
> 
# or 
> dpois(0,5)
[1] 0.006737947

2. What’s the probability of the machine malfunctioning three times next week?

l <- 3.4
x <- 3
e <- exp(1)
ans <- ((e^(-l))*l^x)/factorial(x)

> l <- 3.4
> x <- 3
> e <- exp(1)
> ans <- ((e^(-l))*l^x)/factorial(x)
>     
> ans
[1] 0.2186172
> 

위의 계산 대신 아래와 같은 function을 이용하는 것이 보통이다.

> dpois(x=3, lambda=3.4)
[1] 0.2186172

마찬가지로 적어도 3번까지 고장나는 경우는 0, 1, 2, 3을 포함하므로

> sum(dpois(c(0:3), lambda=3.4))
[1] 0.5583571
> 

3. What’s the expectation and variance of the machine malfunctions?
\begin{eqnarray*} E(X) & = & \lambda \\ Var(X) & = & \lambda \\ & = & 3.4 \end{eqnarray*}

mean and variance of Poisson distribution

\begin{eqnarray*} X \sim Po(3.4) \\ Y \sim Po(2.3) \end{eqnarray*}

위의 조건일 때, Popcorn 기계와 coffee 기계가 한 주일 동안 고장나지 않을 확률을 구하려면 아래를 말한다.
\begin{eqnarray*} P(X + Y = 0) \end{eqnarray*}

여기서 X + Y의 분포는 아래와 같다.
\begin{eqnarray*} X + Y \sim (\lambda_{x} + \lambda_{y}) \end{eqnarray*}

lambda의 합은 5.7이고 (아래 참조), 결국 lambda가 5.7일 때 X=0의 확률(probability)를 구하는 문제이므로 0.003
\begin{eqnarray*} \lambda_{X} + \lambda_{Y} & = & 3.4 + 2.3 \\ & = & 5.7 \\ \end{eqnarray*}

$$X + Y \sim Po(5.7)$$

\begin{eqnarray*} P(X + Y = 0) & = & \frac {e^{- \lambda} \lambda^{r}} {r!} \\ & = & \frac {e^{-5.7} 5.7^{0}}{0!} \\ & = & e^{-5.7} \\ & = & 0.003 \end{eqnarray*}

우선 위의 문제를 binomial distribution 문제로 생각하면 답은
\begin{eqnarray*} P(r=15) & = & _{100}C_{15} * 0.1^{15} * 0.99^{85}\\ \end{eqnarray*}
라고 볼 수 있다.

\begin{eqnarray} X & \sim & B(n, p) \\ X & \sim & Po(\lambda) \end{eqnarray}

Poisson distribution을 대신 사용할 수 있으려면, B(n, p)와 Po(lambda)가 유사해야 한다. 두 distribution의 기대값과 분산값을 살펴보면,

• B(n, p)의 경우 E(X) = np
• Po(lambda)의 경우 E(X) = lambda 이고

• Var(X) = npq 이고
• Var(lambda) = lambda 이다.

따라서, 둘의 성격이 비슷하기 위해서는 npq 와 np가 같아야 한다. 따라서 q는 1이어야 하는데, 현실적으로 1일 수는 없으므로 1에 가깞고, n이 충분히 크다면 둘의 성격이 비슷해질 수 있다고 판단한다. 따라서,

• 만약 n이 충분히 크고
• p가 작으면 (q가 크면)
• $X \sim B(n, p)$와 $Y \sim Po(np)$는 비슷할 것이다.
• 보통은 n > 50인 경우, p = 0.1 보다 작은 경우가 위에 해당한다.

> dbinom(x=15, 100, 0.1)
[1] 0.03268244
> choose(100, 15)
[1] 2.533385e+17
> a <- choose(100, 15)
> b <- .1^15
> c <- .9^85
> a*b*c
[1] 0.03268244
> 

위가 답이긴 하지만 limited calculator 로는
x ~ b (n, p)이고
b(100, 0.1)이므로
n*p = 10 = lambda
따라서 Pois 분포로 보는 답은
lambda = 10 일때 P(r=15)값을 구하는 문제로

\begin{eqnarray*} P(r = 15) & = & e^{-10} * \frac {10^{15}}{15!} \\ & = & 0.0347180 \end{eqnarray*}

> dpois(x=15, lambda=10)
[1] 0.03471807
> 

> dbinom(x=5, 50, 0.05)
[1] 0.06584064
> 

> dpois(x=5, lambda = 2.5)
[1] 0.06680094
> 

> n <- 50
> p <- .05
> q <- 1-p
> x <- 5
> np <- n*p
# Poisson distribution
> e <- exp(1)
> lambda <- np
> lambda
[1] 2.5
> a <- e^(-lambda)
> b <- lambda^x
> c <- factorial(x)
> a*b/c
[1] 0.06680094
>