statistics

  # see: http://rosettacode.org/wiki/Verify_distribution_uniformity/Chi-squared_test#ruby
  def gamma_inc_q(a, x) do
    a1 = a-1
    a2 = a-2

    f0 = fn(t) -> Math.pow(t, a1) * Math.exp(-t) end
    #def f0( t ):
    #    return t**a1*math.exp(-t)
 
    df0 = fn(t) -> (a1-t) * Math.pow(t, a2) * Math.exp(-t) end
    #def df0(t):
    #    return (a1-t)*t**a2*math.exp(-t)
 
    #y = a1
    #while f0(y)*(x-y) >2.0e-8 and y < x: y += .3
    y = gamma_inc_q_red(f0, x, a1)

    if y > x do 
      y = x 
    end
 
    h = 3.0e-4
    n = Math.to_int(y/h)
    h = y/n
    hh = 0.5*h
   
    # gamax = h * sum( f0(t)+hh*df0(t) for t in ( h*j for j in xrange(n-1, -1, -1)))
    gamma_range = for j <- n-1..0, do: h * j
    gamma_series = for t <- gamma_range, do: f0.(t) + hh * df0.(t)
    gammax = h * Enum.sum(gamma_series)
 
    gammax/gamma_spounge(a)
  end

  defp gamma_inc_q_red(f, x, y) do
    if f.(y)*(x-y) > 2.0e-8 and y < x do
      gamma_inc_q_red(f, x, y+0.3)
    else
      y
    end
  end
   

  """
  A = 12
  k1_factrl = 1.0
  coef = [Math.sqrt(2.0*Math.pi)]
  COEF = (1...A).each_with_object(coef) do |k,c|
    c << Math.exp(A-k) * (A-k)**(k-0.5) / k1_factrl
    k1_factrl *= -k
  end
  """
  def gamma_spounge(z) do
    coef = [2.5066282746310002, 198580.06271387744, -696538.0071538023, 984524.6972004091, -719481.3805463575, 290262.7541092609, -64035.01601592932, 7201.864420765038, -354.97463894564885, 5.661005637674728, -0.01474384952133102, 7.490856008760596e-07]
    # accm = (1...A).inject(COEF[0]){|res,k| res += COEF[k] / (z+k)}
    a = z+12
    gamma_spounge_red(coef, z) * Math.exp(-a) * Math.pow(a,(z+0.5)) / z
  end

  defp gamma_spounge_red([h|t], z) do
    gamma_spounge_red(t, z, 1, h)
  end
  defp gamma_spounge_red([], z, k, acc) do
    acc
  end
  defp gamma_spounge_red([h|t], z, k, acc) do
    # accm = (1...A).inject(COEF[0]){|res,k| res += COEF[k] / (z+k)}
    acc = acc + h / (z+k)
    k = k + 1
    gamma_spounge_red(t, z, k, acc)
  end

  """
  def gamma_spounge_py(z) do
    c = None
    #  global c
      a = 12
   
      if c is None:
         k1_factrl = 1.0
         c = []
         c.append(math.sqrt(2.0*math.pi))
         for k in range(1,a):
            c.append( math.exp(a-k) * (a-k)**(k-0.5) / k1_factrl )
            k1_factrl *= -k
   
      accm = c[0]
      for k in range(1,a):
          accm += c[k] / (z+k)
      accm *= math.exp( -(z+a)) * (z+a)**(z+0.5)
      return accm/z;
  end
  """