Sensitivity analysis for metrology applications

Description

Performs a sensitivity analysis for metrology applications

Usage

sensitivityMet(model, x, y, nboot, method,conf)

Arguments

Details

If method = "SRRC", the function uses src function from the package sensitivity, with the option rank=TRUE to compute SRRC sensitivity indices. This method needs the specification of both the input sample X and either the vector of model responses y either the measurement model as a R function in model.

If method = "Sobol", the function uses sobol2007 function from the package sensitivity. The input sample x is divided into two samples of equal dimensions. Xdesign is returned, containing the design for the computation of Sobol indices. The user should evaluate Xdesign with the computational code and provide the corresponding output values using tell function. Details on the computation of Sobol indices are given in Sobol et al. (2007).

Both methods are applicable whether the measurement model is an explicit function (defined in model as a R function) or an external code.

The argument nboot is required in order to evaluate a confidence interval with a specified confidence level for the sensitivity indices. Default is 0 : in this case no bootstrap replicates and no confidence intervals are computed.

The argument conf defines the confidence level for the bootstrap confidence intervals. Default is 0.95.

Value

sensitivityMet returns a list with the following components :

If method = "SRRC" :

If method = "Sobol" :

Author(s)

Alexandre Allard alexandre.allard@lne.fr

References

I.M. Sobol, S. Tarantola, D. Gatelli, S.S. Kucherenko and W. Mauntz, 2007, Estimating the approximation errors when fixing unessential factors in global sensitivity analysis, Reliability Engineering and System Safety, 92, 957-960.

A. Saltelli, P. Annoni, I. Azzini, F. Campolongo, M. Ratto and S. Tarantola, 2010, Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index, Computer Physics Communications 181, 259-270.

Examples

rm(list=ls())

##Simulate the input sample
M=10000
Xmass <- data.frame(X1 = rnorm(M, 100, 5e-5),
                    X2 = rnorm(M, 0.001234, 2e-5),
                    X3 = runif(M, 1100, 1300),
                    X4 = runif(M,7000000,9000000),
                    X5 = runif(M,7950000,8050000))#Data-frame

#Define the measurement model (GUM-S1, 9.3)
calibMass <-function(x){
  return(((x[,1]+x[,2])*(1+(x[,3]-1200)*(1/x[,4]-1/x[,5]))-100)*1e3)
}

#####   Use SRRC with a model function    #####
#Apply sensitivityMet function to evaluate the associated SRRC indices
S_SRRC=sensitivityMet(model=calibMass,x=Xmass, nboot=100, method="SRRC", conf=0.95)
##Print the results
#First order indices
S_SRRC$S1

#####   Use Sobol with a computational code   #####
#Creation of the design for the computation of Sobol sensitivity indices
S_Sobol=sensitivityMet(model=NULL,x=Xmass,y=NULL, nboot=100, method="Sobol", conf=0.95)

#Obtain the design of experiment to submit to the code
XDesign=S_Sobol$SI$X

#Run the computational code with XDesign as a sample of the input quantities
#We use calibMass function (see GUM-S1) as an example
YDesign=calibMass(XDesign)

#Run the Sobol indices calculations with the outputs of the code
S_Sobol$SI=tell(x=S_Sobol$SI,y=YDesign)

##Print the results
#First order indices
S_Sobol$SI$S
#Total order indices
S_Sobol$SI$T

Results

R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(ATmet)
Loading required package: DiceDesign
Loading required package: lhs
Loading required package: metRology

Attaching package: 'metRology'

The following objects are masked from 'package:base':

    cbind, rbind

Loading required package: msm
Loading required package: sensitivity
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/ATmet/sensitivityMet.Rd_%03d_medium.png", width=480, height=480)
> ### Name: sensitivityMet
> ### Title: Sensitivity analysis for metrology applications
> ### Aliases: sensitivityMet
> ### Keywords: Sobol SRRC metrology
> 
> ### ** Examples
> 
> rm(list=ls())
> 
> ##Simulate the input sample
> M=10000
> Xmass <- data.frame(X1 = rnorm(M, 100, 5e-5),
+                     X2 = rnorm(M, 0.001234, 2e-5),
+                     X3 = runif(M, 1100, 1300),
+                     X4 = runif(M,7000000,9000000),
+                     X5 = runif(M,7950000,8050000))#Data-frame
> 
> #Define the measurement model (GUM-S1, 9.3)
> calibMass <-function(x){
+   return(((x[,1]+x[,2])*(1+(x[,3]-1200)*(1/x[,4]-1/x[,5]))-100)*1e3)
+ }
> 
> #####   Use SRRC with a model function    #####
> #Apply sensitivityMet function to evaluate the associated SRRC indices
> S_SRRC=sensitivityMet(model=calibMass,x=Xmass, nboot=100, method="SRRC", conf=0.95)
> ##Print the results
> #First order indices
> S_SRRC$S1
      original          bias  std. error    min. c.i.   max. c.i.
X1  0.64741215 -5.393327e-04 0.005287171  0.637050709 0.659484931
X2  0.24973074  1.294820e-04 0.008408837  0.234493745 0.266799689
X3  0.04127499  1.975976e-03 0.009719719  0.019222364 0.057583148
X4  0.01154153  1.646541e-03 0.010034573 -0.009592797 0.030634053
X5 -0.01488014 -9.827833e-05 0.006695390 -0.028317870 0.001123308
> 
> #####   Use Sobol with a computational code   #####
> #Creation of the design for the computation of Sobol sensitivity indices
> S_Sobol=sensitivityMet(model=NULL,x=Xmass,y=NULL, nboot=100, method="Sobol", conf=0.95)
> 
> #Obtain the design of experiment to submit to the code
> XDesign=S_Sobol$SI$X
> 
> #Run the computational code with XDesign as a sample of the input quantities
> #We use calibMass function (see GUM-S1) as an example
> YDesign=calibMass(XDesign)
> 
> #Run the Sobol indices calculations with the outputs of the code
> S_Sobol$SI=tell(x=S_Sobol$SI,y=YDesign)
> 
> ##Print the results
> #First order indices
> S_Sobol$SI$S
      original         bias std. error   min. c.i.   max. c.i.
X1  0.08088285 -0.011275601 0.22415268 -0.36948726  0.52361625
X2  0.14306779  0.006409264 0.07753200 -0.03143572  0.29911947
X3 -0.25438425 -0.035563938 0.22492043 -0.72014786  0.20125807
X4 -0.51641446  0.012645286 0.22785716 -1.01185632 -0.04543019
X5  0.01225329  0.001246010 0.01101565 -0.01337042  0.03264698
> #Total order indices
> S_Sobol$SI$T
      original         bias std. error   min. c.i.  max. c.i.
X1  0.80120147  0.012169954 0.22096912  0.36267547 1.23562172
X2 -0.00431594 -0.006331359 0.07724665 -0.15983171 0.17033656
X3  0.77239645  0.031313042 0.22727708  0.34690063 1.26132174
X4  1.00167840 -0.014999963 0.22724290  0.53588698 1.53658004
X5 -0.01347975 -0.001314116 0.01132727 -0.03412102 0.01324717
> 
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>