I have the following code for a Monte Carlo simulation:
function (muC1,sigmaC1,muC2,sigmaC2,Msim=10000) {
# a grand loop simulation program for the lifetime of
# a simple system of components. This system has
# two branches. The first branch fails when C1 fails.
# The second branch fails as soon as either of C2 or C3 fails.
#
# get alpha's and lambda's for Gamma lifetime RVs from the
# mu's and sigma's.
# Note C2 has the same distribution as C3, i.e.
# alpha3=alpha2, lambda3=lambda2
lambda1 = muC1/sigmaC1^2
alpha1 = muC1*lambda1
lambda2 = muC2/sigmaC2^2
alpha2 = muC2*lambda2
#
# initialize simulation summary holders
count.C1fail = 0 # number of times C1 fails before branch2
system.lifetime = rep(NA,Msim)
#
# begin grand loop
for (i in 1:Msim) {
C1 = rgamma(1,alpha1,lambda1)
C2 = rgamma(1,alpha2,lambda2)
C3 = rgamma(1,alpha2,lambda2)
branch2 = min(C2,C3)
system.lifetime[i] = max(C1,branch2)
if (C1 < branch2) count.C1fail = count.C1fail+1
} # end grand loop
#
# final summary calculations and wrapup
PC1fail = count.C1fail/Msim
hist(system.lifetime,main="Simulated System Lifetimes")
meanL = mean(system.lifetime)
stddevL = sd(system.lifetime)
MOEL95pct = 1.96*stddevL/sqrt(Msim)
out = list(muC1,sigmaC1,muC2,sigmaC2,Msim,PC1fail,meanL,MOEL95pct)
names(out) = c("muC1","sigmaC1","muC2","sigmaC2","Msim","PC1fail","meanL","MOEL95pct")
out
}
This simulation works perfectly for what I need, but for an assignment I need to run the same simulation (muC1 = 100, sigmaC1 = 20, muC2 = 80, sigmaC2 = 40, Msim = 10000) WITHOUT a grand loop. Any tips?
You can eliminate the loop by just drawing all your random values at once.
C1 <- rgamma(Msim, alpha1, lambda1)
C2 <- rgamma(Msim, alpha2, lambda2)
C3 <- rgamma(Msim, alpha2, lambda2)
branch2 <- pmin(C2, C3)
system.lifetime <- pmax(C1, branch2)
count.C1fail <- sum(C1 < branch2)
Then since we want to operate on vectors rather that just two values, we switch from min and max to pmin and pmax (the vectorized versions of the function). Finally, to count events in vectors we often just sum() over TRUE/FALSE values because TRUE is converted to 1 and FALSE 0 so we just get the count of TRUE events.
Doing it this way also means that it's not necessary to initialize count.C1fail or system.lifetime ahead of time.