HW4
Ali Harb
October 14, 2018
Libraries
Required libraries to perform the assignment
library(knitr)
library(DT)
library(cowplot)
library(stringr)
library(ggplot2)Q3- P19.1
Using Monte Carlo simulation, write an algorithm to calculate an approximation to pi by considering the number of random points selected inside the quarter circle.
ppi <- function(iterations) {
x <- runif(iterations, 0, 1)
y <- runif(iterations, 0, 1)
Q <- x^2 + y^2 < 1
estimate <- (sum(Q)/iterations) * 4
return (estimate)
}
results<-data.frame(Points = numeric(), A_Estimate = numeric())
for(i in seq(1,5, by=1)){
points<-10^i
temp<-ppi(points)
results[i,1]<-points
results[i,2]<-temp
}
datatable(results,options = list(pageLength = 5),class = 'cell-border stripe')Q1 -P194
Use the middle-square method to generate for the following
10 random numbers using x0 D 1009.
20 random numbers using x0 D 653217.
15 random numbers using x0 D 3043.
Comment about the results of each sequence. Was there cycling? Did each sequence degenerate rapidly?
middleSquare <- function(x,len) {
sq <- x^2
len_sq <- len* 2
sq <- as.character(stringr::str_pad(as.character(sq), len_sq, pad = "0"))
start <- len_sq/2 - len/2 + 1
end <- len_sq/2 + len/2
return (as.numeric(substring(sq, start, end)))
}
GetMiddelSquare<-function(Var, nTimes){
Val=Var
len<-nchar(as.character(Val))
for(i in 1:nTimes){
Val[i+1]<-middleSquare(Val[i],len)
}
return(Val)
}a.10 random numbers using x0 D 1009
GetMiddelSquare(1009,10)## [1] 1009 180 324 1049 1004 80 64 40 16 2 0- 20 random numbers using x0 D 653217
GetMiddelSquare(653217,20)## [1] 653217 692449 485617 823870 761776 302674 611550 993402 847533 312186
## [11] 460098 690169 333248 54229 940784 74534 555317 376970 106380 316704
## [21] 301423- 15 random numbers using x0 D 3043.
GetMiddelSquare(3043,15)## [1] 3043 2598 7496 1900 6100 2100 4100 8100 6100 2100 4100 8100 6100 2100
## [15] 4100 8100Comment about the results of each sequence. Was there cycling? Did each sequence degenerate rapidly?
The sequence generated shows that small seed number will produce a very short data points before it get to zero. The number in the sequence does not cycle.
It generate a non-replicated sequence of numbers and the digit and the values are not approaching zero.
It often regenerates the same numbers and cycle some numbers in the sequence.
Q4-P201
Horse Race–Construct and perform a Monte Carlo simulation of a horse race. You can be creative and use odds from the newspaper, or simulate the Mathematical Derby with the entries and odds shown in following table.
## # Mathematical Derby
## | Entry's name | Odds |
## |-----------------|:-------------:|
## |Euler's Folly | 7-1 |
## |Leapin' Leibniz | 5-1 |
## |Newton Lobell | 9-1 |
## |Count Cauchy | 12-1 |
## |Pumped up Poisson| 4-1 |
## |Loping L'Hôpital | 35-1 |
## |Steamin' Stokes | 15-1 |
## |Dancin' Dantzig | 4-1 |Construct and perform a Monte Carlo simulation of 1000 horse races. Which horse won the most races? Which horse won the fewest races? Do these results surprise you? Provide the tallies of how many races each horse won with your output.
horse <- c("Euler's_Folly","Leapin_Leibniz","Newton_Lobell","Count_Cauchy","Pumped_up Poisson","Loping_L'Hopital","Steamin_Stokes","Dancin_Dantzig")
probs <- c(1/7,1/5,1/9,1/12,1/4,1/35,1/15,1/4)
nRaces=1000
mydf <- data.frame(number = 1:8,horses=horse, prob=round(probs,3))
simulation <- sample(x = mydf$horses, size = nRaces, prob = mydf$prob, replace = T)
r=simulation
SimulationVal=0
for(i in 1:length(summary(r))){
SimulationVal[i]<-as.numeric(summary(r)[[i]][[1]])
}
Dispaly <- data.frame(Won = summary(r))
Dispaly$horse <- rownames(Dispaly)
rownames(Dispaly) <- NULL
Dispaly## Won horse
## 1 71 Count_Cauchy
## 2 229 Dancin_Dantzig
## 3 131 Euler's_Folly
## 4 174 Leapin_Leibniz
## 5 22 Loping_L'Hopital
## 6 102 Newton_Lobell
## 7 223 Pumped_up Poisson
## 8 48 Steamin_StokestheSum<-sum(SimulationVal)
for(i in 1:length(SimulationVal)){
SimulationVal[i]<-(SimulationVal[i]/theSum)*100
}
mydf$SimulationVal<-SimulationVal
#print(SimulationVal)
ggplot(mydf, aes(x = mydf$horses, y = mydf$SimulationVal, fill=mydf$SimulationVal)) +
geom_bar(stat = "identity",colour="blue") +
labs(color = "Number of gears")+
xlab("Horses") +
ylab("Frequency") +
ggtitle("Distribution") +
geom_text(aes(label = sprintf("%.2f%%",mydf$SimulationVal)),
vjust = 1.5,color="green")+ theme(axis.text.x=element_text(angle=90,vjust=0.75))
Q3-P211
In many situations, the time T between deliveries and the order quantity Q is not xed. Instead, an order is placed for a speci c amount of gasoline. Depending on how many orders are placed in a given time interval, the time to ll an order varies. You have no reason to believe that the performance of the delivery operation will change. Therefore,
you have examined records for the past 100 deliveries and found the following lag times, or extra days, required to ll your order:
Lag time Number of occurrences
(in days)
5 20
6 13
7 2
Total: 100
## # Mathematical Derby
## |Lag time (in days) | Number of occurrences |
## |---------------------|:---------------------------:|
## |2 | 10 |
## |3 | 25 |
## |4 | 30 |
## |5 | 20 |
## |6 | 13 |
## |7 | 2 |
## | Total: | 100 |Construct a Monte Carlo simulation for the lag time submodel. If you have a handheld calculator or computer available, test your submodel by running 1000 trials and com- paring the number of occurrences of the various lag times with the historical data.
lag.times <- as.factor(2:7)
lag.prob <- c(0.1, 0.25, 0.3, 0.2, 0.13, 0.02)
trials = 1000
lag.simulation <- sample(x = lag.times, size = trials, prob = lag.prob, replace = T)
lag.results <- as.data.frame(table(lag.simulation))
lag.results <- cbind(lag.results, lag.prob)
new.names <- c("Lag_Time", "Simulated", "Actual")
colnames(lag.results) <- new.names
lag.results$Simulated <- lag.results$Simulated / trials
lag.results$Difference <- abs(lag.results$Actual - lag.results$Simulated)
datatable(lag.results)Q2 -P222
If the waiting time of a person is the time the person stands in a queuethe time from arrival at the lobby until entry into an available elevatorwhat are the average and maximum times a person waits in a queue?
arrivalL <- 0:12
arrivalL <- arrivalL * 10 + 15
arrivalU<- arrivalL + 9
arrivalU[length(arrivalU)] <- 145
arrivalMid <- (arrivalL + arrivalU) / 2
arrivalOccur <- c(11, 35, 42, 61, 108, 193, 240, 207, 150, 85, 44, 21, 3)
unloadL <- 0:8
unloadL <- unloadL * 5 + 45
unloadU <- unloadL + 4
unloadU[length(unloadU)] <- 90
unloadMid <- (unloadL+ unloadU) / 2
unloadOccur <- c(20, 54, 114, 103, 156, 221, 250, 171, 109)probs
arrivalProbs <- arrivalOccur / sum(arrivalOccur)
arrivalCum <- cumsum(arrivalProbs)
unloadProbs <- unloadOccur / sum(unloadOccur)
unloadCum <- cumsum(unloadProbs)model
arrivalLinear <- approxfun(arrivalCum, arrivalMid, method="linear")
unloadLinear <- approxfun(unloadCum, unloadMid, method="linear")ArrivalLoad <- cbind(arrivalCum,arrivalMid)
ArrivalUnload <- cbind(unloadCum,unloadMid)
str(as.data.frame(ArrivalLoad))## 'data.frame': 13 obs. of 2 variables:
## $ arrivalCum: num 0.00917 0.03833 0.07333 0.12417 0.21417 ...
## $ arrivalMid: num 19.5 29.5 39.5 49.5 59.5 ...library(cowplot)
p1 <- ggplot(as.data.frame(ArrivalLoad), aes(x = arrivalCum, y = arrivalMid)) +
geom_line() +
geom_point(shape = 1, fill = "white", size = 2, stroke = 1) +
theme(legend.position="none") +
geom_point(aes(colour=factor(arrivalLinear(arrivalCum))))
p2 <- ggplot(as.data.frame(ArrivalUnload), aes(x = unloadCum, y = unloadMid)) +
geom_line() +
geom_point(shape = 1, fill = "white", size = 2, stroke = 1) +
theme(legend.position="none") +
geom_point(aes(colour=factor(arrivalLinear(unloadCum))))
#plot(unloadCum, unloadMid, main="Unload Data")
#points(unloadCum, unloadLinear(unloadCum), type="l")
plot_grid(p1, p2, labels = "AUTO")
arrivalLm <- lm(arrivalMid ~ poly(arrivalCum, 3))
unloadLm <- lm(unloadMid ~ poly(unloadCum, 3))
par(mfrow=c(1,2))
plot(arrivalCum, arrivalMid, main="Arrival Data")
points(arrivalCum, arrivalLinear(arrivalCum), type="l")
SimData <- data.frame(arrivalCum=seq(from=0, to=1, by=0.01))
y <- predict(arrivalLm, newdata=SimData)
points(SimData$arrivalCum, y, type="l", col="red")
plot(unloadCum, unloadMid, main="Unload Data")
points(unloadCum, unloadLinear(unloadCum), type="l")
newData <- data.frame(unloadCum=seq(from=0, to=1, by=0.01))
y <- predict(unloadLm, newdata=newData)
points(newData$unloadCum, y, type="l", col="red")