require File.expand_path(File.join(File.dirname(__FILE__), '..', '..', 'fathom'))
class Fathom::Distributions::Gaussian
  extend Fathom::Distributions::SharedMethods
  class << self
    def rng
      @rng ||= GSL::Rng.alloc(GSL::Rng::MT19937_1999, Kernel.rand(100_000))
    end
    
    def rand(sd)
      rng.gaussian(sd)
    end
    
    def inverse_cdf(opts={})
      mean = opts[:mean]
      sd = opts[:sd]
      sd ||= opts[:std]
      sd ||= opts[:standard_deviation]
      lower = opts.fetch(:lower, true)
      lower = false if opts[:upper]
      confidence_interval = opts.fetch(:confidence_interval, 0.05)
      value = lower ? GSL::Cdf.gaussian_Pinv(confidence_interval, sd) : GSL::Cdf.gaussian_Qinv(confidence_interval, sd)
      value + mean
    end
    alias :lower_bound :inverse_cdf

    def upper_bound(opts={})
      inverse_cdf(opts.merge(:lower => false))
    end

    def interval_values(opts={})
      confidence_interval = opts.fetch(:confidence_interval, 0.9)
      bound = (1 - confidence_interval) / 2.0
      [lower_bound(opts.merge(:confidence_interval => bound)), upper_bound(opts.merge(:confidence_interval => bound))]
    end
    
    # If only I had the background to explain what this is....
    # I want to know how many standard deviations are expressed by the confidence interval
    # I can then divide  the range by this number to get the standard deviation
    def standard_deviations_under(confidence_interval)
      GSL::Cdf.gaussian_Qinv((1 - confidence_interval) / 2) * 2
    end
    
    
  end
end