a test code for variance-based sensitivity analysis

test_variance_based_sensitivity.rb

require 'pp'

def f(x_vec)
  a0 = 10.0
  vec = [0.0, 1.0, 2.0, 3.0]
  y = a0 + vec.zip(x_vec).map {|c,x| c * x }.inject(:+)
  y += 2.0 * x_vec[0] * x_vec[3]
  y
end

def gaussian(mean, stddev, rand)
  theta = 2 * Math::PI * rand.call
  rho = Math.sqrt(-2 * Math.log(1 - rand.call))
  scale = stddev * rho
  x = mean + scale * Math.cos(theta)
  y = mean + scale * Math.sin(theta)
  return x, y
end

def make_random_vector
  x1, x2 = gaussian(0.0, 1.0, Kernel.method(:rand))
  x3, x4 = gaussian(0.0, 1.0, Kernel.method(:rand))
  [x1, x2, x3, x4]
end

def make_swapped_vector(a, b, i)
  ab_i = a.dup
  ab_i[i] = b[i]
  ab_i
end

def calc_variance(arr)
  ave = arr.inject(:+) / arr.size
  var = arr.map {|x| (x - ave)**2 }.inject(:+) / arr.size
  var
end

num_sample = 10000
dimension = 4
first_order_index = Array.new(dimension, 0.0)
total_effect_index = Array.new(dimension, 0.0)

f_vec = []
num_sample.times do |j|
  vec_a, vec_b = make_random_vector, make_random_vector

  fa = f(vec_a)
  fb = f(vec_b)
  dimension.times do |i|
    vec_ab_i = make_swapped_vector(vec_a, vec_b, i)
    first_order_index[i] += fb * ( f(vec_ab_i) - fa )
    total_effect_index[i] += ( fa - f(vec_ab_i) )**2
  end

  f_vec += [fa,fb]
end

var = calc_variance(f_vec)
first_order_index.map! {|s| s / num_sample / var }
total_effect_index.map! {|s| s / (2*num_sample) / var }

pp first_order_index, total_effect_index