module Statistics3::Base

Constants

B0
B1
B10
B12
B14
B16
B2
B4
B6
B8
LOG_2PI
N
SQ2PI

Public Instance Methods

binX_(n, p, x) click to toggle source

discrete distributions

# File lib/statistics3/base.rb, line 518
def binX_(n, p, x); bindist(n, p, x); end
bin_x(n, p, x) click to toggle source
# File lib/statistics3/base.rb, line 519
def bin_x(n, p, x); bindist(n, 1.0 - p, n - x);  end
bindens(n, p, x) click to toggle source
# File lib/statistics3/base.rb, line 403
def bindens(n, p, x)
  p = p.to_f
  q = 1.0 - p
  combi(n, x) * p**x * q**(n - x)
end
bindist(n, p, x) click to toggle source
# File lib/statistics3/base.rb, line 409
def bindist(n, p, x)
  (0..x).inject(0.0) do |s, k|
    s + bindens(n, p, k)
  end
end
chi2X_(n, x) click to toggle source

Returns the integral of Chi-squared distribution with n degrees of freedom over [0, x].

# File lib/statistics3/base.rb, line 462
def chi2X_(n, x); chi2dist(n, x); end
chi2_x(n, x) click to toggle source

Returns the integral of Chi-squared distribution with n degrees of freedom over [x, Infty).

# File lib/statistics3/base.rb, line 459
def chi2_x(n, x); 1.0 - chi2dist(n, x); end
chi2dens(n, x) click to toggle source
# File lib/statistics3/base.rb, line 148
def chi2dens(n, x)
  if n == 1
    1.0/Math.sqrt(2 * Math::PI * x) * Math::E**(-x/2.0)
  elsif n == 2
    0.5 * Math::E**(-x/2.0)
  else
    n = n.to_f
    n2 = n/2
    x = x.to_f
    1.0 / 2**n2 / gamma(n2) * x**(n2 - 1.0) * Math.exp(-x/2.0)
  end
end
chi2dist(n, x) click to toggle source

Returns the integral of Chi-squared distribution with n degrees of freedom over [0, x].

# File lib/statistics3/base.rb, line 189
def chi2dist(n, x); 1.0 - q_chi2(n, x); end
combi(n, x) click to toggle source
# File lib/statistics3/base.rb, line 397
def combi(n, x)
  raise RangeError if n < 0 || x < 0
  x = n - x if x*2 > n
  perm(n, x) / perm(x, x)
end
fX_(n1, n2, x) click to toggle source

Returns the integral of F-distribution with n1 and n2 degrees of freedom over [0, x].

# File lib/statistics3/base.rb, line 506
def fX_(n1, n2, x); fdist(n1, n2, x); end
f_x(n1, n2, x) click to toggle source

Returns the integral of F-distribution with n1 and n2 degrees of freedom over [x, Infty).

# File lib/statistics3/base.rb, line 503
def f_x(n1, n2, x); 1.0 - fdist(n1, n2, x); end
fdist(n1, n2, f) click to toggle source

Returns the integral of F-distribution with n1 and n2 degrees of freedom over [0, x].

# File lib/statistics3/base.rb, line 378
def fdist(n1, n2, f); 1.0 - q_f(n1, n2, f); end
gamma(x) click to toggle source
# File lib/statistics3/base.rb, line 53
def gamma(x)
  if (x < 0.0)
    return Math::PI / (Math.sin(Math.PI * x) * Math.exp(loggamma(1 - x))) #/
  end
  Math.exp(loggamma(x))
end
loggamma(x) click to toggle source

Gamma function

# File lib/statistics3/base.rb, line 34
def loggamma(x)
  v = 1.0
  while (x < N)
    v *= x
    x += 1.0
  end
  w = 1.0 / (x * x)
  ret = B16 / (16 * 15)
  ret = ret * w + B14 / (14 * 13)
  ret = ret * w + B12 / (12 * 11)
  ret = ret * w + B10 / (10 *  9)
  ret = ret * w + B8  / ( 8 *  7)
  ret = ret * w + B6  / ( 6 *  5)
  ret = ret * w + B4  / ( 4 *  3)
  ret = ret * w + B2  / ( 2 *  1)
  ret = ret / x + 0.5 * LOG_2PI - Math.log(v) - x + (x - 0.5) * Math.log(x)
  ret
