=begin
  I want to noodle around a bit with basic probabilities, odds, that sort of thing.  Here's the example:
  
  Upon being awakened by the sound of a burglar alary, what is your degree of belief that a burglary attempt took place?
  
  Supporting information:
  
  * There is a 95% chance that an attempted burglary will trigger the alarm system, P(alarm|burglary) = 95%
  * There is a 1% chance that the alarms will be triggered by non-burglary attempts, p(alarm|no burglary) = 1%
  * There is a 1/10_000 chance of a particular home being burglarized, generally, P(burglary) = 10^-4
  
  O(burglary|alarm) = L(alarm|burglary)O(burglary)
  
  P(burglary|alarm) = O(burglary|alarm) / 1 + O(burglary|alarm)
  
=end

require 'mathn'

def likelihood(effect, cause)
  effect / cause
end
alias :l :likelihood

def odds(effect, cause=nil)
  cause ? ( likelihood(cause, effect) * odds(effect) ) : (effect / (1 - effect))
end
alias :o :odds

# The same as the odds above, but only in the case of a single hypothesis
def prior_odds(e)
  e / ( 1 - e )
end

def likelihood_ratio(e, h)
  p(e,h) / p(e, 1 - h)
end

def probability(e, h)
  (e * h) / h
end

def product_rule(*e)
end
alias :p :probability

class Array
  
  def givens
    @givens ||= {}
  end
  
  def give(sym, array=Array.new(self.size, 1/self.size))
    self.givens[sym] = array
  end
  
  def given(sym, index)
    self.givens[sym][index]
  end
end

=begin
  Since Rational is the automatic choice for probabilistic data, 
  and since I don't want to override how mathn infers numbers, 
  I am adding some baggage to Rational:
  
  * It still reduces to the LCD
  * It keeps track of all events, 
    so that I can keep a new event proportional to old ones
  * It has an add_event (add) which takes a true or false value
    true values, records that a condition was found
    
  This really only works for binary data, but this is a noodle file.
  
=end
class Rational < Numeric
  class << self
    alias :orig_reduce :reduce
    def reduce(num, den=1)
      val = orig_reduce(num, den)
      val.total_events = den
      val.positive_events = num
      val
    end
  end
  
  def total_events(val=nil)
    @total_events ||= 0
    @total_events = val if val
    @total_events
  end
  alias :events :total_events
  alias :total :total_events
  
  def total_events=(val)
    total_events(val)
  end
  
  def positive_events(val=nil)
    @positive_events ||= 0
    @positive_events = val if val
    @positive_events
  end
  alias :positive :positive_events
  
  def positive_events=(val)
    positive_events(val)
  end
  
  # Uses the Rational constructor to calculate the lowest common denominator
  def add_event(positive=true)
    num = positive ? self.positive_events + 1 : self.positive_events
    den = self.total_events + 1
    other = Rational(num, den)
    @numerator, @denominator, @positive_events, @total_events = other.numerator, other.denominator, num, den
    self
  end
  alias :add :add_event
end


class A
  class << self
    include GSL

    def r
      @r ||= Rng.alloc
    end

    def pdf
      r = self.r.gaussian
      pdf = Ran.gaussian_pdf(r)
      puts r, pdf
    end
  end
end