Assignment2

library(discreteRV)

## 
## Attaching package: 'discreteRV'

## The following object is masked from 'package:base':
## 
##     %in%

library(data.table)

Assignment 2: Jonatan Sanchez, Fernando Garcia, Raul Garcia

Binomial distribution

Exercises done by hand

2- 0,314

4- Expected = 1,376, Variance = 1,25, Standard deviation = 1,12

5- 0,237

6- 0,0014

8- 0,763

Done in R

#1
experiments= 1e5;
n = 9;
p = 0.649;

##simulation
table(rbinom(experiments, n, p))

## 
##     0     1     2     3     4     5     6     7     8     9 
##     4   145   988  4293 12018 22137 26910 21443  9999  2063

##estimated PMF
xd = 1;
pmf = c(dbinom(0, n, p))
while(xd <= n){
  pmf = c(pmf, dbinom(xd, n, p))
  xd = xd + 1;
}
plot(x= 0:9 ,y = pmf, type="l", xlab = "x", ylab = "P(x)", main = "PMF")

##estimated CDF
xd = 1;
cdf = c(pbinom(0, n, p))
while(xd <= n){
  cdf = c(cdf, pbinom(xd, n, p))
  xd = xd + 1;
}
cumsum = cumsum(pmf);
plot(c(0,0:9), c(0,cumsum), type="s", xlab = "x", ylab = "F(x)", main = "CDF")

#2
dbinom(7, size = 9, prob= 0.649) + dbinom(8, size = 9, prob=0.649) 

## [1] 0.3145223

#3
experiments2 = 20000;
nF = 16;
pC = 1 - 0.914;
##simulation
table(rbinom(experiments2, nF, pC))

## 
##    0    1    2    3    4    5    6    7 
## 4703 7152 5054 2173  714  175   24    5

##estimated PMF
xd = 1;
pmf = c(dbinom(0, nF, pC))
while(xd <= nF){
  pmf = c(pmf, dbinom(xd, nF, pC))
  xd = xd + 1;
}
plot(x= 0:16 ,y = pmf, type="l", xlab = "x", ylab = "P(x)", main = "PMF")

##estimated CDF
xd = 1;
cdf = c(pbinom(0, nF, pC))
while(xd <= nF){
  cdf = c(cdf, pbinom(xd, nF, pC))
  xd = xd + 1;
}
cumsum = cumsum(pmf);
plot(c(0,0:16), c(0,cumsum), type="s", xlab = "x", ylab = "F(x)", main = "CDF")

#4
Ex = nF*pC
Ex2 = mean(rbinom(experiments, nF, pC))
Vx = Ex*(1-pC)
Vx2 = var(rbinom(experiments, nF, pC))
Devx = sqrt(Vx)
Devx2 = sd(rbinom(experiments, nF, pC))

Ex

## [1] 1.376

Ex2

## [1] 1.37428

Vx

## [1] 1.257664

Vx2

## [1] 1.254637

Devx

## [1] 1.121456

Devx2

## [1] 1.124302

#5
dbinom(0, size = nF, prob = pC)

## [1] 0.2372134

#6
1 - (dbinom(0, size = nF, prob = pC) + dbinom(1, size = nF, prob = pC) + dbinom(2, size = nF, prob = pC) + dbinom(3, size = nF, prob = pC) + dbinom(4, size = nF, prob = pC) + dbinom(5, size = nF, prob = pC))

## [1] 0.001515806

#7
xd = 1;
pmf = c(dbinom(0, nF, pC))
while(xd <= nF){
  pmf = c(pmf, dbinom(xd, nF, pC))
  xd = xd + 1;
}
plot(x= 0:16 ,y = pmf, type="l", xlab = "x", ylab = "P(x)", main = "PMF")

#8
1 - dbinom(0, size = nF, prob = pC)

## [1] 0.7627866

Geometric distribution

Javi <- 0.413 
Python <- 0.649
n <- 200000
table <- data.table(PMF=dgeom(rgeom(n,Javi),Javi),CDF=pgeom(rgeom(n,Javi),Javi))
table # 1

