class Random

# File lib/randomext.rb, line 184
def F(r1, r2)
  f = r2 / r1.to_f
  f*chi_square(r1)/chi_square(r2)
end

# File lib/randomext.rb, line 301
def bernoulli(p)
  (rand < p) ? 1 : 0
end

# File lib/randomext.rb, line 133
def beta(alpha, beta)
  y1 = gamma(alpha); y2 = gamma(beta)
  y1/(y1+y2)
end

# File lib/randomext.rb, line 31
def cauthy(loc, scale)
  loc + scale*standard_cauthy()
end

# File lib/randomext.rb, line 169
def chi_square(r)
  if r == 1
    standard_normal ** 2
  elsif r > 1
    gamma(r*0.5, 2)
  else
    raise ArgumentError, "Random#chi_square:r (degree of freedom) must be >= 1"
  end
end

# File lib/randomext.rb, line 142
def dirichlet(*as)
  if as.any?{|a| a <= 0.0}
    raise ArgumentError, "Random#dirichlet: parameters must be positive"
  end

  ys = as.map{|a| gamma(a) }
  sum = ys.inject(0.0, &:+)
  ys.map{|y| y/sum }
end

# File lib/randomext.rb, line 60
def exponential(scale=1.0)
  if scale < 0.0
    raise ArgumentError, "Random#exponential: scale parameter must be positive"
  end
  scale * standard_exponential
end

# File lib/randomext.rb, line 111
def gamma(shape, scale=1.0)
  if scale <= 0.0
    raise ArgumentError, "Random#gamma: scale parameter must be positive"
  end
  
  case
  when shape <= 0.0
    raise ArgumentError, "Random#gamma: shape parameter must be positive"
  when shape > 1.0
    scale * _gamma(shape)
  when shape == 1.0
    exponential(scale)
  when shape < 1.0
    scale*_gamma(shape+1)*rand_open_interval**(1.0/shape)
  end
end

# File lib/randomext.rb, line 309
def geometric(theta)
  if theta <= 0.0 || theta >= 1.0
    raise ArgumentError, "Random#geometric: theta should be in (0, 1)"
  end
  
  d= -1/(Math.log(1-theta))
  (d * standard_exponential).floor + 1
end

# File lib/randomext.rb, line 102
def gumbel(loc=0.0, scale=1.0)
  loc - scale * Math.log(standard_exponential)
end

# File lib/randomext.rb, line 72
def laplace(loc=0.0, scale=1.0)
  sign = rand(2) == 1 ? 1 : -1
  loc + sign*scale*standard_exponential
end

# File lib/randomext.rb, line 50
def levy(loc=0.0, scale=1.0)
  begin z = standard_normal.abs; end until z > 0
  loc + scale/z**2
end

# File lib/randomext.rb, line 248
def logistic(mu, theta)
  u = rand_open_interval
  mu + theta*log(u/(1-u))
end

# File lib/randomext.rb, line 22
def lognormal(mu=0.0, sigma=1.0)
  Math.exp(normal(mu, sigma))
end

# File lib/randomext.rb, line 352
def logseries(theta)
  if theta <= 0 || 1 <= theta
    raise ArgumentError, "Random#logseries: theta must be in (0, 1)"
  end
  q = 1 - theta
  v = rand_open_interval
  if v >= theta
    1
  else
    u = rand_open_interval
    (log(v)/log(1-q**u)).ceil
  end
end

# File lib/randomext.rb, line 338
def negative_binomial(r, theta)
  if r <= 0.0
    raise ArgumentError, "Random#negative_binomial: r must be positive"
  end
  if theta <= 0.0 && theta >= 1.0
    raise ArgumentError, "Random#negative_binomial: theta must be in (0, 1)"
  end
  poisson(gamma(r, 1/theta - 1))
end

# File lib/randomext.rb, line 258
def non_central_chi_square(r, lambda)
  if lambda < 0.0
    raise ArgumentError, "Random#non_central_chi_square: lambda must be positive"
  end
  if !r.integer? || r <= 0
    raise ArgumentError, "Random#non_central_chi_square: r must be positive integer"
  end
  j = poisson(lambda/2)
  chi_square(r + 2*j)
end

# File lib/randomext.rb, line 274
def non_central_t(r, lambda)
  if lambda == 0.0
    raise ArgumentError, "Random#non_central_t: lambda must not be 0"
  end

  if r == 1
    z = standard_normal + lambda
    w = standard_normal.abs
    z/w
  elsif r == 2
    z = standard_normal + lambda
    w = standard_exponential
    z/Math.sqrt(w)
  elsif r > 2
    d = Math.sqrt(r/2.0)
    z = standard_normal + lambda
    w = _gamma(r/2.0)
    d*z/Math.sqrt(w)
  else
    raise ArgumentError, "Random#non_central_t: r must be positive"
  end
end

# File lib/randomext.rb, line 9
def normal(mean=0.0, sd=1.0)
  mean + standard_normal()*sd
end

# File lib/randomext.rb, line 236
def pareto(a, b=1.0)
  if a <= 0 || b <= 0
    raise ArgumentError, "Random#pareto: parameters a and b must be positive"
  end
  b * (1.0 - rand)**(-1/a)
end

# File lib/randomext.rb, line 323
def planck(a, b)
  if a <= 0 || b <= 0
    raise ArgumentError, "Random#planck: parameters must be positive"
  end
  
  y = _gamma(a+1)
  j = zeta(a+1)
  b*y/j
end

# File lib/randomext.rb, line 158
def power(shape, a, b)
  if shape <= 0 || a >= b
    raise ArgumentError, "Random#power: shape must be positive, and b should be greater than a"
  end
  
  a + (b-a)*(rand_open_interval**(1/shape))
end

# File lib/randomext.rb, line 369
def rand_open_interval
  begin; x = rand; end until x != 0.0
  x
end

# File lib/randomext.rb, line 81
def rayleigh(sigma=1.0)
  sigma*Math.sqrt(2*standard_exponential)
end

# File lib/randomext.rb, line 39
def standard_cauthy
  y1 = standard_normal()
  begin; y2 = standard_normal(); end until y2 != 0.0
  return y1/y2
end

# File lib/randomext.rb, line 193
def t(r)
  if r ==1
    standard_cauthy
  elsif r == 2
    standard_normal/Math.sqrt(exponential(1))
  elsif r > 2
    rdiv2 = r/2.0
    Math.sqrt(rdiv2)*standard_normal/Math.sqrt(_gamma(rdiv2))
  else
    raise ArgumentError, "Random#t: r (degree of freedom) must be >= 1"
  end
end

# File lib/randomext.rb, line 213
def wald(mean, shape)
  if shape <= 0.0
    raise ArgumentError, "Random#wald: shape parameter must be positive"
  end
  p = mean**2
  q = p/(2*shape)
  z = standard_normal
  return mean if z == 0.0
  v = mean + q*z**2
  x1 = v + Math.sqrt(v**2-p)
  return x1 if rand*(x1 + mean) <= mean
  return p/x1
end

# File lib/randomext.rb, line 90
def weibull(g, mu=1.0)
  if g <= 0
    raise ArgumentError, "Random#weibull: shape parameter must be positive"
  end
  mu * standard_exponential**(1.0/g)
end