SciRuby/distribution

module Distribution
  # Calculate cdf and inverse cdf for Multivariate Distribution.
  module NormalMultivariate
    class << self
      # Returns multivariate cdf distribution
      # * a is the array of lower values
      # * b is the array of higher values
      # * s is an symmetric positive definite covariance matrix
      def cdf(aa, bb, sigma, epsilon = 0.0001, alpha = 2.5, max_iterations = 100) # :nodoc:
        fail "Doesn't work yet"
        a = [nil] + aa
        b = [nil] + bb
        m = aa.size
        sigma = sigma.to_gsl if sigma.respond_to? :to_gsl

        cc = GSL::Linalg::Cholesky.decomp(sigma)
        c = cc.lower
        intsum = 0
        varsum = 0
        n = 0
        d = Array.new(m + 1, nil)
        e = Array.new(m + 1, nil)
        f = Array.new(m + 1, nil)
        (1..m).each do|i|
          d[i] = 0.0 if a[i].nil?
          e[i] = 1.0 if b[i].nil?
        end
        d[1] = uPhi(a[1].quo(c[0, 0])) unless d[1] == 0
        e[1] = uPhi(b[1].quo(c[0, 0])) unless e[1] == 1
        f[1] = e[1] - d[1]

        error = 1000
        begin
            w = (m + 1).times.collect { |_i| rand * epsilon }
            y = []
            (2..m).each do |i|
              y[i - 1] = iPhi(d[i - 1] + w[i - 1] * (e[i - 1] - d[i - 1]))
              sumc = 0
              (1..(i - 1)).each do |j|
                sumc += c[i - 1, j - 1] * y[j]
              end

              d[i] = uPhi((a[i] - sumc).quo(c[i - 1, i - 1])) unless a[i].nil?
              # puts "sumc:#{sumc}"

              unless b[i].nil?
                # puts "e[#{i}] :#{c[i-1,i-1]}"
                e[i] = uPhi((b[i] - sumc).quo(c[i - 1, i - 1]))
              end
              f[i] = (e[i] - d[i]) * f[i - 1]
            end
            intsum += intsum + f[m]
            varsum += f[m]**2
            n += 1
            error = alpha * Math.sqrt((varsum.quo(n) - (intsum.quo(n))**2).quo(n))
          end while(error > epsilon && n < max_iterations)

        f = intsum.quo(n)
        # p intsum
        # puts "f:#{f}, n:#{n}, error:#{error}"
        f
      end

      def iPhi(pr)
        Distribution::Normal.p_value(pr)
      end

      def uPhi(x)
        Distribution::Normal.cdf(x)
      end
    end
  end
end