Data 604 - HW #3
Brian Kreis
September 18, 2017
Chapter 3 of Kelton
Problem #1
I thought this was a pretty efficient and straight forward way to adjust the number of dice https://stackoverflow.com/questions/14820203/simulating-rolling-two-dice
#a
x = c(1, 2, 3, 4, 5, 6)
p = c(1/6, 1/6, 1/6, 1/6, 1/6, 1/6)
#replicated 5 times
replicate(5, mean(rowSums(replicate(2, sample(x, 500, replace=T, prob = p)))))## [1] 6.840 6.786 7.072 7.064 7.024#b
#change the probabilities of the die to load them
p = c(1/8, 1/8, 1/8, 1/8, 1/8, 3/8)
#replicated 5 times
replicate(5, mean(rowSums(replicate(2, sample(x, 500, replace=T, prob = p)))))## [1] 8.004 8.180 8.390 8.156 8.282#c
p = c(1/6, 1/6, 1/6, 1/6, 1/6, 1/6)
data <- data.frame(rowSums(replicate(2, sample(x, 10000, replace=T, prob = p))))
colnames(data) <- c("rolls")
library(dplyr)
library(tidyr)
df <- data %>% count(rolls)
df <- df %>% mutate(calc_prob = n/(sum(df$n))) %>%
mutate(expected = (c(1, 2, 3, 4, 5, 6, 5, 4, 3, 2, 1)/36)) %>%
gather (type, probability, calc_prob:expected)
library (ggplot2)
ggplot(df, aes(x = rolls, y = probability, fill = type)) +
geom_histogram(stat = "identity", position = "dodge") +
ggtitle("Expected versus Calculated") +
labs(y= "Probability", x = "Rolled Value") +
scale_x_continuous(breaks=c(2:12))
Problem #5
a <- 4.5
b <- 6.7
n <- 50
#50 random values distributed continuously between a and b
samp <- data.frame(runif(50, a, b))
colnames(samp) <- c("samples")
#caluclate the mean and sd
stDev <- sd(samp$samples)
meanV <- mean(samp$samples)
#
#calculate the density
samp <- samp %>%
mutate(density = (b-a) * dnorm(samples, mean=5.8, sd =2.3))
samp## samples density
## 1 5.342030 0.3741067
## 2 4.880747 0.3523042
## 3 5.459709 0.3774432
## 4 4.585160 0.3319127
## 5 5.383179 0.3753817
## 6 6.070317 0.3789705
## 7 4.768871 0.3451127
## 8 5.668491 0.3809737
## 9 4.894681 0.3531517
## 10 6.457469 0.3663203
## 11 6.276078 0.3735091
## 12 4.928762 0.3551786
## 13 5.078721 0.3632868
## 14 6.087756 0.3786221
## 15 6.073886 0.3789010
## 16 6.579793 0.3602833
## 17 5.922535 0.3810558
## 18 4.897022 0.3532931
## 19 4.955876 0.3567434
## 20 6.080023 0.3787792
## 21 4.905843 0.3538228
## 22 5.660787 0.3808986
## 23 6.531762 0.3627642
## 24 5.917276 0.3811012
## 25 5.029968 0.3607989
## 26 5.271017 0.3716367
## 27 5.426523 0.3765991
## 28 5.167421 0.3674338
## 29 5.193769 0.3685691
## 30 4.620478 0.3345762
## 31 5.065274 0.3626151
## 32 4.664031 0.3377806
## 33 5.892730 0.3812870
## 34 5.429015 0.3766651
## 35 5.196533 0.3686856
## 36 5.413029 0.3762340
## 37 6.510118 0.3638358
## 38 5.508725 0.3785492
## 39 6.598185 0.3592964
## 40 4.791565 0.3466258
## 41 5.265341 0.3714247
## 42 4.699701 0.3403369
## 43 5.366121 0.3748672
## 44 6.137034 0.3775219
## 45 5.892436 0.3812889
## 46 5.963262 0.3806368
## 47 6.301148 0.3726453
## 48 4.680301 0.3389543
## 49 5.624122 0.3804829
## 50 6.605564 0.3588947#Monte Carlo Estimate
mean(samp$density)## [1] 0.3660432#Exact
pnorm(b, mean=5.8, sd=2.3) - pnorm(a, mean=5.8, sd=2.3)## [1] 0.3662509me <- qnorm(0.95) * (stDev/(sqrt(n)))
lower <- mean(samp$density) - me
upper <- mean(samp$density) + me
paste(lower, upper)## [1] "0.223762614200363 0.508323717102806"me <- qnorm(0.90) * (stDev/(sqrt(n)))
lower <- mean(samp$density) - me
upper <- mean(samp$density) + me
paste(lower, upper)## [1] "0.25518839993191 0.476897931371259"# Our confidence intervals cover the exact integral after multiple triesProblem 17
oats = floor(runif(90, 0, 10))
peas = floor(runif(90, 0, 8))
beans = floor(runif(90, 0, 14))
barley = floor(runif(90, 0, 11))
df <- data.frame(oats,peas, beans, barley)
df$cost <- 1.05*oats + 3.17*peas + 1.99 *beans + 0.95*barley
df$revenue <- 1.29*oats + 3.76*peas + 2.23 *beans + 1.65*barley
df$profit <- df$revenue - df$cost
df## oats peas beans barley cost revenue profit
## 1 1 7 10 7 49.79 61.46 11.67
## 2 7 4 11 3 44.77 53.55 8.78
## 3 7 2 1 1 16.63 20.43 3.80
## 4 5 1 5 4 22.17 27.96 5.79
## 5 1 0 5 1 11.95 14.09 2.14
## 6 2 5 5 5 32.65 40.78 8.13
## 7 3 5 11 4 44.69 53.80 9.11
## 8 1 7 13 7 55.76 68.15 12.39
## 9 7 2 8 3 32.46 39.34 6.88
## 10 7 7 8 1 46.41 54.84 8.43
## 11 5 1 12 8 39.90 50.17 10.27
## 12 5 5 6 9 41.59 53.48 11.89
## 13 1 1 1 5 10.96 15.53 4.57
## 14 9 5 3 10 40.77 53.60 12.83
## 15 3 4 7 1 30.71 36.17 5.46
## 16 6 7 5 8 46.04 58.41 12.37
## 17 2 5 7 0 31.88 36.99 5.11
## 18 0 3 4 6 23.17 30.10 6.93
## 19 3 7 2 3 32.17 39.60 7.43
## 20 9 3 9 2 38.77 46.26 7.49
## 21 5 7 8 1 44.31 52.26 7.95
## 22 6 6 7 6 44.95 55.81 10.86
## 23 7 3 8 10 42.28 54.65 12.37
## 24 5 0 10 8 32.75 41.95 9.20
## 25 0 2 12 9 38.77 49.13 10.36
## 26 2 0 13 6 33.67 41.47 7.80
## 27 1 0 5 10 20.50 28.94 8.44
## 28 0 3 5 5 24.21 30.68 6.47
## 29 7 5 13 6 54.77 66.72 11.95
## 30 8 5 5 6 39.90 50.17 10.27
## 31 7 6 4 6 40.03 50.41 10.38
## 32 2 7 9 9 50.75 63.82 13.07
## 33 9 2 10 10 45.19 57.93 12.74
## 34 1 2 2 2 13.27 16.57 3.30
## 35 1 3 7 3 27.34 33.13 5.79
## 36 3 0 13 3 31.87 37.81 5.94
## 37 8 4 0 1 22.03 27.01 4.98
## 38 7 1 5 10 29.97 40.44 10.47
## 39 4 0 5 8 21.75 29.51 7.76
## 40 0 3 13 0 35.38 40.27 4.89
## 41 1 5 5 0 26.85 31.24 4.39
## 42 1 1 10 8 31.72 40.55 8.83
## 43 1 5 10 4 40.60 48.99 8.39
## 44 6 5 3 7 34.77 44.78 10.01
## 45 6 5 4 1 31.06 37.11 6.05
## 46 0 6 13 8 52.49 64.75 12.26
## 47 4 0 10 4 27.90 34.06 6.16
## 48 3 6 7 1 37.05 43.69 6.64
## 49 7 7 2 3 36.37 44.76 8.39
## 50 8 4 1 1 24.02 29.24 5.22
## 51 9 6 2 6 38.15 48.53 10.38
## 52 4 4 2 4 24.66 31.26 6.60
## 53 4 4 7 1 31.76 37.46 5.70
## 54 7 2 10 1 34.54 40.50 5.96
## 55 2 0 0 10 11.60 19.08 7.48
## 56 3 7 5 5 40.04 49.59 9.55
## 57 2 0 13 9 36.52 46.42 9.90
## 58 7 4 2 5 28.76 36.78 8.02
## 59 4 0 6 7 22.79 30.09 7.30
## 60 9 3 2 2 24.84 30.65 5.81
## 61 4 2 6 5 27.23 34.31 7.08
## 62 9 3 9 6 42.57 52.86 10.29
## 63 7 3 9 0 34.77 40.38 5.61
## 64 7 0 9 5 30.01 37.35 7.34
## 65 7 3 4 4 28.62 35.83 7.21
## 66 0 4 10 5 37.33 45.59 8.26
## 67 3 3 11 10 44.05 56.18 12.13
## 68 8 2 10 10 44.14 56.64 12.50
## 69 3 7 2 4 33.12 41.25 8.13
## 70 1 6 9 2 39.88 47.22 7.34
## 71 4 3 6 1 26.60 31.47 4.87
## 72 7 4 2 1 24.96 30.18 5.22
## 73 1 2 11 5 34.03 41.59 7.56
## 74 9 1 4 3 23.43 29.24 5.81
## 75 4 3 13 1 40.53 47.08 6.55
## 76 3 4 12 8 47.31 58.87 11.56
## 77 3 0 2 8 14.73 21.53 6.80
## 78 4 4 12 10 50.26 63.46 13.20
## 79 1 7 11 8 52.73 65.34 12.61
## 80 7 5 4 1 32.11 38.40 6.29
## 81 4 5 6 1 32.94 38.99 6.05
## 82 4 4 12 8 48.36 60.16 11.80
## 83 9 0 0 1 10.40 13.26 2.86
## 84 9 0 7 1 24.33 28.87 4.54
## 85 3 3 12 1 37.49 43.56 6.07
## 86 3 6 2 7 32.80 42.44 9.64
## 87 9 5 12 8 56.78 70.37 13.59
## 88 0 2 4 2 16.20 19.74 3.54
## 89 0 6 1 3 23.86 29.74 5.88
## 90 3 4 7 6 35.46 44.42 8.96#revenue, cost, profit
paste(sum(df$revenue), sum(df$cost), sum(df$profit))## [1] "3779.24 3050.45 728.79"