##                PMF       CDF
##      1: 0.24243100 0.4130000
##      2: 0.24243100 0.6554310
##      3: 0.41300000 0.4130000
##      4: 0.14230700 0.4130000
##      5: 0.14230700 0.8812722
##     ---                     
## 199996: 0.08353421 0.7977380
## 199997: 0.24243100 0.4130000
## 199998: 0.41300000 0.4130000
## 199999: 0.24243100 0.6554310
## 200000: 0.24243100 0.4130000

table <- data.table(PMF=dgeom(rgeom(n,Python),Python),CDF=pgeom(rgeom(n,Python),Python)) # Using n/2 because it's 100k
table # 4

##                PMF       CDF
##      1: 0.07995745 0.6490000
##      2: 0.22779900 0.6490000
##      3: 0.22779900 0.8767990
##      4: 0.22779900 0.6490000
##      5: 0.22779900 0.6490000
##     ---                     
## 199996: 0.64900000 0.6490000
## 199997: 0.22779900 0.9567564
## 199998: 0.22779900 0.6490000
## 199999: 0.22779900 0.9981300
## 200000: 0.22779900 0.8767990

table <- data.table(PMF=dgeom(rgeom(n,Python),Python),CDF=pgeom(rgeom(n,Python),Python)) # 7
table # 7

##                 PMF   CDF
##      1: 0.079957449 0.649
##      2: 0.227799000 0.649
##      3: 0.009850838 0.649
##      4: 0.009850838 0.649
##      5: 0.227799000 0.649
##     ---                  
## 199996: 0.227799000 0.649
## 199997: 0.649000000 0.649
## 199998: 0.649000000 0.649
## 199999: 0.079957449 0.649
## 200000: 0.227799000 0.649

#3
pmf = dgeom(0:10, Python)
plot(pmf, type="l", xlab = "x", ylab = "P(x)", main = "PMF")

#2
cumsum = cumsum(pmf)
plot(c(0,0:10*2), c(0,cumsum), type="s", xlab = "x", ylab = "F(x)", main = "CDF")

#5
pmf = dgeom(0:20, (1-Python))
cumsum = cumsum(pmf)
plot(c(0,0:20*2), c(0,cumsum), type="s", xlab = "x", ylab = "F(x)", main = "CDF")

#6
pmf = dgeom(0:10, (1-Javi))
plot(pmf, type="l", xlab = "x", ylab = "P(x)", main = "PMF")

Hypergeometric distribution

Exercises done by hand

1- It will always be 0 probability (0%) because if there are 82 and 75 of them are engineers, there will always be more than between 68 and 74 engineers (atleast before the decision of who goes to the BOD).

4- 0,575

5- 0,00000449

7- 0,000004518

8- 0,92

Using R

dhyper(7, 7, 75, 6) # Exercise 1 

## [1] 0

plot(0:6, dhyper(0:6, 7, 75, 6), type="h", main="PMF economists on BOD")# Exercise 2

plot(0:6, phyper(0:6, 7, 75, 6, ),type="s",main="CDF economists on BOD")# Exercise 2

plot(0:75, dhyper(0:75, 7, 75, 6, ),type="h",main="PMF engineers on BOD")# Exercise 3

plot(0:75, phyper(0:75, 7, 75, 6, ),type="s",main="CDF engineers on BOD")# Exercise 3

dhyper(0, 7, 75, 6)# Exercise 4

## [1] 0.5750471

dhyper(5, 7, 75, 6) # Exercise 5 

## [1] 4.497921e-06

plot(0:6, phyper(0:6, 7, 75, 6, ),type="s",main="CDF economists on BOD")# Exercise 6

dhyper(5, 7, 75, 6) # Exercise 7 

## [1] 4.497921e-06

1-dhyper(6, 7, 75, 6) # Exercise 8

## [1] 1

Joint distribution

Exercise results by hand

1- Marginal distribution X1 = 0.34

2- Marginal distribution Y1 = 0.26

3- Expected value of X1 = 0.49

4- Expected value of Y1 = 0.19

5- Variance of X1 = 2.31

6- Variances of Y1 = 0.85

Done in R

