COMMunication
RESearch.NET

This page is read only. You can view the source, but not change it. Ask your administrator if you think this is wrong.
====== Geometric Binomial and Poisson Distributions ======
정리 
\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*}

===== Geometric Distributions =====

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

<WRAP box help>
It’s time to exercise your probability skills. The probability of Chad making a successful run down the slopes is 0.2 for any given trial (assume trials are independent). What’s the probability he’ll need two trials? What’s the probability he’ll make a successful run down the slope in one or two trials? Remember, when he’s had his first successful run, he’s going to stop.

Hint: You may want to draw a probability tree to help visualize the problem.
</WRAP>

{{:b:head_first_statistics:pasted:20191030-004453.png?400}}

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

{{:b:head_first_statistics:pasted:20191030-005359.png}}

| X  | P(X=x)  |
| 1  | 0.2  |
| 2  | 0.8 * 0.2 = 0.16  |
| 3  | 0.8 * 0.8 * 0.2 = 0.128 |
| 4  | 0.8 * 0.8 * 0.8 * 0.2 =  0.1024  |
| . . .  | . . . . .  |

| X  | P(X=x)  | Power of 0.8  | Power of 0.2  |
| 1  | 0.2  | 0  | 1  |
| 2  | 0.8 * 0.2   | 1  | 1  |
| 3  | 0.8<sup>2</sup> * 0.2   | 2  | 1  |
| 4  | 0.8<sup>3</sup> * 0.2   | 3  | 1  |
| 5  | 0.8<sup>4</sup> * 0.2   | 4  | 1  |
| r  | . . . . .  | r - 1  | 1  |

$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.
  - 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}} $

<code>
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))
</code>

<code>
> 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))
</code>

{{:b:head_first_statistics:pasted:20191030-023820.png}}

r번 시도한 다음에 성공을 얻을 확률
첫 번째 성공을 얻을 때까지 r번 이상 시도를 해야하는 확률
$$ P(X > r) = q^{r} $$

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

Solution. 
  * 21번째 성공 + 22번째 + 23번째 + . . . . 
  * 위는 구할 수 없음
  * 따라서
  * 1 - (1번째 성공 + 2번째 성공 + . . . + 20번째 성공)
  * 그런데 이것은 
  * 20번까지는 실패하는 확률 = $q^{r} $ 과 같다
<code>
p <- .2
q <- 1-p
n <- 19
s <- dgeom(x = 0:n, prob = p)
# 20번째까지 성공할 확률
sum(s)
# 따라서 아래는 20번 이후에 성공할 확률
1-sum(s)
## 혹은 (교재가 이야기하는) 20번까지 실패하는 확률
q^20
</code>

<code>
> 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
> 
</code>
{{:b:head_first_statistics:pasted:20191031-004954.png}}

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

혹은 1번째 성공 + 2번째 성공 + . . . + r 번째 성공으로 구해도 된다
<code>
# r = 20 이라고 하면 
p <- .2
q <- 1-p
n <- 19
s <- dgeom(x = 0:n, prob = p)
sum(s)
</code>


{{:b:head_first_statistics:pasted:20191031-005143.png}}

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

{{:b:head_first_statistics:pasted:20191031-005411.png}}

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

