Module: ErrorBootstrapping

Defined in:
lib/ml_ratiosolve/error_bootstrapping.rb

Overview

Methods for using parametric bootstrapping to estimate confidence intervals for the ML ratio estimation.

Author:

  • Colin J. Fuller

Class Method Summary collapse

Class Method Details

.bootstrap_ci(results, level) ⇒ Array

Calculate a bootstrapped confidence interval from the output of #estimate_with_gen_data.

Parameters:

  • results (Array) —

    An array of hashes as output from #estimate_with_gen_data

  • level (Float) —

    A number between 0 and 1 that is the level of the confidence interval. For instance, a value of 0.95 will lead to calculation of the bounds on the central 95% of values.

Returns:

  • (Array) —

    an array of two NMatrix objects, the lower bound on the

    CI and the upper bound on the CI



112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
 
# File 'lib/ml_ratiosolve/error_bootstrapping.rb', line 112

def bootstrap_ci(results, level)
  means = results.map { |e| e[:mu].to_a }
  means = N[*means]
  size_i = means.shape[1]
  size_sim = means.shape[0]
  ci_lower = N.new([size_i], 0.0)
  ci_upper = N.zeros_like(ci_lower)
  size_i.times do |i|
    means_i = means[0...size_sim, i].to_a
    means_i.flatten!
    means_i = means_i.select { |e| e.finite? }
    means_i.sort!
    lower_ci_index = ((1.0-level)/2.0 * means_i.length).ceil
    upper_ci_index = ((1.0 - (1.0-level)/2.0) * means_i.length).floor
    ci_lower[i] = means_i[lower_ci_index]
    ci_upper[i] = means_i[upper_ci_index]
  end
  [ci_lower, ci_upper]
end

.estimate_with_gen_data(n_gen, parameters, x, n_iter, tol = nil) ⇒ Array

Re-estimate distribution parameters by generating simulated data a number of times and performing the iterative estimation in MLRatioSolve

Parameters:

  • n_gen (Numeric) —

    number of datasets to simulate

  • parameters (Hash) —

    A hash containing the parameters from the estimation run on the simulated data.

  • x (NMatrix) —

    The original experimental data

  • n_iter (Numeric) —

    max number of iterations per estimate

  • tol=nil (Numeric) —

    if non-nil, the iterations will terminate early if the absolute change in the likelihood between two successive iterations is less than this

Returns:

  • (Array) —

    An array containing n_gen hashes, each of which is the result of the estimation run on one simulated dataset.



91
92
93
94
95
96
97
98
 
# File 'lib/ml_ratiosolve/error_bootstrapping.rb', line 91

def estimate_with_gen_data(n_gen, parameters, x, n_iter, tol=nil)
  result = Array.new(n_gen) { nil }
  n_gen.times do |i|
    xi = gen_data(parameters, x)
    result[i] = MLRatioSolve.do_iters_with_start(n_iter, xi, parameters[:gamma], tol)
  end
  result
end

.gen_data(parameters, x) ⇒ NMatrix

Generate a set of simulated data consisting of random numbers drawn from the distributions with the supplied parameters

  1. Any skipped data points (as returned by MLRatioSolve::skip_indices)

will be set to 0 here.

Parameters:

  • parameters (Hash) —

    A hash containing the mean, variance, and scale parameters formatted like the output from MLRatioSolve::do_iters_with_start

  • x (NMatrix) —

    The original experimental data

Returns:

  • (NMatrix) —

    A matrix of simulated data with the same dimensions as



58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
 
# File 'lib/ml_ratiosolve/error_bootstrapping.rb', line 58

def gen_data(parameters, x)
  mu = parameters[:mu]
  sig2 = parameters[:sig2]
  gamma = parameters[:gamma]
  size_i = x.shape[0]
  size_n = x.shape[1]

  sim_data_out = N.zeros_like(x)

  size_i.times do |i|
    size_n.times do |n|
      next if MLRatioSolve.skip_indices.include?([i,n])
      sim_data_out[i,n] = randnorm(mu[i]/gamma[n], sig2[i]/gamma[n]**2)
    end
  end

  sim_data_out
end

.randnorm(mu, s2) ⇒ Float

Generate a random normal variate using the Box-Muller transform

Parameters:

  • mu (Numeric) —

    The desired mean

  • s2 (Numeric) —

    The desired variance

Returns:

  • (Float) —

    A random variate from the normal distribution with supplied parameters



40
41
42
43
44
 
# File 'lib/ml_ratiosolve/error_bootstrapping.rb', line 40

def randnorm(mu, s2)
  ru0 = rand
  ru1 = rand
  ([Math.sqrt(-2.0 * Math.log(ru0)) * Math.cos(2*Math::PI*ru1), Math.sqrt(-2.0 * Math.log(ru0)) * Math.sin(2*Math::PI*ru1)].sample)*Math.sqrt(s2) + mu
end