pXY=c(.12 , .09, .05 , .12 , .09, .05 ,.12, .08, .05, .11, .08, .05) 
#the vector of joint probabilities from the rat-drug example
pXYm =matrix(pXY,3,4) # the matrix of joint probabilities
px=rowSums(pXYm) 
py=colSums(pXYm) 
c(.12 , .09, .05 , .12 , .09, .05 ,.12, .08, .05, .11, .08, .05) #the vector of joint probabilities from the rat-drug example

##  [1] 0.12 0.09 0.05 0.12 0.09 0.05 0.12 0.08 0.05 0.11 0.08 0.05

matrix(pXY,3,4) # the matrix of joint probabilities

##      [,1] [,2] [,3] [,4]
## [1,] 0.12 0.12 0.12 0.11
## [2,] 0.09 0.09 0.08 0.08
## [3,] 0.05 0.05 0.05 0.05

rowSums(pXYm) # the marginal PMFs

## [1] 0.47 0.34 0.20

colSums(pXYm) # the marginal PMFs

## [1] 0.26 0.26 0.25 0.24

# so X1 = 0.34 and Y1 = 0.26
Eu <- sum(px * pXYm)
Ev <- sum(py * pXYm)
Euv <- sum(px * py * pXYm)

## Warning in px * py: longitud de objeto mayor no es múltiplo de la longitud
## de uno menor

sum(px * pXYm)

## [1] 0.3765

sum(py * pXYm)

## [1] 0.2551

sum(px * py * pXYm)

## Warning in px * py: longitud de objeto mayor no es múltiplo de la longitud
## de uno menor

## [1] 0.094572

Euv - Eu * Ev

## [1] -0.00147315

Poisson distribution

Done in R

#1
ppois(26, lambda = 71, lower.tail = FALSE)

## [1] 1

#2
pmf = dpois(0:60*2, lambda = 71)
cumsum = cumsum(pmf)
plot(c(0,0:60*2), c(0,cumsum), type="s", xlab = "x", ylab = "F(x)", main = "CDF")

#3
ppois(55, lambda = 71, lower.tail = FALSE)

## [1] 0.9708053

#4
pmf = dpois(0:1440*2, lambda = 71*24)
cumsum = cumsum(pmf)
plot(c(0,0:1440*2), c(0,cumsum), type="s", xlab = "x", ylab = "F(x)", main = "CDF")

#5
sum(dpois(18:70, lambda = 71))

## [1] 0.4842169

#6
plot(dpois(0:60*2, lambda = 71), type="l", xlab = "x", ylab = "P(x)", main = "PMF")

#7
##Simulation
table(rpois(experiments2, lambda = 71))

## 
##  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58 
##   2   4   3   1   6   5  13  15  31  36  46  73  87 119 149 191 218 300 
##  59  60  61  62  63  64  65  66  67  68  69  70  71  72  73  74  75  76 
## 354 412 462 554 625 692 779 842 839 897 910 985 958 993 882 860 845 795 
##  77  78  79  80  81  82  83  84  85  86  87  88  89  90  91  92  93  94 
## 721 598 534 512 448 399 320 304 248 193 163 146 103  75  70  45  43  25 
##  95  96  97  98  99 100 101 102 103 105 108 110 
##  20  10  13  10   5   3   1   4   1   1   1   1

data.table(PMF = dpois(0:60*2, lambda = 71), CDF = ppois(0:60*2, lambda = 71))