Reminding . . . [[:b/head_first_statistics/using_discrete_probability_distributions#expected_value_discrete_prob|Expected value in discrete probability]]
$E(X) = \sum x*P(X=x)$

| textbook  | x  | P(X = x)  | xP(X = x)  | xP(X ≤ x): \\  $E(X) = \sum (x*P(X=x))$  |
| r code  | trial  | ''px <- q^(trial-1)*p''  | ''npx <- trial*(q^(trial-1))*p''  |  ''plex <- cumsum(trial*(q^(trial-1))*p)''  |
|   |   | ''px''   | ''npx <- trial*px''   | ''plex <- cumsum(npx)''   |
|   |   | x번째 (trial번째) \\ 성공할 확률  | x번째의 기대치 \\ (주사위 경우처럼)  | 그 x번째까지 성공할 \\ 확률에 대한 기대값  |



<code>
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)
</code>

<code>
> 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
> 
</code>

<code>
> 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")
</code>
{{:b:head_first_statistics:pasted:20191031-010726.png}}
{{:b:head_first_statistics:pasted:20191031-010811.png}}
위에서 $X \sim \text{Geo}(p)$ 일때, 기대값은 $E(X) = \displaystyle \frac{1}{p}$ 임을 알 수 있다. 아래는 그 증명이다.

===== Proof =====
[[:Mean and Variance of Geometric Distribution#Mean|기하분포의 기대값(expected value)]]


===== Variance proof =====
[[:Mean and Variance of Geometric Distribution#Variance|기하분포의 분산 증명]]
아래는 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*}

<code>
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")
</code>

<code>
> 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")
> 
</code>
{{:b:head_first_statistics:pasted:20191031-062004.png}}
{{:b:head_first_statistics:pasted:20191031-062021.png}}

위에서 보듯이 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에서 계산해 본 것
<code>
q/(p^2)
[1] 20
> 
</code>
===== Sum up =====
\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]]
===== e.g., =====
<WRAP box>
The probability that another snowboarder will make it down the slope without falling over is 0.4. Your job is to play like you’re the snowboarder and work out the following probabilities for your slope success.

  - The probability that you will be successful on your second attempt, while failing on your first.
  - The probability that you will be successful in 4 attempts or fewer.
  - The probability that you will need more than 4 attempts to be successful.
  - The number of attempts you expect you’ll need to make before being successful.
  - The variance of the number of attempts.
</WRAP>

$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}}$





====== Binomial Distributions ======

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

{{:b:head_first_statistics:pasted:20191030-035316.png}}
{{:b:head_first_statistics:pasted:20191030-035452.png}}

| x  | P(X=x)                    | power of .75  | power of .25  |
| 0  | 0.75 * 0.75 * 0.75        | 3  | 0  |
| 1  | 3 * (0.75 * 0.75 * 0.25)  | 2  | 1  |
| 2  | 3 * (0.75 * 0.25 * 0.25)  | 1  | 2  |
| 3  | 0.25 * 0.25 * 0.25        | 0  | 3  |
{{:b:head_first_statistics:pasted:20191030-040346.png}}

$$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}$$

  - You’re running 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.
  - There are a finite number of trials. (note that this is different from that of geometric distribution)

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)$$

===== Expectation and Variance of =====

{{:b:head_first_statistics:pasted:20191104-013651.png}}

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

\begin{eqnarray*}
E(X) & = & \sum{n*p(x)} \\
& = & (1*p)+(0*q) \\
& = & 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*}

===== e.g., =====
<WRAP box>
In the latest round of Who Wants To Win A Swivel Chair, there are 5 questions. The probability of
getting a successful outcome in a single trial is 0.25
  - What’s the probability of getting exactly two questions right?
  - What’s the probability of getting exactly three questions right? 
  - What’s the probability of getting two or three questions right? 
  - What’s the probability of getting no questions right?
  - What are the expectation and variance?
</WRAP>

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

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

Ans 3. 
<code>
ans1 + ans2 
</code>

<code>> ans1 + ans2 
[1] 0.3515625
</code>

Ans 4. 
<code>
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
</code>

<code>> 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
> </code>

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

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

===== Another way to see E(X) and Var(X) =====
==== Bernoulli Distribution ====
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}$

참고.
| x  | 0  | 1  |
| p(x)  | q = (1-p)  | p  | 

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

<WRAP box>
Bernoulli distribution expands to binomial distribution, geometric distribution, etc.
Binomial distribution = The distribution of number of success in n independent Bernoulli trials.
Geometric distribution = The distribution of number of trials to get the first success in independent Bernoulli trials.
</WRAP>

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

