#---------------------------------------------------------------------------------#
#-------------------Algorithm FORM : old version of V. Moutoussamy----------------#
#---------------------------------------------------------------------------------#
FORMv0 <- function(f, u.dep, inputDist, N.calls, eps = 1e-7,
Method = "HLRF", IS = FALSE, q = 0.5, copula = "unif"){
# Few error messages :
if (mode(f) != "function") stop("f is not a function")
if( mode(inputDist) != "list") stop("inputDist is not a list")
if (missing(inputDist)) stop("inputDist is missing")
if (missing(N.calls)) stop("N.calls is missing")
ndim <- length(inputDist)
#-------------------------------------------------------------------------------#
#--------------- Create the list whith paramters of the input ------------------#
transformtionToInputSpace <- function(inputDist){
InputDist <- list()
InputDist <- inputDist
for(i in 1:ndim){
nparam <- length(inputDist[[i]][[2]])
for(j in 1:nparam){
InputDist[[i]]$q <- paste("q",InputDist[[i]][[1]], sep = "");
InputDist[[i]]$p <- paste("p",InputDist[[i]][[1]], sep = "");
InputDist[[i]]$d <- paste("d",InputDist[[i]][[1]], sep = "");
InputDist[[i]]$r <- paste("r",InputDist[[i]][[1]], sep = "");
}
}
return(InputDist)
}
InputDist <- transformtionToInputSpace(inputDist)
#------ Transformation to the Gaussian space ------#
Gtransf <- function(X){
ndim <- length(InputDist)
XU <- pnorm(X)
for(i in 1:ndim){
XU[i] <- do.call(InputDist[[i]]$q, c(list(XU[i, drop = FALSE]), InputDist[[i]][[2]]))
}
return(f(XU))
}
#------ Gradient------#
GRADfct <- function(func,X){
GG <- 0
df <- 0
RES <- list()
dfunc <- 0
RES[[1]] <- 0
f.X <- func(X)
RES.2 <- X
RES.3 <- f.X
for (i in 1:ndim) {
dX <- X
dX[i] <- dX[i] + 1e-6
f.dX <- func(dX)
dfunc[i] <- (f.dX-f.X)/1e-6
RES.2 <- rbind(RES.2, dX)
RES.3 <- rbind(RES.3, f.dX)
}
RES[[2]] <- RES.2
RES[[3]] <- RES.3
RES[[1]] <- dfunc
return(RES)
}
#----------------------------------------------------------------------#
x <- NULL # contains the input design
y <- NULL # contains f(x)
dy <- NULL # contains all the derivatives of f(x)
List.res <- list() # List of the output
u.new <- 0
cp <- 0 # Numbers of calls to f
fact.imp <- 0 # Importance factor
N.FORM <- floor(N.calls*q) # Number of calls reserved to FORM
u.ans <- Gtransf(u.dep)
cp <- cp + 1
x <- u.dep
y <- u.ans
res <- GRADfct(Gtransf, u.dep)
cp <- cp + ndim + 1
x <- rbind(x,res[[2]])
y <- rbind(y,res[[3]])
dy <- rbind(dy,t(res[[1]]))
g <- res[[1]]
err <- 1 #distance of two tested points
#----------------------------------Algorithme HLRF :----------------- #
if(Method == "HLRF"){
a <- g/(sqrt(g%*%g))
while( (err > eps) & (cp <= (N.FORM - (2 + ndim)) ) ){ # modif BIS
u.new <- (t(u.dep)%*%a)*a - u.ans*a/sqrt(g%*%g)
err <- sqrt(t(u.dep-u.new)%*%(u.dep-u.new))
u.dep <- u.new
res <- GRADfct(Gtransf, u.dep)
cp <- cp + ndim + 1
x <- rbind(x,res[[2]])
y <- rbind(y,res[[3]])
dy <- rbind(dy,t(res[[1]]))
u.ans <- res[[3]][1]
g <- res[[1]]
a <- ifelse( g == 0,c(0,0), g/(sqrt(g%*%g)) )
}
B <- -a%*%u.dep #indices of reliability
P <- pnorm(-B) #failure probability
fact.imp <- a^2
List.res <- list(P, B, cp, u.dep, fact.imp, x, y, dy);
names(List.res) <- c("pf", "beta", "compt.f", "design.point", "fact.imp", "x", "y", "dy")
if(IS == FALSE){return(List.res)}
}
#-----------------------------End of the algorithm HLRF :--------------------------#
#-----------------------------Algorithm Abdo-Rackwitz :--------------------------#
if(Method == "AR"){
lambda <- 2*(u.ans - (t(g)%*%u.dep))/(g%*%g)
while( (err > eps) & (cp <= (N.FORM - (2 + ndim)) ) ){ # modif BIS
u.new <- -lambda*g/2
err <- sqrt(t(u.dep-u.new)%*%(u.dep-u.new))
u.dep <- u.new
res <- GRADfct(Gtransf, u.dep)
cp <- cp + ndim + 1
x <- rbind(x, res[[2]])
y <- rbind(y, res[[3]])
dy <- rbind(dy,t(res[[1]]))
u.ans <- res[[3]][1]
g <- res[[1]]
lambda <- 2*(u.ans - (t(g)%*%u.dep))/(g%*%g)
if(cp >= N.FORM){break}
}
a <- g/(sqrt(g%*%g))
B <- -a%*%u.dep
P <- pnorm(-B)
fact.imp <- a^2
List.res <- list(P, B, cp, u.dep, fact.imp, x, y, dy);
names(List.res) <- c("pf", "beta", "compt.f", "design.point", "fact.imp", "x", "y", "dy")
if(IS == FALSE){return(List.res)}
}
#------------------------End of the algorithm of Abdo-Rackwitz----------------------#
#---------------------------Importance Sampling :---------------------------------#
if(IS == TRUE){
alpha <- 0.05
N.IS <- N.calls - cp #Number of calls reserved to the importance sampling
norm.u <- u.dep%*%u.dep
P.temp <- 0
P <- 0
VAR <- 0
s <- 0
Z <- numeric(N.IS)
v.temp <- matrix(rnorm(N.IS*ndim,0,1),ncol=ndim)
v.temp <- v.temp + matrix(rep(u.dep,N.IS),ncol=ndim,byrow=TRUE) #One simulate points around u.dep
zz <- v.temp - matrix(rep(u.dep,N.IS), ncol = ndim, byrow = TRUE)
s <- exp( u.dep%*%u.dep/2 - as.numeric(v.temp%*%u.dep) )
test.v.temp <- apply(v.temp, MARGIN = 1, Gtransf)
x <- rbind(x, v.temp)
y <- rbind(y, matrix(test.v.temp,ncol=1))
indic.Def <- 1*(test.v.temp <=0)
P.temp <- indic.Def*s
VAR <- sum(P.temp^2)
cp <- cp + N.IS #Numbers of call to the failure fonction inside the "sapply"
P <- sum(P.temp)/N.IS #failure probability
B <- -qnorm(P) #indices of reliability
VAR <- (sum(P.temp^2)/N.IS - P^2)/(N.IS-1) #standard error of the estimator
IC.inf <- P - qnorm(1-alpha/2)*sqrt(VAR) #low bound of confidence interval
IC.sup <- P + qnorm(1-alpha/2)*sqrt(VAR) #high bound of confidence interval
Inter.conf <- c(IC.inf,IC.sup)
fact.imp <- a^2
List.res <- list(P, B, cp, u.dep, fact.imp, VAR, Inter.conf, x, y, dy)
names(List.res) <-c("pf", "beta", "compt.f", "design.point", "fact.imp", "variance", "conf", "x", "y", "dy");
return(List.res)
}
#----------------------------End Importance Sampling----------------------------#
}
#---------------------------------------------------------------------------------#
#----------------------------- End of algorithm ---------------------------------#
#------------------------------------FORM-----------------------------------------#
#---------------------------------------------------------------------------------#
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.