##              PMF          CDF
##  1: 1.462486e-31 1.462486e-31
##  2: 3.686197e-28 3.791496e-28
##  3: 1.548510e-25 1.639541e-25
##  4: 2.602013e-23 2.838296e-23
##  5: 2.342276e-21 2.634577e-21
##  6: 1.311935e-19 1.523060e-19
##  7: 5.010199e-18 6.009300e-18
##  8: 1.387715e-16 1.721442e-16
##  9: 2.914780e-15 3.743776e-15
## 10: 4.801767e-14 6.393493e-14
## 11: 6.369922e-13 8.803616e-13
## 12: 6.950385e-12 9.984387e-12
## 13: 6.347263e-11 9.491255e-11
## 14: 4.922546e-10 7.674295e-10
## 15: 3.282349e-09 5.344226e-09
## 16: 1.901876e-08 3.239908e-08
## 17: 9.664674e-08 1.726049e-07
## 18: 4.342212e-07 8.147630e-07
## 19: 1.737229e-06 3.432841e-06
## 20: 6.228573e-06 1.299502e-05
## 21: 2.012708e-05 4.446129e-05
## 22: 5.892020e-05 1.382357e-04
## 23: 1.569856e-04 3.925081e-04
## 24: 3.823016e-04 1.022498e-03
## 25: 8.542475e-04 2.454265e-03
## 26: 1.757658e-03 5.449710e-03
## 27: 3.341008e-03 1.123765e-02
## 28: 5.884704e-03 2.159805e-02
## 29: 9.631426e-03 3.882609e-02
## 30: 1.468603e-02 6.550916e-02
## 31: 2.091307e-02 1.040953e-01
## 32: 2.787488e-02 1.563116e-01
## 33: 3.485052e-02 2.225766e-01
## 34: 4.095139e-02 3.015955e-01
## 35: 4.531079e-02 3.903025e-01
## 36: 4.729020e-02 4.842169e-01
## 37: 4.663339e-02 5.781405e-01
## 38: 4.351702e-02 6.670133e-01
## 39: 3.848584e-02 7.466952e-01
## 40: 3.230222e-02 8.144844e-01
## 41: 2.576511e-02 8.692806e-01
## 42: 1.955464e-02 9.114195e-01
## 43: 1.413869e-02 9.422856e-01
## 44: 9.750086e-03 9.638457e-01
## 45: 6.419825e-03 9.782225e-01
## 46: 4.040242e-03 9.873841e-01
## 47: 2.432735e-03 9.929691e-01
## 48: 1.402816e-03 9.962292e-01
## 49: 7.753943e-04 9.980530e-01
## 50: 4.111890e-04 9.990318e-01
## 51: 2.093741e-04 9.995360e-01
## 52: 1.024515e-04 9.997857e-01
## 53: 4.821301e-05 9.999045e-01
## 54: 2.183664e-05 9.999589e-01
## 55: 9.525658e-06 9.999830e-01
## 56: 4.004908e-06 9.999932e-01
## 57: 1.623933e-06 9.999974e-01
## 58: 6.354796e-07 9.999990e-01
## 59: 2.401389e-07 9.999996e-01
## 60: 8.768217e-08 9.999999e-01
## 61: 3.095279e-08 1.000000e+00
##              PMF          CDF

#8
dpois(65, lambda = 71*24)

## [1] 0

Normal distribution

Percentile computation

##Calculating 3sd for each side

vector <- c()
mean <- 6493
deviation <- 1890
for (i in 0:300){
  result <- 6493-(0.01*i)*1890
  vector[300-i] <- result
}
for (i in 1:300){
  result <- 6493+(0.01*i)*1890
  vector[i+300] <- result
}
print(vector)