==== Expectation and Variance value ====
\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의 기대값과 분산값을 구해낸 것이다. 

==== Proof of E and Var from Bernoulli 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*}


==== From a scratch (Proof of Binomial Expected Value) ====
[[:Mean and Variance of Binomial Distribution|이항분포에서의 기댓값과 분산에 대한 수학적 증명]], Mathematical proof of Binomial Distribution Expected value and Variance
====== Poisson Distribution ======
$$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]] 문서를 참조.

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

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

===== What does the Poisson distribution look like? =====

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

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

<code>
> 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")
> </code>

위에서 언급한 

\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)과 같은 케이스). 


[{{:b:head_first_statistics:pasted:20191107-095627.png|Figure 1. lambda=30}}]

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

<code>
> 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")
> </code>

[{{:b:head_first_statistics:pasted:20191107-095447.png|Figure 2. lambda = 3}}]

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

<WRAP box>
Your job is to play like you’re the popcorn machine and say what the probability is of you malfunctioning a  particular number of times next week. Remember, the mean number of times you break down in a week is 3.4.

  - What’s the probability of the machine not malfunctioning next week?
  - What’s the probability of the machine malfunctioning three times next week?
  - What’s the expectation and variance of the machine malfunctions?
</WRAP>

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

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

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

<code>
> e^(-3.4)
[1] 0.03337327
> 
</code>

__2. What’s the probability of the machine malfunctioning three times next week?__
<code>
l <- 3.4
x <- 3
e <- exp(1)
ans <- ((e^(-l))*l^x)/factorial(x)
</code>

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

위의 계산 대신 아래와 같은 function을 이용하는 것이 보통이다.
<code>
> dpois(x=3, lambda=3.4)
[1] 0.2186172
</code>

마찬가지로 적어도 3번까지 고장나는 경우는 0, 1, 2, 3을 포함하므로 
<code>
> sum(dpois(c(0:3), lambda=3.4))
[1] 0.5583571
> 
</code>

__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]]
===== Two Poisson distribution cases =====

\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*}

===== Broken Cookies case =====
<WRAP box>
The Case of the Broken Cookies

Kate works at the Statsville cookie factory, and her job is to make sure that boxes of cookies meet the factory’s strict rules on quality control. Kate know that **the probability that a cookie is broken is 0.1**, and her boss has asked her to find the probability that there will be 15 broken cookies in a box of 100 cookies. "It’s easy," he says. "Just use the binomial distribution where n is 100, and p is 0.1."

Kate picks up her calculator, but when she tries to calculate 100!, her calculator displays an error because the number is too big. “Well,” says her boss, "you’ll just have to calculate it manually. But I’m going home now, so have a nice night."

Kate stares at her calculator, wondering what to do. Then she smiles. "Maybe I can leave early tonight, after all."
Within a minute, Kate’s calculated the probability. She’s managed to find the probability and has managed to avoid calculating 100! altogether. She picks up her coat and walks out the door.

**How did Kate find the probability so quickly, and avoid the error on her calculator?**
</WRAP>

\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 보다 작은 경우가 위에 해당한다.


<code>
> 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
> 
</code>
위가 답이긴 하지만 limited calculator 로는 
x ~ b (n, p)이고
b(100, 0.1)이므로 
n*p = 10 = lambda 
따라서
<code>
> dpois(x=15, lambda=10)
[1] 0.03471807
> 
</code>


<WRAP box>
A student needs to take an exam, but hasn’t done any revision for it. He needs to guess the answer to each question, and the probability of getting a question right is 0.05. There are 50 questions on the exam paper. What’s  the probability he’ll get 5 questions right? Use the Poisson approximation to the binomial distribution to find out.

만약에 binomial distribution 으로 계산을 한다면
<code>
> dbinom(x=5, 50, 0.05)
[1] 0.06584064
> 
</code>

Poisson distribution을 이용하라고 한다. . . 

