set.seed(0)


## 3.1:
##
n_tosses = 100
tosses = sample( c(0, 1), n_tosses, replace=TRUE ) ## 0=> tails; 1 => heads
n_heads = cumsum( tosses )
p_head = n_heads/(1:n_tosses)
plot( 1:n_tosses, p_head, type='p', xlab='number of tosses', ylab='p_head (empirical)' )
abline(h=0.5, col='green')
grid()

p_0 = 0.5
sigma_d = sqrt(n_tosses*p_0*(1-p_0))
print(sigma_d)


## 3.2:
##
n_tosses = 100
tosses = sample( seq( 1, 6 ), n_tosses, replace=TRUE )
n_sixes = cumsum(tosses==6)
p_six = n_sixes/(1:n_tosses)
plot(1:n_tosses, p_six, type='p', xlab='number of tosses', ylab='p_six(empirical)')
abline(h=1/6, col='green')
grid()

p_0 = 1/6
sigma_d = sqrt(n_tosses*p_0*(1-p_0))
print(sigma_d)


## 3.3:
##
print(dbinom( 0:2, 2, 0.02 ))

## 3.4:
##
print(dbinom( 0:2, 2, 0.01 ))


## 3.5:
##
print(dbinom( 0:3, 3, 0.08 ))


## 3.6:
##
print(dhyper(0:1, 1, 7, 2))


## 3.7:
##
print(dhyper(0:1, 1, 7, 3))


## 3.8:
##
print(dhyper(0:2, 2, 6, 2))


## 3.9:
##
print(dhyper(0:2, 2, 6, 3))


## 3.10:
##
n = 400
p_0 = 0.1
m = n*p_0
v = n*p_0*(1-p_0)
sigma = sqrt(v)
print(sprintf('mean= %f; sigma= %f; 3*sigma= %f', m, sigma, 3*sigma) )


## 3.11:
##
Cs = choose(10, c(4, 6))
Ps = factorial(10)/factorial(10-4)
print(c(Cs, Ps))


## 3.12:
##
denom = choose(52, 4)
probs = c(13, 13*4*48)/denom
print(probs)


## 3.13:
##
p_a = ppois(1, 1.2)
p_b = 1-ppois(2, 1.2)
print(c(p_a, p_b))


## 3.14:
##
p_a = ppois(3, 3.6)
p_b = 1-ppois(4, 3.6)
print(c(p_a, p_b))


## 3.15:
##
p_a = ppois(2, 2.0)
p_b = dpois(2, 2.0)
p_c = 1-ppois(1, 2.0)
print(c(p_a, p_b, p_c))


## 3.16:
##
p_a = ppois(2, 5.4)
p_b = 1-ppois(7, 5.4)
print(c(p_a, p_b))


## 3.17:
##
p_a = dhyper(0, 100, 1000-100, 5)
p_c = dbinom(0, 5, 0.1)
print(c(p_a, p_c))


## 3.18:
##
p_a = dbinom(0, 10, 0.01)
p_c = dpois(0, 0.1)
print(c(p_a, p_c))


## 3.19:
##
p_a = dbinom(0:2, 100, 0.002)
print(p_a)
p_c = dpois(0:2, 100*0.002)
print(p_c)


## 3.20 - 3.22:
##
c0s = c(3.6, 2.0, 5.4)
sigma_c = sqrt(c(3.6, 2.0, 5.4))
ll = ceiling(c0s + 3*sigma_c) ## the 3 sigma limits
print(ll)
print(1-ppois(ll-1, c0s))


## 3.23:
##
K = 100
N = 1000
n = 5
m = n * (K/N) ## expressions for the mean of hypergeometric RV
v = n * (K/N) * ((N-K)/N) * ((N-n)/(N-1)) ## expressions for the variance of hypergeometric RV
s = sqrt(v)
ll = ceiling(m + 3*s)
print(1-phyper(ll-1, 100, 1000-100, 5))



## 3.24:
##
print(dbinom(0:4, 4, 0.2))


## 3.25:
##
print(dbinom(0:5, 5, 0.1))


## 3.26:
##
print(dhyper(0:2, 2, 7-2, 2))


## 3.27:
##
print(dhyper(0:3, 3, 9-3, 3))