##   [1]   841.9   860.8   879.7   898.6   917.5   936.4   955.3   974.2
##   [9]   993.1  1012.0  1030.9  1049.8  1068.7  1087.6  1106.5  1125.4
##  [17]  1144.3  1163.2  1182.1  1201.0  1219.9  1238.8  1257.7  1276.6
##  [25]  1295.5  1314.4  1333.3  1352.2  1371.1  1390.0  1408.9  1427.8
##  [33]  1446.7  1465.6  1484.5  1503.4  1522.3  1541.2  1560.1  1579.0
##  [41]  1597.9  1616.8  1635.7  1654.6  1673.5  1692.4  1711.3  1730.2
##  [49]  1749.1  1768.0  1786.9  1805.8  1824.7  1843.6  1862.5  1881.4
##  [57]  1900.3  1919.2  1938.1  1957.0  1975.9  1994.8  2013.7  2032.6
##  [65]  2051.5  2070.4  2089.3  2108.2  2127.1  2146.0  2164.9  2183.8
##  [73]  2202.7  2221.6  2240.5  2259.4  2278.3  2297.2  2316.1  2335.0
##  [81]  2353.9  2372.8  2391.7  2410.6  2429.5  2448.4  2467.3  2486.2
##  [89]  2505.1  2524.0  2542.9  2561.8  2580.7  2599.6  2618.5  2637.4
##  [97]  2656.3  2675.2  2694.1  2713.0  2731.9  2750.8  2769.7  2788.6
## [105]  2807.5  2826.4  2845.3  2864.2  2883.1  2902.0  2920.9  2939.8
## [113]  2958.7  2977.6  2996.5  3015.4  3034.3  3053.2  3072.1  3091.0
## [121]  3109.9  3128.8  3147.7  3166.6  3185.5  3204.4  3223.3  3242.2
## [129]  3261.1  3280.0  3298.9  3317.8  3336.7  3355.6  3374.5  3393.4
## [137]  3412.3  3431.2  3450.1  3469.0  3487.9  3506.8  3525.7  3544.6
## [145]  3563.5  3582.4  3601.3  3620.2  3639.1  3658.0  3676.9  3695.8
## [153]  3714.7  3733.6  3752.5  3771.4  3790.3  3809.2  3828.1  3847.0
## [161]  3865.9  3884.8  3903.7  3922.6  3941.5  3960.4  3979.3  3998.2
## [169]  4017.1  4036.0  4054.9  4073.8  4092.7  4111.6  4130.5  4149.4
## [177]  4168.3  4187.2  4206.1  4225.0  4243.9  4262.8  4281.7  4300.6
## [185]  4319.5  4338.4  4357.3  4376.2  4395.1  4414.0  4432.9  4451.8
## [193]  4470.7  4489.6  4508.5  4527.4  4546.3  4565.2  4584.1  4603.0
## [201]  4621.9  4640.8  4659.7  4678.6  4697.5  4716.4  4735.3  4754.2
## [209]  4773.1  4792.0  4810.9  4829.8  4848.7  4867.6  4886.5  4905.4
## [217]  4924.3  4943.2  4962.1  4981.0  4999.9  5018.8  5037.7  5056.6
## [225]  5075.5  5094.4  5113.3  5132.2  5151.1  5170.0  5188.9  5207.8
## [233]  5226.7  5245.6  5264.5  5283.4  5302.3  5321.2  5340.1  5359.0
## [241]  5377.9  5396.8  5415.7  5434.6  5453.5  5472.4  5491.3  5510.2
## [249]  5529.1  5548.0  5566.9  5585.8  5604.7  5623.6  5642.5  5661.4
## [257]  5680.3  5699.2  5718.1  5737.0  5755.9  5774.8  5793.7  5812.6
## [265]  5831.5  5850.4  5869.3  5888.2  5907.1  5926.0  5944.9  5963.8
## [273]  5982.7  6001.6  6020.5  6039.4  6058.3  6077.2  6096.1  6115.0
## [281]  6133.9  6152.8  6171.7  6190.6  6209.5  6228.4  6247.3  6266.2
## [289]  6285.1  6304.0  6322.9  6341.8  6360.7  6379.6  6398.5  6417.4
## [297]  6436.3  6455.2  6474.1  6493.0  6511.9  6530.8  6549.7  6568.6
## [305]  6587.5  6606.4  6625.3  6644.2  6663.1  6682.0  6700.9  6719.8
## [313]  6738.7  6757.6  6776.5  6795.4  6814.3  6833.2  6852.1  6871.0
## [321]  6889.9  6908.8  6927.7  6946.6  6965.5  6984.4  7003.3  7022.2
## [329]  7041.1  7060.0  7078.9  7097.8  7116.7  7135.6  7154.5  7173.4
## [337]  7192.3  7211.2  7230.1  7249.0  7267.9  7286.8  7305.7  7324.6
## [345]  7343.5  7362.4  7381.3  7400.2  7419.1  7438.0  7456.9  7475.8
## [353]  7494.7  7513.6  7532.5  7551.4  7570.3  7589.2  7608.1  7627.0
## [361]  7645.9  7664.8  7683.7  7702.6  7721.5  7740.4  7759.3  7778.2
## [369]  7797.1  7816.0  7834.9  7853.8  7872.7  7891.6  7910.5  7929.4
## [377]  7948.3  7967.2  7986.1  8005.0  8023.9  8042.8  8061.7  8080.6
## [385]  8099.5  8118.4  8137.3  8156.2  8175.1  8194.0  8212.9  8231.8
## [393]  8250.7  8269.6  8288.5  8307.4  8326.3  8345.2  8364.1  8383.0
## [401]  8401.9  8420.8  8439.7  8458.6  8477.5  8496.4  8515.3  8534.2
## [409]  8553.1  8572.0  8590.9  8609.8  8628.7  8647.6  8666.5  8685.4
## [417]  8704.3  8723.2  8742.1  8761.0  8779.9  8798.8  8817.7  8836.6
## [425]  8855.5  8874.4  8893.3  8912.2  8931.1  8950.0  8968.9  8987.8
## [433]  9006.7  9025.6  9044.5  9063.4  9082.3  9101.2  9120.1  9139.0
## [441]  9157.9  9176.8  9195.7  9214.6  9233.5  9252.4  9271.3  9290.2
## [449]  9309.1  9328.0  9346.9  9365.8  9384.7  9403.6  9422.5  9441.4
## [457]  9460.3  9479.2  9498.1  9517.0  9535.9  9554.8  9573.7  9592.6
## [465]  9611.5  9630.4  9649.3  9668.2  9687.1  9706.0  9724.9  9743.8
## [473]  9762.7  9781.6  9800.5  9819.4  9838.3  9857.2  9876.1  9895.0
## [481]  9913.9  9932.8  9951.7  9970.6  9989.5 10008.4 10027.3 10046.2
## [489] 10065.1 10084.0 10102.9 10121.8 10140.7 10159.6 10178.5 10197.4
## [497] 10216.3 10235.2 10254.1 10273.0 10291.9 10310.8 10329.7 10348.6
## [505] 10367.5 10386.4 10405.3 10424.2 10443.1 10462.0 10480.9 10499.8
## [513] 10518.7 10537.6 10556.5 10575.4 10594.3 10613.2 10632.1 10651.0
## [521] 10669.9 10688.8 10707.7 10726.6 10745.5 10764.4 10783.3 10802.2
## [529] 10821.1 10840.0 10858.9 10877.8 10896.7 10915.6 10934.5 10953.4
## [537] 10972.3 10991.2 11010.1 11029.0 11047.9 11066.8 11085.7 11104.6
## [545] 11123.5 11142.4 11161.3 11180.2 11199.1 11218.0 11236.9 11255.8
## [553] 11274.7 11293.6 11312.5 11331.4 11350.3 11369.2 11388.1 11407.0
## [561] 11425.9 11444.8 11463.7 11482.6 11501.5 11520.4 11539.3 11558.2
## [569] 11577.1 11596.0 11614.9 11633.8 11652.7 11671.6 11690.5 11709.4
## [577] 11728.3 11747.2 11766.1 11785.0 11803.9 11822.8 11841.7 11860.6
## [585] 11879.5 11898.4 11917.3 11936.2 11955.1 11974.0 11992.9 12011.8
## [593] 12030.7 12049.6 12068.5 12087.4 12106.3 12125.2 12144.1 12163.0

Exercises done by hand:

1- 0,00427

2- ~1 (0,999.. rounded to 1)

3- 92,43

4- 91,41

5- 0,9984

6- X ~ Norm(92,3.4)

7- 89,97

8- 0,1928

Using R

1-pnorm(101,92,3.4) # Exercise 1

## [1] 0.004059761

1-pnorm(78.7,92,3.4) # Exercise 2 

## [1] 0.9999542

qnorm(0.128, 92, 3.4,lower.tail = TRUE) # Exercise 3 

## [1] 88.13795

qnorm(0.171,92,3.4) # Exercise 4 

## [1] 88.76925

pnorm(102,92,3.4)  # Exercise 5 

## [1] 0.9983652

qnorm(0.275,92,3.4) # Exercise 7

## [1] 89.96762

pnorm(89.1,92,3.4)-pnorm(83,92,3.4)  # Exercise 8 

## [1] 0.1927862