mercredi 5 septembre 2018

How to define a function for a global sensitivity analysis in R?

I want to run my model stochastically and determine the total order sensitivity indices using sobol (sensitivity package in R). In my understanding, I need to randomly sample all parameters and generate a dataframe. However, I struggle to understand how I use the values within the dataframe in the actual model and thereby compute the model response to the different values. Below I wrote a small test function.

Error: 'Error in x$model(x$X, ...) : argument "RC1" is missing, with no default'

library(sensitivity)  
n<- 1000
p<- rnorm(n, mean = 579, sd = 21)     
r <- runif(n, min = 0.03, max = 0.07)   
RC1<- runif(n, min = 500, max = 2500)   
RC2<- runif(n, min = 750, max = 3000)   
RC3<- runif(n, min = 3000, max = 7000)   
C1 <- runif(n, min = 1500, max = 2320)   
C2 <- runif(n, min = 1790, max = 2490)   
C3 <- runif(n, min = 1560, max = 2420)   
m1 <- runif(n, min = 100, max = 200)    
m2 <- runif(n, min = 50, max = 100)      
m3 <- runif(n, min = 50, max = 100)      
fb_a1<- runif(n, min = 8, max = 11)        
fb_a2<- runif(n, min = 6, max = 7) 
fb_a3<- runif(n, min = 3, max = 5)         
fb_y1<- rnorm(n, mean = 5, sd = 1)         
fb_y2<- rnorm(n, mean = 7, sd = 2)        
fb_y3<- rnorm(n, mean = 8, sd = 1)

X1 <- data.frame(p, r, RC1, RC2, RC3, C1, C2, C3, m1, m2, m3, fb_a1, fb_a2, fb_a3, fb_y1, fb_y2, fb_y3, Xf_decline, Xf_cost)
X2 <- data.frame(p, r, RC1, RC2, RC3, C1, C2, C3, m1, m2, m3, fb_a1, fb_a2, fb_a3, fb_y1, fb_y2, fb_y3, Xf_decline, Xf_cost)

# Should the function entail all parameter in brackets or the dataframe?
testfun <- function(p,r,RC1,RC2,RC3,C1,C2,C3,m1,m2,m3,fb_a1, fb_a2, fb_a3){
                    ARC<- vector('numeric')                               
                    Y<- list()                                         
                    Pi<- rep( list(vector('numeric')), 3 )         
                    ## Zip specific parameters
                    RC<- c(RC1, RC2, RC3)  # Error occurs here 
                    C<- c(C1, C2, C3)
                    m<- c(m1, m2, m3)
                    fb_a<- c(fb_a1, fb_a2, fb_a3)
                    fb_y<- c(fb_y1, fb_y2, fb_y3)

                    for(c in 1:4){
                      for(i in 1:3){ 
                        ARC[i]    <- (RC[i] * r) / (1 - (1+r) ** (-m[i])) 
                        Y[[i]]    <- approx(c(1, fb_a[i], m[i]-5, m[i]), 
                                            c(0, fb_y[i], fb_y[i], fb_y[i]), 
                                            method="linear", xout=1:m[i])
                        for(a in 1:m[i]){
                          Pi[[i]][a]<- ((p * Y[[i]]$y[a]) - C[i] - ARC[i])
                        }
                      }
                    }
                    return(Pi) }

Y  <- sobolEff(model = testfun, X1=X1, X2=X2, order=0, nboot = 0, conf=0.95)




Aucun commentaire:

Enregistrer un commentaire