data604_assign3
Ali Harb
September 21, 2017
Libraries
library(psych)3.5 Problem
3.2.1 - 50 throws of two dice
TwoDiceRGN<-function(n){
dice1<-sample(1:6, n, replace = T)
dice2<-sample(1:6, n, replace = T)
sumDice<-dice1+dice2
dicedf<-cbind(dice1,dice2)
dicedf<-cbind(dicedf,sumDice)
return(dicedf)
}
n=50
mydf<-TwoDiceRGN(n)
SumOfTwoDice<-mydf[,3]
describe(SumOfTwoDice)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 50 6.28 2.7 6 6.22 2.97 2 12 10 0.22 -0.84 0.38barplot(table(SumOfTwoDice))
1) Extend the simulation of throwing two dice
a) 500 throws of two dice
n=500
mydf<-TwoDiceRGN(n)
SumOfTwoDice<-mydf[,3]
describe(SumOfTwoDice)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 500 7.13 2.46 7 7.12 2.97 2 12 10 0.01 -0.72 0.11barplot(table(SumOfTwoDice))
b) Load the dice by changing the probabiliies of the faces to be something other than uniform at 1/16 for each face.
TwoDiceRGNProb <- function(n,prob1,prob2){
dice1<-sample(1:6, n, replace = T,prob =prob1)
dice2<-sample(1:6, n, replace = T,prob =prob2)
sumDice<-dice1+dice2
dicedf<-cbind(dice1,dice2)
dicedf<-cbind(dicedf,sumDice)
return(dicedf)
}
p1 <- c(2/8,1/10,3/10,4/16,4/10,3/8)
p2 <- c(3/8,5/10,4/10,2/8,3/10,4/10)
n=500
mydf<-TwoDiceRGNProb(n,p1,p2)
SumOfTwoDice<-mydf[,3]
describe(SumOfTwoDice)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 500 7.23 2.5 7 7.25 2.97 2 12 10 -0.1 -0.62 0.11par(mar=c(1,1,1,1))
hist(SumOfTwoDice)
c) Use @RISK or another Excel add-in for static Monte Carlo spreadsheet simulation, to make 10,000 throws of the pair of fair dice, and compare your results to the true probabilities of getting the sum equal to each of 2, 3,…, 12 as well as the true expected value of 7.
n=10000
df<-0
mydf<-TwoDiceRGN(n)
SumOfTwoDice<-mydf[,3]
Prop<-c(1/36, 2/36, 3/36, 4/36, 5/36, 6/36, 5/36, 4/36, 3/36, 2/36, 1/36)
Prop<-round(Prop,4)
R_SimUlation<-table(SumOfTwoDice)/length(SumOfTwoDice)
PropDiffrences<-Prop-R_SimUlation
df<-cbind(Prop,R_SimUlation)
df<-cbind(df,PropDiffrences)
df## Prop R_SimUlation PropDiffrences
## 2 0.0278 0.0265 0.0013
## 3 0.0556 0.0590 -0.0034
## 4 0.0833 0.0822 0.0011
## 5 0.1111 0.1110 0.0001
## 6 0.1389 0.1367 0.0022
## 7 0.1667 0.1648 0.0019
## 8 0.1389 0.1393 -0.0004
## 9 0.1111 0.1095 0.0016
## 10 0.0833 0.0865 -0.0032
## 11 0.0556 0.0569 -0.0013
## 12 0.0278 0.0276 0.0002As the result shown, the simulation will behave accordingly with true probability and produce a normal distribution. The values defer slightly at the tails and these differences get smaller and smaller as we get close to the center of the distribution.
5) We will performed the Monte Carlo Integration for the following function h(x) to be the probability density function of a normal variable with mean miu and standard deviation sigma > 0.
MyFun <- function(x, a, b, u, sigma){
q <- 1/(sqrt(2*pi*(sigma^2)))
e <- -((x-u)^2)/(2*(sigma^2))
result <- (b-a)*q*exp(e)
return(result)
}
getConfidence<-function(interval,intervalValue,meanVal,sigma,n){
confidence<-round(qnorm(intervalValue)*sigma/sqrt(n),4)
lowerTail<-round(meanVal-confidence,4)
upperTail<-round(meanVal+confidence,4)
results<-data.frame(interval,confidence,lowerTail,upperTail)
return(results)
}
n <- 50
a <- 4.5
b <- 6.7
u <- 5.8
sigma <- 2.3
rand <- runif(n, a, b)
p <- sapply(rand, MyFun, a=a, b=b, u=u, sigma=sigma)
monte_est <- mean(p)
True_est <- pnorm(b, mean=u, sd=sigma)-pnorm(a,mean=u, sd=sigma)
esd<- sd(p)
df_est<-data.frame(monte_est,True_est,esd)
colnames(df_est)<-c("Monte Carlo Estimate","Exact Estimate","SD Estimate")
confidence90<-getConfidence(90,0.95,monte_est,esd,n)
df<-as.data.frame(confidence90)
confidence95<-getConfidence(95,0.975,monte_est,esd,n)
df<-rbind(df,confidence95)
confidence99<-getConfidence(99,0.995,monte_est,esd,n)
df<-rbind(df,confidence99)
df## interval confidence lowerTail upperTail
## 1 90 0.0036 0.3619 0.3691
## 2 95 0.0043 0.3612 0.3698
## 3 99 0.0057 0.3598 0.3712The result shows a very small variation within an acceptable range of fluctuation
17) For the goods and demand limits, create a spreadsheet that will simulate costs, revenue and profit for a 90 day season. Assume that demand is always met.
Given:
products:oats,peas,beans,barley
Buying prices:1.05,3.17,1.99, 0.95
sells prices: 1.29,3.76,2.23, 1.65
min demand: 0 ,0 ,0 ,0
max demand: 10 ,8 ,14 ,14
Find:
- selling simulation for 90 days
for one day and the whole season simulation: total cost, total revenue, total profit
- cost<-sum(product cost)
- revenue<-sum(product sells)
- profit<-sum(evenue-cost)
options(width = 300)
Product.info<- data.frame(
"Pruduct" = c("Oats","Peas","Beans","Barley"),
"Cost" = c(1.05,3.17,1.99,0.95),
"Sale" = c(1.29,3.76,2.23,1.65),
"Max_Demand"=c(10,8,14,11)
)
days<-90
ProductVarRGN<-function(dataframe,days){
dailySimu<-data.frame(matrix(NA,nrow=days, ncol=0))
n <-nrow(dataframe)
for( i in 1:n){
pdf<-ceiling(runif(days, min=0, max=dataframe[i ,"Max_Demand"]))
dailySimu<-cbind(dailySimu,pdf)
}
colnames(dailySimu)<-c("Oats","Peas","Beans","Barley")
return(dailySimu)
}
getdailyProfit<-function(RGNdf,df,days){
dailySimuCost<-data.frame(matrix(NA,nrow=days, ncol=0))
dailySimuSale<-data.frame(matrix(NA,nrow=days, ncol=0))
for(j in 1:length(RGNdf)){
dailyCost<-RGNdf[j]*df[j,"Cost"]
dailySimuCost<-cbind(dailySimuCost,dailyCost)
dailyRevenue<-RGNdf[j]*df[j,"Sale"]
dailySimuSale<-cbind(dailySimuSale,dailyRevenue)
}
costSum<-c()
saleSum<-c()
Profit<-c()
n<-nrow(RGNdf)
for(i in 1:n){
costSum[i]<-sum(dailySimuCost[i,])
saleSum[i]<-sum(dailySimuSale[i,])
}
Profit<-saleSum-costSum
mydf<-data.frame(dailySimuCost,costSum,dailySimuSale,saleSum,Profit)
colnames(mydf)<-c("Oats_Cx","Peas_Cx","Beans_Cx","Barley_Cx","Cost",
"Oats_Rx","Peas_Rx","Beans_Rx","Barley_Rx","Revenue",
"Profit")
return(mydf)
}
RGN<-ProductVarRGN(Product.info,days)
DailySeasonProfit<-getdailyProfit(RGNdf=RGN,df=Product.info,days)
totalSeasonCost<-sum(DailySeasonProfit[,5])
totalSeasonRevenue<-sum(DailySeasonProfit[,10])
totalSeasonProfit<-sum(DailySeasonProfit[,11])
SimulationResult<-data.frame("Variables" = c("Total_Revnue","Total_Cost","Total_Profit"),
"Amount" = c(totalSeasonRevenue,totalSeasonCost,totalSeasonProfit)
)
Product.info## Pruduct Cost Sale Max_Demand
## 1 Oats 1.05 1.29 10
## 2 Peas 3.17 3.76 8
## 3 Beans 1.99 2.23 14
## 4 Barley 0.95 1.65 11head(DailySeasonProfit,15)## Oats_Cx Peas_Cx Beans_Cx Barley_Cx Cost Oats_Rx Peas_Rx Beans_Rx Barley_Rx Revenue Profit
## 1 4.20 3.17 5.97 10.45 23.79 5.16 3.76 6.69 18.15 33.76 9.97
## 2 10.50 15.85 27.86 8.55 62.76 12.90 18.80 31.22 14.85 77.77 15.01
## 3 10.50 25.36 11.94 2.85 50.65 12.90 30.08 13.38 4.95 61.31 10.66
## 4 3.15 9.51 17.91 7.60 38.17 3.87 11.28 20.07 13.20 48.42 10.25
## 5 3.15 15.85 1.99 6.65 27.64 3.87 18.80 2.23 11.55 36.45 8.81
## 6 2.10 22.19 15.92 1.90 42.11 2.58 26.32 17.84 3.30 50.04 7.93
## 7 4.20 12.68 3.98 7.60 28.46 5.16 15.04 4.46 13.20 37.86 9.40
## 8 10.50 9.51 11.94 1.90 33.85 12.90 11.28 13.38 3.30 40.86 7.01
## 9 5.25 12.68 11.94 8.55 38.42 6.45 15.04 13.38 14.85 49.72 11.30
## 10 9.45 22.19 9.95 10.45 52.04 11.61 26.32 11.15 18.15 67.23 15.19
## 11 3.15 19.02 17.91 10.45 50.53 3.87 22.56 20.07 18.15 64.65 14.12
## 12 10.50 6.34 23.88 6.65 47.37 12.90 7.52 26.76 11.55 58.73 11.36
## 13 9.45 15.85 9.95 7.60 42.85 11.61 18.80 11.15 13.20 54.76 11.91
## 14 3.15 3.17 19.90 5.70 31.92 3.87 3.76 22.30 9.90 39.83 7.91
## 15 6.30 6.34 27.86 6.65 47.15 7.74 7.52 31.22 11.55 58.03 10.88SimulationResult## Variables Amount
## 1 Total_Revnue 4488.06
## 2 Total_Cost 3576.51
## 3 Total_Profit 911.55barplot(SimulationResult$Amount, names = SimulationResult$Variables,
xlab = "Variables", ylab = "Amount",
main = "Walther's Sells Output")