Запуск симуляции без использования зацикливания в R - PullRequest
Новая
       363

Запуск симуляции без использования зацикливания в R

0 голосов
/ 18 февраля 2020

У меня есть следующий код для симуляции Монте-Карло: 

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 
}

Эта симуляция отлично работает для того, что мне нужно, но для задания мне нужно запустить ту же симуляцию (muC1 = 100, sigmaC1 = 20, muC2 = 80, sigmaC2 = 40, Msim = 10000) БЕЗ большого l oop. Какие-нибудь советы?

1 Ответ

1 голос
/ 19 февраля 2020

Вы можете устранить l oop, просто нарисовав все свои случайные значения сразу. 

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)

Тогда, поскольку мы хотим оперировать векторами, а не двумя значениями, мы переключаемся с min и max до pmin и pmax (векторизованные версии функции). Наконец, для подсчета событий в векторах мы часто просто sum() превышаем значения TRUE / FALSE, потому что TRUE преобразуется в 1 и FALSE 0, поэтому мы просто получаем количество событий TRUE.

Выполнение таким образом также означает, что нет необходимости инициализировать count.C1fail или system.lifetime заранее.

...