## CSUEB STATS 6310 - Stochastic Processes and Simulation - Spring 2017 - Prof B.Trumbo
## Homework Gui Larangeira 

## PSGS 4.6


set.seed(1212)
m = 100000;  lam1 = 1/5;  lam2 = 1/4     # constants
x1 = rexp(m, lam1);  x2 = rexp(m, lam2)  # simulation
v = pmin(x1, x2)                         # m-Vector of minimums
mean(x2==v)                   # x2==v is logical vector where TRUE if v, 
## [1] 0.55762
                              # the minimum, corresponds to x2 (assumed female). 
                              # Mean is the probability as usual.
                                       
5/9                           # Compare to the theoretical value                                    
## [1] 0.5555556
#b)
## Example 4.1
#  Scenario 1

set.seed(1212)
m = 100000;  lam1 = 1/5;  lam2 = 1/4     # constants
x1 = rexp(m, lam1);  x2 = rexp(m, lam2)  # simulation
v = pmin(x1, x2)                         # m-Vector of minimums
mean(v)                                  # approximates P{V > 5}
## [1] 2.22507
# Compare to theoretical 
20/9
## [1] 2.222222
#c)

## --
#  Scenario 2

set.seed(1212)
m = 100000;  lam1 = 1/5;  lam2 = 1/5
x1 = rexp(m, lam1);  x2 = rexp(m, lam2)
t = x1 + x2
mean(t > 5)
## [1] 0.73618
# We have to wait to begin being served (x1~Exp(1/5)) and to finish being served (x2~Exp(1/5))
# So it is equivalent to being served by the slower teller twice.

## PSGS 4.8

#a)
set.seed(12)
m = 100000                 # iterations
                           # SYSTEM A
n = 5;  lam = 1/4          # constants: system parameters 
x = rexp(m*n, lam)         # vector of all simulated data
DTA = matrix(x, nrow=m)    # each row a simulated system
w = apply(DTA, 1, max)     # lifetimes of m systems
mean(w)                    # approximates E(W)
## [1] 9.128704
mean(w > 12)                # approximates P{W > 12}
## [1] 0.22413
quantile(w, .9)            # approximates 90th percentile
##      90% 
## 15.51431
set.seed(12)               # SYSTEM B
nb = 4;  lamb = 1/5        # constants: system parameters (SYSTEM B) 
xb = rexp(m*nb, lamb)        # vector of all simulated data
DTAb = matrix(xb, nrow=m)    # each row a simulated system
wb = apply(DTAb, 1, max)     # lifetimes of m systems
mean(wb)                    # approximates E(W)
## [1] 10.4161
mean(wb > 12)                # approximates P{W > 12}
## [1] 0.31568
quantile(wb, .9)            # approximates 90th percentile
##      90% 
## 18.24178
# System B is more reliable than System A by all metrics. Higher mean leaftime, P{W>12} and 90% quantile.

ecdf = (1:m)/m
w.sort = sort(w); wb.sort = sort(wb)      # order statistics for system simulated above
plot(w.sort, ecdf, type="l", xlim=c(0,30), col="blue",xlab="Years")
lines(wb.sort, ecdf, col="red",type="l")
legend(22,0.2,c("System A","System B"),lty=c(1,1), lwd=c(2.5,2.5),col=c("blue","red"))