end
newton_a(y, ini, epsilon = 1.0e-6, limit = 30) { |prev| ... } click to toggle source

Newton approximation

# File lib/statistics3/base.rb, line 19
def newton_a(y, ini, epsilon = 1.0e-6, limit = 30)
  x = ini
  limit.times do |i|
    prev = x
    f, df = yield(prev)
    x = (y - f)/df + prev
    if (x - prev).abs < epsilon
      return x
    end
  end
  $stderr.puts("Warning(newton approximation): over limit")
  x
end
normal__X_(z) click to toggle source

Returns the integral of normal distribution over [0, x].

# File lib/statistics3/base.rb, line 434
def normal__X_(z); normaldist(z) - 0.5; end
normal___x(z) click to toggle source

Returns the integral of normal distribution over [x, Infty).

# File lib/statistics3/base.rb, line 437
def normal___x(z); 1.0 - normaldist(z); end
normaldist(z) click to toggle source

Returns the integral of normal distribution over (-Infty, x].

# File lib/statistics3/base.rb, line 114
def normaldist(z)
  p_nor(z)
end
normalxXX_(z) click to toggle source

Returns the integral of normal distribution over (-Infty, x].

# File lib/statistics3/base.rb, line 431
def normalxXX_(z); normaldist(z); end
normalx__x(z) click to toggle source

Returns the integral of normal distribution over (-Infty, -x] + [x, Infty).

# File lib/statistics3/base.rb, line 440
def normalx__x(z); 2.0 - normaldist(z) * 2.0; end
p_nor(z) click to toggle source

normal-distribution

(-\infty, z]
# File lib/statistics3/base.rb, line 62
def p_nor(z)
  if z < -12 then return 0.0 end
  if z > 12 then return 1.0 end
  if z == 0.0 then return 0.5 end
  
  if z > 0.0
    e = true
  else
    e = false
    z = -z
  end
  z = z.to_f
  z2 = z * z
  t = q = z * Math.exp(-0.5 * z2) / SQ2PI
  
  3.step(199, 2) do |i|
    prev = q
    t *= z2 / i
    q += t
    if q <= prev
      return(e ? 0.5 + q : 0.5 - q)
    end
  end
  e ? 1.0 : 0.0
end
p_t(df, t) click to toggle source

t-distribution ([1]) (-infty, x]

# File lib/statistics3/base.rb, line 196
def p_t(df, t)
  c2 = df.to_f / (df + t * t);
  s = Math.sqrt(1.0 - c2)
  s = -s if t < 0.0
  p = 0.0;
  i = df % 2 + 2
  while i <= df
    p += s
    s *= (i - 1) * c2 / i
    i += 2
  end
  if df & 1 != 0
    0.5+(p*Math.sqrt(c2)+Math.atan(t/Math.sqrt(df)))/Math::PI
  else
    (1.0 + p) / 2.0
  end
end
pchi2(n, y) click to toggle source

[x, infty) Pr([x, infty)) = y -> x

# File lib/statistics3/base.rb, line 163
  def pchi2(n, y)
    if n == 1
      w = pnorm(1 - y/2) # = pnormal___x(y/2)
      w * w
    elsif n == 2
#      v = (1.0 / y - 1.0) / 33.0
#      newton_a(y, v) {|x| [q_chi2(n, x), -chi2dens(n, x)] }
      -2.0 * Math.log(y)
    else
      eps = 1.0e-5
      v = 0.0
      s = 10.0
      loop do
        v += s
        if s <= eps then break end
        if (qe = q_chi2(n, v) - y) == 0.0 then break end
        if qe < 0.0
          v -= s
          s /= 10.0 #/
        end
      end
      v
    end
  end
pchi2X_(n, y) click to toggle source

Return the P-value of the corresponding integral.

# File lib/statistics3/base.rb, line 470
def pchi2X_(n, y); pchi2dist(n, y); end
pchi2_x(n, y) click to toggle source

Return the P-value of the corresponding integral.