$ X \sim B(50, 0.05) $ 일 때, $P(X=5)$를 구하는 것. 이 때의 기대값 E(X)는 $ E(X) = np = 50 * .05 = 2.5 $ 이므로 위의 문제는 

\begin{eqnarray*}
X & \sim & Po(\lambda) \\
X & \sim & Po(2.5)
\end{eqnarray*}

일 때, $P(X=5)$를 구하는 것과 같다.

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

수식을 따르면, 
\begin{eqnarray*}
P(X = 5) & = & \frac {e^{-2.5} * 2.5^{5}}{5!} \\
& = & 0.067
\end{eqnarray*}

<code>
> 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
> 
</code>

</WRAP>

===== Exercise =====
<WRAP box>
Here are some scenarios. Your job is to say which distribution each of them follows, say what the expectation and variance are, and find any required probabilities. 

1. A man is bowling. The probability of him knocking all the pins over is 0.3. If he has 10 shots, what’s the probability he’ll knock all the pins over less than three times?

Binomial distribution 을 이용한다면, 
\begin{eqnarray*}
X & \sim & B(n, p) \\
X & \sim & B(10, 0.3) 
\end{eqnarray*}

\begin{eqnarray*}
E(X) & = & np \\
& = & 10 * 0.3 \\
& = & 3
\end{eqnarray*}

\begin{eqnarray*}
Var(X) & = & npq \\
& = & 10 * 0.3 * 0.7 \\
& = & 2.1
\end{eqnarray*}

r을 이용한다면 ''pbinom'' 을 이용한다.

<code>
> pbinom(q=2, 10, 0.3)
[1] 0.3827828
> </code>

손으로 계산을 한다고 하면, 
$P(X=0), P(X=1), P(X=2)$를 구한 후 모두 더하여 P(X < 3)을 구한다. 

\begin{eqnarray*}
P(X = 0) & = & {10 \choose 0} * 0.3^0 * 0.7^10 \\
& = & 1 * 1 * 0.028 \\
& = & 0.028 
\end{eqnarray*}

\begin{eqnarray*}
P(X = 1) & = & {10 \choose 1} *0.3^1 * 0.7^9 \\
& = & 10 * 0.3 * 0.04035 \\
& = & 0.121
\end{eqnarray*}

\begin{eqnarray*}
P(X = 2) & = & {10 \choose 2} * 0.3^2 * 0.7^8 \\
& = & 45 * 0.09 * 0.0576 \\
& = & 0.233
\end{eqnarray*}

\begin{eqnarray*}
P(X<3) & = & P(X=0) + P(X=1) + P(X=2) \\
& = & 0.028 + 0.121 + 0.233 \\
& = & 0.382
\end{eqnarray*}
</WRAP>

<WRAP box>
2. On average, 1 bus stops at a certain point every 15 minutes. What’s the probability that no buses will turn up in a single 15 minute interval?

위는 Poisson distribution 문제이므로 기대값과 분산값은 각각 lambda 값인 1 (15분마다 1대씩 버스가 온다고 한다)

\begin{eqnarray*}
P(X=0) & = & \frac {e^{-1}{1^0}}{0!} \\
& = & \frac {e^-1 * 1}{1} \\
& = & .368 
\end{eqnarray*}
</WRAP>

<WRAP box>
3. 20% of cereal packets contain a free toy. What’s the probability you’ll need to open fewer than 4 cereal packets before finding your first toy?

이는 geometric distribution 문제이므로, 

$$X \sim Geo(.2)$$

$P(X \le 3)$ 을 구하는 문제이므로 

\begin{eqnarray*}
P(X \le 3) & = & 1 - q^r \\
& = & 1 - 0.8^{3} \\
& = & 1 - 0.512 \\
& = & 0.488
\end{eqnarray*}

기대값과 분산은 각각 $1/p$, $q/p^2$ 이므로 $5$와 $20$.
</WRAP>