title: Why are my estimated and theoretical results for sobol sensitivity analysis different?
body: <p>I am working on the sobol sensitivity analysis. I am trying to compute the first order effect and total effect indices in both an estimated and theoretical way.</p>
  
  <p>Firstly, I computed the estimated values by following the steps in Wikipedia "Variance-based sensitivity analysis". </p>
  
  <p><img src="https://i.stack.imgur.com/kVOeY.png" alt="These are the estimated equations"></p>
  
  <p>Here is the code: </p>
  
  <pre><code>set.seed(123)
  
  x_1&lt;-ceiling(1000*runif(1000))  ## generate x1 from range [1,1000]
  x_2&lt;-ceiling(100*runif(1000))   ## generate x2 from range [1,100]
  x_3&lt;-ceiling(10*runif(1000))    ## generate x3 from range [1,10]
  x&lt;-cbind(x_1,x_2,x_3)           ## combine as one matrix
  
  A&lt;-matrix(x[1:500,],ncol=3)     ## divide this one matrix into two 
  B&lt;-matrix(x[501:1000,],ncol=3)  
  
  AB1&lt;-cbind(B[,1],A[,-1])        ## replace the first column of sample A by the first column of sample B
  AB2&lt;-cbind(A[,1],B[,2],A[,3])   ## replace the second column of sample A by the second column of sample B 
  AB3&lt;-cbind(A[,-3],B[,3])        ## replace the third column of sample A by the third column of sample B 
  
  trial&lt;-function(x){             ## define the trial function: 
      #x1+x2*x3^2 
      x[,1]+(x[,2])*(x[,3])^2
  }
  
  Y_A&lt;-trial(A)                 ## the output of A
  Y_B&lt;-trial(B)                 ## the output of B
  
  Y_AB1&lt;-trial(AB1)             ## the output of AB1
  Y_AB2&lt;-trial(AB2)             ## the output of AB2
  Y_AB3&lt;-trial(AB3)             ## the output of AB3
  
  Y&lt;-matrix(cbind(Y_A,Y_B),ncol=1)  ## the matrix of total outputs
  
  S1&lt;-mean(Y_B*(Y_AB1-Y_A))/var(Y)    ## first order effect of x1 
  St1&lt;-(sum((Y_A-Y_AB1)^2)/(2*500))/var(Y) ## total order effect of x1
  
  S2&lt;-mean(Y_B*(Y_AB2-Y_A))/var(Y)     ## first order effect of x2
  St2&lt;-(sum((Y_A-Y_AB2)^2)/(2*500))/var(Y)  ## total order effect of x2
  
  S3&lt;-mean(Y_B*(Y_AB3-Y_A))/var(Y)     ## first order effect of x3
  St3&lt;-(sum((Y_A-Y_AB3)^2)/(2*500))/var(Y) ## total order effect of x3
  
  S&lt;-matrix(c(S1,St1,S2,St2,S3,St3),nrow=2,ncol=3)  ## define the results 
  matrix
  
  rownames(S)&lt;-list("first order","total")
  colnames(S)&lt;-list("X1","X2","X3")
  S                                               ## print result
  </code></pre>
  
  <p>The results are:</p>
  
  <pre><code>                    X1        X2        X3
  first order 0.01734781 0.2337758 0.5261082
  total       0.01523861 0.4078107 0.6471387
  </code></pre>
  
  <p>And then, I wanted to compute the theoretical results to validate the above estimated results. I follow the equations in the picture below:</p>
  
  <p><img src="https://i.stack.imgur.com/yGAeD.png" alt="enter image description here"></p>
  
  <p>In order to compute the partial variance, I approached it in two ways: monte carlo integration and doing the integration by hand ( I just don't trust computer...) </p>
  
  <p>Here is the code:</p>
  
  <pre><code>pimc=function(n){                      ## monte carlo integration
  a=runif(n)
    g=(a+mean(x_2)*(mean(x_3)^2))^2
    pimc=mean(g)
    pimc
  
  }
  
  pimc(10000)/var(Y)
                                                                    ## 
  # doing integration by hand
  (1/3+mean(x_2)*(mean(x_3))^2+((mean(x_2))^2)*((mean(x_3))^4))/var(Y)  
  ## first order effect x1
  (mean(x_1)^2+mean(x_1)*mean(x_3)^2+(1/3)*mean(x_3)^4)/var(Y)          
  ## first order effect x2
  (mean(x_1)^2+(2/3)*(mean(x_1)*mean(x_2))+(1/5)*(mean(x_2)^2))/var(Y)  
  ## first order effect x3
  </code></pre>
  
  <p>The results are the following:</p>
  
  <pre><code>       [,1]
  [1,] 0.4414232        ## first order indice of x1 computed by monte carlo integration
          [,1]
  [1,] 0.4414238          ## first order indice of x1 computed by hand
         [,1]
  [1,] 0.0505044         ## first order indice of x2 computed by hand
         [,1]   
  [1,] 0.05086958         ## first order indice of x3 computed by hand
  </code></pre>
  
  <p>The results are constant by using monte-carlo integration and hand solving. <strong>But they are quite different from the estimated results</strong>. Why is this? </p>
  
  <p>When we do partial integral, we only look at one parameter, and we treat other parameters as constant numbers. So I chose to take the mean of each parameter as the constant number. Is it wrong? </p>