# File lib/statistics3/base.rb, line 467
def pchi2_x(n, y); pchi2dist(n, 1.0 - y); end
pchi2dist(n, y) click to toggle source

Returns the P-value of chi2dist().

# File lib/statistics3/base.rb, line 192
def pchi2dist(n, y); pchi2(n, 1.0 - y); end
perm(n, x = n) click to toggle source

discrete distributions

# File lib/statistics3/base.rb, line 386
def perm(n, x = n)
  raise RangeError if n < 0 || x < 0
  r = 1
  while x >= 1
    r *= n
    n -= 1
    x -= 1
  end
  r
end
pf(q, n1, n2) click to toggle source

[x, infty)

# File lib/statistics3/base.rb, line 328
def pf(q, n1, n2)
  if(q < 0.0 || q > 1.0 || n1 < 1 || n2 < 1)
    $stderr.printf("Error : Illegal parameter in pf()!\n")
    return 0.0
  end
  
  if n1 <= 240 || n2 <= 240
    eps = 1.0e-5
    if(n2 == 1) then eps = 1.0e-4 end
    fw = 0.0
    s = 1000.0
    loop do
      fw += s
      if s <= eps  then return fw end
      if (qe = q_f(n1, n2, fw) - q) == 0.0 then return fw end
      if qe < 0.0
        fw -= s
        s /= 10.0 #/
      end
    end
  end
  
  eps = 1.0e-6
  qn = q
  if q < 0.5 then qn = 1.0 - q
    u = pnorm(qn)
    w1 = 2.0 / n1 / 9.0
    w2 = 2.0 / n2 / 9.0
    w3 = 1.0 - w1
    w4 = 1.0 - w2
    u2 = u * u
    a = w4 * w4 - u2 * w2
    b = -2. * w3 * w4
    c = w3 * w3 - u2 * w1
    d = b * b - 4 * a * c
    if(d < 0.0)
      fw = pfsub(a, b, 0.0)
    else
      if(a.abs > eps)
        fw = pfsub(a, b, d)
      else
        if(b.abs > eps) then return -c / b end
        fw = pfsub(a, b, 0.0)
      end
    end
    fw * fw * fw
  end
end
pfX_(n1, n2, x) click to toggle source

Return the P-value of the corresponding integral.

# File lib/statistics3/base.rb, line 515
def pfX_(n1, n2, x); pfdist(n1, n2, x); end
pf_x(n1, n2, x) click to toggle source

Return the P-value of the corresponding integral.

# File lib/statistics3/base.rb, line 512
def pf_x(n1, n2, x); pfdist(n1, n2, 1.0 - x); end
pfdist(n1, n2, y) click to toggle source

Returns the P-value of fdist().

# File lib/statistics3/base.rb, line 381
def pfdist(n1, n2, y); pf(1.0 - y, n1, n2); end
pfsub(x, y, z) click to toggle source

inverse of F-distribution ([2])

# File lib/statistics3/base.rb, line 323
def pfsub(x, y, z)
  (Math.sqrt(z) - y) / x / 2.0
end
pnorm(qn) click to toggle source

inverse of normal distribution ([2]) Pr( (-infty, x] ) = qn -> x

# File lib/statistics3/base.rb, line 90
def pnorm(qn)
  b = [1.570796288, 0.03706987906, -0.8364353589e-3,
       -0.2250947176e-3, 0.6841218299e-5, 0.5824238515e-5,
       -0.104527497e-5, 0.8360937017e-7, -0.3231081277e-8,
       0.3657763036e-10, 0.6936233982e-12]
  
  if(qn < 0.0 || 1.0 < qn)
    $stderr.printf("Error : qn <= 0 or qn >= 1  in pnorm()!\n")
    return 0.0;
  end
  qn == 0.5 and return 0.0
  
  w1 = qn
  qn > 0.5 and w1 = 1.0 - w1
  w3 = -Math.log(4.0 * w1 * (1.0 - w1))
  w1 = b[0]
  1.upto 10 do |i|
    w1 += b[i] * w3**i;
  end
  qn > 0.5 and return Math.sqrt(w1 * w3)
  -Math.sqrt(w1 * w3)
