module ErrorBootstrapping

Calculate a bootstrapped confidence interval from the output of estimate_with_gen_data. @param results [Array] An array of hashes as output from

#estimate_with_gen_data

@param 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.

@return [Array] an array of two NMatrix objects, the lower bound on the

CI and the upper bound on the CI

# 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

Re-estimate distribution parameters by generating simulated data a number of times and performing the iterative estimation in MLRatioSolve @param n_gen [Numeric] number of datasets to simulate @param parameters [Hash] A hash containing the parameters from the estimation run on the simulated data. @param x [NMatrix] The original experimental data @param n_iter [Numeric] max number of iterations per estimate @param 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

@return [Array] An array containing n_gen hashes, each of which is the

result of the estimation run on one simulated dataset.

# 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

Generate a set of simulated data consisting of random numbers drawn from the distributions with the supplied parameters @param parameters [Hash] A hash containing the mean, variance, and scale

parameters formatted like the output from 
MLRatioSolve::do_iters_with_start

@param x [NMatrix] The original experimental data

@return [NMatrix] A matrix of simulated data with the same dimensions as

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

will be set to 0 here.

# 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

Generate a random normal variate using the Box-Muller transform @param mu [Numeric] The desired mean @param s2 [Numeric] The desired variance

@return [Float] A random variate from the normal distribution with

supplied parameters

# 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