Hello,
I found an error in you z function. The code below is working for me.
M <- rnorm(n = 10, mean = 10, sd = 1)
S <- rnorm(n = 10, mean = 1, sd = 1)
Nrp <- 10
Nvr <- 10
Nsm <- 10
MCcrude<-function(M,S,P=NULL,DS=NULL,Nsm,Nrp){
Nvr=length(M)
pf=NULL
#pf <- numeric(Nrp)
for(k1 in seq(1:Nrp)){
X=NULL
#X <- numeric(Nrp)
for(j1 in 1:Nvr){
X=cbind(X,rnorm(Nsm,mean= M[j1],sd = S[j1]))
#message(X)
#X[j1] <- rnorm(Nsm,mean= M[j1],sd =S[j1])
}
z <- 30*X[,2]*X[,1]-X[,2]*X[,1]^2-30*X[,3]*X[,1]-X[,3]*X[,1] # error
pf=rbind(pf,length(which(z<=0))/Nsm) #probability of failure
#pf[Nrp] <- sum(z<=0)/Nsm
}
return(pf)
}
MCcrude(M = M, S = S, Nsm = 20, Nrp = 20)
Results:
> MCcrude(M = M, S = S, Nsm = 20, Nrp = 20)
[,1]
[1,] 0.85
[2,] 0.95
[3,] 0.95
[4,] 0.95
[5,] 1.00
[6,] 1.00
[7,] 0.90
[8,] 0.95
[9,] 0.95
[10,] 1.00
[11,] 0.95
[12,] 0.95
[13,] 0.95
[14,] 0.90
[15,] 0.95
[16,] 0.95
[17,] 0.90
[18,] 0.90
[19,] 1.00
[20,] 1.00
There were 50 or more warnings (use warnings() to see the first 50)