end
pnormal__X_(y) click to toggle source

Return the P-value of the corresponding integral.

# File lib/statistics3/base.rb, line 448
def pnormal__X_(y); pnormalxXX_(y + 0.5); end
pnormal___x(y) click to toggle source

Return the P-value of the corresponding integral.

# File lib/statistics3/base.rb, line 451
def pnormal___x(y); pnormalxXX_(1.0 - y); end
pnormaldist(y) click to toggle source

Returns the P-value of normaldist(x).

# File lib/statistics3/base.rb, line 119
def pnormaldist(y)
  pnorm(y)
end
pnormalxXX_(z) click to toggle source

Return the P-value of the corresponding integral.

# File lib/statistics3/base.rb, line 445
def pnormalxXX_(z); pnormaldist(z); end
pnormalx__x(y) click to toggle source

Return the P-value of the corresponding integral.

# File lib/statistics3/base.rb, line 454
def pnormalx__x(y); pnormalxXX_(1.0 - y/2.0); end
poissonX_(m, x) click to toggle source
# File lib/statistics3/base.rb, line 521
def poissonX_(m, x); poissondist(m, x); end
poisson_x(m, x) click to toggle source
# File lib/statistics3/base.rb, line 522
def poisson_x(m, x); 1.0 - poissondist(m, x-1); end
poissondens(m, x) click to toggle source
# File lib/statistics3/base.rb, line 415
def poissondens(m, x)
  return 0.0 if x < 0
  m = m.to_f
  m ** x * Math::E ** (-m) / perm(x)
end
poissondist(m, x) click to toggle source
# File lib/statistics3/base.rb, line 421
def poissondist(m, x)
  (0..x).inject(0.0) do |s, k|
    s + poissondens(m, k)
  end
end
pt(q, n) click to toggle source
# File lib/statistics3/base.rb, line 240
def pt(q, n)
  q = q.to_f
  if(q < 1.0e-5 || q > 1.0 || n < 1)
    $stderr.printf("Error : Illigal parameter in pt()!\n")
    return 0.0
  end
  
  if(n <= 5) then return ptsub(q, n) end
  if(q <= 5.0e-3 && n <= 13) then return ptsub(q, n) end
  
  f1 = 4.0 * (f = n.to_f)
  f5 = (f4 = (f3 = (f2 = f * f) * f) * f) * f
  f2 *= 96.0
  f3 *= 384.0
  f4 *= 92160.0
  f5 *= 368640.0
  u = pnormaldist(1.0 - q / 2.0)
  
  w0 = (u2 = u * u) * u
  w1 = w0 * u2
  w2 = w1 * u2
  w3 = w2 * u2
  w4 = w3 * u2
  w = (w0 + u) / f1
  w += (5.0 * w1 + 16.0 * w0 + 3.0 * u) / f2
  w += (3.0 * w2 + 19.0 * w1 + 17.0 * w0 - 15.0 * u) / f3
  w += (79.0 * w3 + 776.0 * w2 + 1482.0 * w1 - 1920.0 * w0 - 9450.0 * u) / f4
  w += (27.0 * w4 + 339.0 * w3 + 930.0 * w2 - 1782.0 * w1 - 765.0 * w0 + 17955.0 * u) / f5
  u + w
end
pt__X_(n, y) click to toggle source

Return the P-value of the corresponding integral.

# File lib/statistics3/base.rb, line 495
def pt__X_(n, y); ptdist(n, 0.5 + y); end
pt___x(n, y) click to toggle source

Return the P-value of the corresponding integral.

# File lib/statistics3/base.rb, line 498
def pt___x(n, y); ptdist(n, 1.0 - y); end
ptdist(n, y) click to toggle source

Returns the P-value of tdist().

# File lib/statistics3/base.rb, line 275
def ptdist(n, y)
  if y > 0.5
    pt(2.0 - y*2.0, n)
  else
    - pt(y*2.0, n)
  end
end
ptsub(q, n) click to toggle source

inverse of t-distribution ([2]) (-infty, -q/2] + [q/2, infty)

# File lib/statistics3/base.rb, line 216
def ptsub(q, n)
  q = q.to_f
  if(n == 1 && 0.001 < q && q < 0.01)
    eps = 1.0e-4
  elsif (n == 2 && q < 0.0001)
    eps = 1.0e-4
  elsif (n == 1 && q < 0.001)
    eps = 1.0e-2
  else
    eps = 1.0e-5
  end
  s = 10000.0
  w = 0.0
  loop do
    w += s
    if(s <= eps) then return w end
    if((qe = 2.0 - p_t(n, w)*2.0 - q) == 0.0) then return w end
    if(qe < 0.0)
      w -= s
      s /= 10.0 #/
    end
  end
end
ptxXX_(n, y) click to toggle source

Return the P-value of the corresponding integral.

# File lib/statistics3/base.rb, line 492
def ptxXX_(n, y); ptdist(n, y); end
ptx__x(n, y) click to toggle source

Return the P-value of the corresponding integral.

# File lib/statistics3/base.rb, line 489
def ptx__x(n, y); ptdist(n, 1.0 - y / 2.0); end
q_chi2(df, chi2) click to toggle source

chi-square distribution ([1]) [x, infty)

# File lib/statistics3/base.rb, line 125
def q_chi2(df, chi2)
  chi2 = chi2.to_f
  if (df & 1) != 0
    chi = Math.sqrt(chi2)
    if (df == 1) then return 2 * normal___x(chi); end
    s = t = chi * Math.exp(-0.5 * chi2) / SQ2PI
    k = 3
    while k < df
      t *= chi2 / k;  s += t;
      k += 2
    end
    2 * (normal___x(chi) + s)
  else
    s = t = Math.exp(-0.5 * chi2)
    k = 2
    while k < df
      t *= chi2 / k;  s += t;
      k += 2
    end
    s
  end
end
q_f(df1, df2, f) click to toggle source

F-distribution ([1]) [x, infty)

# File lib/statistics3/base.rb, line 285
def q_f(df1, df2, f)
  if (f <= 0.0) then return 1.0; end
  if (df1 % 2 != 0 && df2 % 2 == 0)
    return 1.0 - q_f(df2, df1, 1.0 / f)
  end
  cos2 = 1.0 / (1.0 + df1.to_f * f / df2.to_f)
  sin2 = 1.0 - cos2
  
  if (df1 % 2 == 0)
      prob = cos2 ** (df2.to_f / 2.0)
      temp = prob
      i = 2
      while i < df1
        temp *= (df2.to_f + i - 2) * sin2 / i
        prob += temp
        i += 2
      end
      return prob
  end
  prob = Math.atan(Math.sqrt(df2.to_f / (df1.to_f * f)))
  temp = Math.sqrt(sin2 * cos2)
  i = 3
  while i <= df1
    prob += temp
    temp *= (i - 1).to_f * sin2 / i.to_f;
    i += 2.0
  end
  temp *= df1.to_f
  i = 3
  while i <= df2
    prob -= temp
    temp *= (df1.to_f + i - 2) * cos2 / i.to_f
    i += 2
  end
  prob * 2.0 / Math::PI
end
t__X_(n, x) click to toggle source

Returns the integral of t-distribution with n degrees of freedom over [0, x].

# File lib/statistics3/base.rb, line 481
def t__X_(n, x); tdist(n, x) - 0.5; end
t___x(n, x) click to toggle source

Returns the integral of t-distribution with n degrees of freedom over [x, Infty).

# File lib/statistics3/base.rb, line 484
def t___x(n, x); 1.0 - tdist(n, x); end
tdist(n, t) click to toggle source

Returns the integral of t-distribution with n degrees of freedom over (-Infty, x].

# File lib/statistics3/base.rb, line 272
def tdist(n, t); p_t(n, t); end
txXX_(n, x) click to toggle source

Returns the integral of t-distribution with n degrees of freedom over (-Infty, x].

# File lib/statistics3/base.rb, line 478
def txXX_(n, x); tdist(n, x); end
tx__x(n, x) click to toggle source

Returns the integral of normal distribution with n degrees of freedom over (-Infty, -x] + [x, Infty).

# File lib/statistics3/base.rb, line 475
def tx__x(n, x); 2.0 - tdist(n, x) * 2.0; end