Class: Random

Inherits:
Object
  • Object
show all
Defined in:
lib/randomext.rb

Defined Under Namespace

Classes: Binomial

Instance Method Summary (collapse)

Instance Method Details

- (Integer) bernoulli(p)

Draw a random sample from a Bernoulli distribution.

Parameters:

  • p (Float) —

    the probability returning 1

Returns:

  • (Integer) —

    a random sample, 0 or 1 



301
302
303
 
# File 'lib/randomext.rb', line 301

def bernoulli(p)
  (rand < p) ? 1 : 0
end

- (Float) beta(alpha, beta)

Draws a random sample from a beta distribution.

Parameters:

  • alpha (Float) —

    a shape parameter (alpha > 0.0)

  • beta (Float) —

    another shape parameter (beta > 0.0)

Returns:

  • (Float) —

    a random sample 



133
134
135
136
 
# File 'lib/randomext.rb', line 133

def beta(alpha, beta)
  y1 = gamma(alpha); y2 = gamma(beta)
  y1/(y1+y2)
end

- (Integer) binomial(n, theta)

Draws a random sample from a binomial distribution

Inverse function method is used.

Parameters:

  • n (Integer) —

    the number of trials (n > 0)

  • theta (Float) —

    success probability (0 < theta < 1)

Returns:

  • (Integer) —

    a random sample in 0..n 



61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
 
# File 'ext/binomial.c', line 61

static VALUE random_binomial_inv(VALUE self, VALUE num, VALUE prob)
{
  int n = NUM2INT(num);
  double theta = NUM2DBL(prob);
  int mode = floor(theta*(n+1));
  int xl = mode;
  int xu = mode+1;
  double pl = binomial_distribution(xl, n, theta);
  double pu = pl*forward_ratio(xl, n, theta);
  double u = rb_random_real(self);

  check_binomial_params(n, theta, "Random#binomial");
  
  for (;xl >=0 || xu <= n;) {
    if (xl >= 0) {
      if (u <= pl)
        return INT2NUM(xl);
      u = u - pl;
      pl *= backward_ratio(xl, n, theta);
      --xl;
    }
    if (xu <= n) {
      if (u <= pu)
        return INT2NUM(xu);
      u = u - pu;
      pu *= forward_ratio(xu, n, theta);
      ++xu;
    }
  }
  
  return INT2FIX(0);
}

- (Float) cauthy(loc, scale)

Draws a random sample from a Cauthy distribution.

Parameters:

  • loc (Float) —

    location parameter

  • scale (Float) —

    scale parameter

Returns:

  • (Float) —

    a random sample 



31
32
33
 
# File 'lib/randomext.rb', line 31

def cauthy(loc, scale)
  loc + scale*standard_cauthy()
end

- (Float) chi_square(r)

Draw a random sample from a chi_square distribution.

Parameters:

  • r (Integer) —

    degree of freedom (r >= 1)

Returns:

  • (Float) —

    a random sample 



169
170
171
172
173
174
175
176
177
 
# 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

- (Array<Float>) dirichlet(*as)

Draws a random sample from a Dirichlet distribution.

Parameters:

  • as (Array<Float>) —

    shape parameters

Returns:

  • (Array<Float>) —

    a random sample in the (K-1)-dimensional simplex (K == as.size) 



142
143
144
145
146
147
148
149
150
 
# 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

- (Float) exponential(scale = 1.0)

Draws a random sample from a exponential distribution.

Inverse function method is used.

Parameters:

  • scale (Float) (defaults to: 1.0) —

    scale parameter (scale > 0)

Returns:

  • (Float) —

    a random sample 



60
61
62
63
64
65
 
# 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

- (Float) F(r1, r2)

Draws a random sample from a F distribution.

Parameters:

  • r1 (Integer) —

    degree of freedom (r1 >= 1)

  • r2 (Integer) —

    degree of freedom (r2 >= 1)

Returns:

  • (Float) —

    a random sample 



184
185
186
187
 
# File 'lib/randomext.rb', line 184

def F(r1, r2)
  f = r2 / r1.to_f
  f*chi_square(r1)/chi_square(r2)
end

- (Float) gamma(shape, scale = 1.0)

Draws a random sample from a gamma distribution

Parameters:

  • shape (Float) —

    shape parameter (shape > 0.0)

  • scale (Float) (defaults to: 1.0) —

    scale parameter (scale > 0.0)

Returns:

  • (Float) —

    a random sample 



111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
 
# 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

- (Integer) geometric(theta)

Draws a random sample from a geometric distribution.

Parameters:

  • theta (Float) —

    the probability of sucess (0 < theta < 1)

Returns:

  • (Integer) —

    a random sample in [1, INFINITY) 



309
310
311
312
313
314
315
316
 
# 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

- (Float) gumbel(loc = 0.0, scale = 1.0)

Draws a random sample from a Gumbel distribution

Parameters:

  • loc (Float) (defaults to: 0.0) —

    location parameter

  • scale (Float) (defaults to: 1.0) —

    scale parameter

Returns:

  • (Float) —

    a random sample 



102
103
104
 
# File 'lib/randomext.rb', line 102

def gumbel(loc=0.0, scale=1.0)
  loc - scale * Math.log(standard_exponential)
end

- (Integer) hypergeometric(N, M, n)

Draws a random sample from a hypergeometric distribution.

Inverse method is used.

Parameters:

  • N (Integer) —

    a population (N >= 0)

  • M (Integer) —

    the number of successes (0 <= M <= N)

  • n (Integer) —

    the number of samples (0 <= n <= N)

Returns:

  • (Integer) —

    a random sample in [max(0, n-(N-M)), min(n, M)] 



34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
 
# File 'ext/hypergeometric.c', line 34

static VALUE random_hypergoemtric_inv(VALUE self, VALUE vN, VALUE vM, VALUE vn)
{
  int N = NUM2INT(vN);
  int M = NUM2INT(vM);
  int n = NUM2INT(vn);
  int ok = (N >= 0) && (M >= 0) && (n >= 0) && (M <= N) && (n <= N);
  int mode, x_min, x_max, xu, xl;
  double pl, pu;
  double u;
  
  if (!ok)
    rb_raise(rb_eArgError,
             "Random#hypergeometric: paramters must be:"
             "(N >= 0) && (M >= 0) && (n >= 0) && (M <= N) && (n <= N)");

  mode = (M+1)*(n+1) / (N+2);
  x_min = MAX2(0, n - (N-M));
  x_max = MIN2(n, M);
  xu = mode;
  pu = hypergeometric_distribution(mode, N, M, n);
  xl = mode-1;
  pl = pu * backward_ratio(mode, N, M, n);
  u = rb_random_real(self);

  for (;x_min <= xl || xu <= x_max;) {
    if (xu <= x_max) {
      if (u <= pu)
        return INT2NUM(xu);
      u -= pu;
      pu *= forward_ratio(xu, N, M, n);
      ++xu;
    }
    if (xl >= x_min) {
      if (u <= pl)
        return INT2NUM(xl);
      u -= pl;
      pl *= backward_ratio(xl, N, M, n);
      --xl;
    }
  }
  
  return INT2NUM(x_min);
}

- (Float) laplace(loc = 0.0, scale = 1.0)

Draws a random sample from a Laplace distribution

Parameters:

  • loc (Float) (defaults to: 0.0) —

    location parameter

  • scale (Float) (defaults to: 1.0) —

    scale parameter

Returns:

  • (Float) —

    a random sample 



72
73
74
75
 
# 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

- (Float) levy(loc = 0.0, scale = 1.0)

Draws a random sample from a Levy distribution.

Parameters:

  • loc (Float) (defaults to: 0.0) —

    location parameter

  • scale (Float) (defaults to: 1.0) —

    scale parameter

Returns:

  • (Float) —

    a random sample 



50
51
52
53
 
# 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

- (Float) logistic(mu, theta)

Draws a random sample from a logistic distribution.

Parameters:

  • mu (Float) —

    the location parameter

  • theta (Float) —

    the scale parameter

Returns:

  • (Float) —

    a random sample 



248
249
250
251
 
# File 'lib/randomext.rb', line 248

def logistic(mu, theta)
  u = rand_open_interval
  mu + theta*log(u/(1-u))
end

- (Float) lognormal(mu = 0.0, sigma = 1.0)

Draws a random sample from a log normal distribution.

The probabilistic mass function of lognormal distribution is defined: 


    1/sqrt(2*PI*sigma**2)*exp(-(log(x)-mu)**2/(2*sigma**2))

Parameters:

  • mu (Float) (defaults to: 0.0) —

    the mean in a normal distribution

  • sigma (Float) (defaults to: 1.0) —

    the standard deviarion in a normal distribution

Returns:

  • (Float) —

    a random sample in (0, INFINITY) 



22
23
24
 
# File 'lib/randomext.rb', line 22

def lognormal(mu=0.0, sigma=1.0)
  Math.exp(normal(mu, sigma))
end

- (Integer) logseries(theta)

Draws a random sample from a log series distribution.

Parameters:

  • theta (Float) —

    a parameter (0 < theta < 1)

Returns:

  • (Integer) —

    a random sample in [0, INFINITY) 



352
353
354
355
356
357
358
359
360
361
362
363
364
 
# 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

- (Integer) negative_binomial(r, theta)

Draws a random sample from a negative binomial distribution.

Parameters:

  • r (Float) —

    s parameter (0 < r)

  • theta (Float) —

    a parameter (0 < theta < 1)

Returns:

  • (Integer) —

    a random sample in [0, INFINITY) 



338
339
340
341
342
343
344
345
346
 
# 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

- (Float) non_central_chi_square(r, lambda)

Draws a random sample from a Non-Central Chi-Square distribution.

Parameters:

  • r (Integer) —

    a parameter (r > 0)

  • lambda (Float) —

    another parameter (lambda > 0)

Returns:

  • (Float) —

    a random sample in (0, INFINITY) 



258
259
260
261
262
263
264
265
266
267
 
# 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

- (Float) non_central_t(r, lambda)

Draws a random sample from a Non-Central t distribution

Parameters:

  • r (Integer) —

    a parameter (r > 0)

  • lambda (Float) —

    another parameter (lambda > 0)

Returns:

  • (Float) —

    a random sample 



274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
 
# 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

- (Float) normal(mean = 0.0, sd = 1.0)

Draws a random sample from a normal(Gaussian) distribution.

Parameters:

  • mean (Float) (defaults to: 0.0) —

    mean

  • sd (Float) (defaults to: 1.0) —

    standard deviarion

Returns:

  • (Float) —

    a random sample 



9
10
11
 
# File 'lib/randomext.rb', line 9

def normal(mean=0.0, sd=1.0)
  mean + standard_normal()*sd
end

- (Float) pareto(a, b = 1.0)

Draws a random sample from a Pareto distribution.

The probabilistic mass function for the distribution is defined as: 


    p(x) = a*b**a/x**(a+1)

Parameters:

  • a (Float) —

    shape parameter (a > 0.0)

  • b (Float) (defaults to: 1.0) —

    scale parameter (b > 0.0)

Returns:

  • (Float) —

    a random sample in [b, INFINITY) 



236
237
238
239
240
241
 
# 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

- (Float) planck(a, b)

Draws a random sample from a Planck distribution.

Parameters:

  • a (Float) —

    shape parameter

  • b (Float) —

    scale parameter

Returns:

  • (Float) —

    a random sample in (0, INFINITY) 



323
324
325
326
327
328
329
330
331
 
# 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

- (Integer) poisson(lambda)

Draws a random sample from a Poisson distribution.

Inverse function method is used.

Parameters:

  • lambda (Float) —

    mean

Returns:

  • (Integer) —

    a random sample in [0, INFINITY) 



29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
 
# File 'ext/poisson.c', line 29

static VALUE random_poisson_inv(VALUE self, VALUE l)
{
  double lambda = NUM2DBL(l);
  int mode, xu, xl;
  double pu, pl, u;
  
  if (lambda <= 0.0)
    rb_raise(rb_eArgError, "Random#poisson: lambda must be positive");

  mode = floor(lambda);
  xu = mode;
  xl = mode - 1;
  pu = poisson_distribution(mode, lambda);
  pl = pu * backward_ratio(xu, lambda);
  u = rb_random_real(self);
  
  for (;;) {
    if (u <= pu)
      return INT2NUM(xu);
    u = u - pu;
    pu *= forward_ratio(xu, lambda);
    ++xu;

    if (xl >= 0) {
      if (u <= pl)
        return INT2NUM(xl);
      u = u - pl;
      pl *= backward_ratio(xl, lambda);
      --xl;
    }
  }
}

- (Float) power(shape, a, b)

Draws a random sample from a power function distribution

Parameters:

  • shape (Float) —

    shape parameter (shape > 0.0)

  • a (Float) —

    lower boundary parameter

  • b (Float) —

    upper boundary parameter (a < b)

Returns:

  • (Float) —

    a random sample 



158
159
160
161
162
163
164
 
# 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

- (Float) rand_open_interval

Draw a sample from the uniform distribution on (0, 1)

Returns:

  • (Float) —

    a random sample in (0, 1) 



369
370
371
372
 
# File 'lib/randomext.rb', line 369

def rand_open_interval
  begin; x = rand; end until x != 0.0
  x
end

- (Float) rayleigh(sigma = 1.0)

Draws a random sample from a Rayleigh distribution

Parameters:

  • sigma (Float) (defaults to: 1.0) —

    scale parameter

Returns:

  • (Float) —

    a random sample 



81
82
83
 
# File 'lib/randomext.rb', line 81

def rayleigh(sigma=1.0)
  sigma*Math.sqrt(2*standard_exponential)
end

- (Float) standard_cauthy

Draws a random sample from the standard Cauthy distribution.

Computed using Polar method from the standard normal distribution.

Returns:

  • (Float) —

    a random sample 



39
40
41
42
43
 
# 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

- (Float) standard_exponential

Draws a random sample from the standard exponential distribution.

Returns:

  • (Float) —

    a random sample 



71
72
73
74
 
# File 'ext/standard_exponential.c', line 71

static VALUE random_standard_exp(VALUE self)
{
  return DBL2NUM(standard_exponential(self));
}

- (Float) standard_normal

Draws a random sample from the standard normal distribution.

Ziggurat method is used for random sampling.

Returns:

  • (Float) —

    a random sample 



86
87
88
89
 
# File 'ext/standard_normal.c', line 86

static VALUE random_standard_normal(VALUE self)
{
  return DBL2NUM(randomext_random_standard_normal(self));
}

- (Float) t(r)

Draws a random sample from a t distribution.

Parameters:

  • r (Integer) —

    degree of freedom (r >= 1)

Returns:

  • (Float) —

    a random sample 



193
194
195
196
197
198
199
200
201
202
203
204
 
# 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

- (Float) vonmises(mu, kappa)

Draws a random sample from a von Mises distribution.

The return value is contained in [-PI, PI].

Parameters:

  • mu (Float) —

    direction parameter (-PI <= mu <= PI)

  • kappa (Float) —

    concentration parameter (kappa > 0)

Returns:

  • (Float) —

    A random sample in [-PI, PI] 



13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
 
# File 'ext/other.c', line 13

static VALUE random_vonmises(VALUE self, VALUE vmu, VALUE vkappa)
{
  double mu = NUM2DBL(vmu);
  double kappa = NUM2DBL(vkappa);
  double s;
  
  if (kappa <= 0)
    rb_raise(rb_eArgError, "Random#vonmises: parameter kappa must be positive");
  
  s = (1 + sqrt(1+4*kappa*kappa))/(2*kappa);
  
  for (;;) {
    double u = rb_random_real(self);
    double z = cos(2*M_PI*u);
    double W = (1-s*z)/(s-z);
    double T = kappa*(s-W);
    double U = rb_random_real(self);
    double V = rb_random_real(self);
    double x, y;
    
    if (V > T*(2-T) && V > T*exp(1-T))
      continue;

    if (U < 0.5)
      y = -acos(W);
    else
      y = acos(W);
    x = y + mu;
    if (x >= M_PI)
      return DBL2NUM(x - M_PI);
    else if (x < -M_PI)
      return DBL2NUM(x + M_PI);
    else
      return DBL2NUM(x);
  }
}

- (Float) wald(mean, shape)

Draws a random sample from a wald distribution.

A wald distribution is also called an inverse Gaussian distribution.

Parameters:

  • mean (Float) —

    mean

  • shape (Float) —

    shape parameter (shape > 0.0)

Returns:

  • (Float) —

    a random sample in (0, INFINITY) 



213
214
215
216
217
218
219
220
221
222
223
224
225
 
# 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

- (Float) weibull(g, mu = 1.0)

Draws a random sample from a Weibull distribution

Parameters:

  • g (Float) —

    shape parameter (g > 0.0)

  • mu (Float) (defaults to: 1.0) —

    scale parameter

Returns:

  • (Float) —

    a random sample 



90
91
92
93
94
95
 
# 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

- (Integer) zeta(s)

Draws a random sample from a zeta distribution.

Parameters:

  • s (Integer) —

    a parameter (s > 0.0)

Returns:

  • (Integer) —

    a random sample in [1, INFINITY) 



100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
 
# File 'ext/other.c', line 100

static VALUE random_zeta(VALUE self, VALUE vs)
{
  double s = NUM2DBL(vs);
  double q = s - 1.0;
  double r = -1.0/q;
  double t = pow(2.0, q);

  for (;;) {
    double u = 1.0 - rb_random_real(self);
    double v = rb_random_real(self);
    int x = floor(pow(u, r));
    double w = pow(1.0 + 1.0/x, q);
    if (v*x*(w-1)/(t-1) <= w/t)
      return INT2NUM(x);
  }
}

- (Integer) zipf_mandelbrot(n, q = 0.0, s = 1.0)

Draws a random sample from a Zipf-Mandelbrot distribution.

In case of q == 0.0, the distribution is called a Zipf distribution.

Parameters:

  • n (Integer) —

    the maximum of return value (n > 0)

  • q (Float) —

    a parameter (q >= 0.0)

  • s (Float) —

    a parameter (s > 0.0)

Returns:

  • (Integer) —

    a random sample in 1..n 



61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
 
# File 'ext/other.c', line 61

static VALUE random_zipf(int argc, VALUE *argv, VALUE self)
{
  VALUE vN, vs, vq;
  int N;
  double s, q;
  double sum;
  int i;
  double u;
  
  rb_scan_args(argc, argv, "12", &vN, &vq, &vs);
  N = NUM2INT(vN);
  s = NIL_P(vs) ? 1.0 : NUM2DBL(vs);
  q = NIL_P(vq) ? 0.0 : NUM2DBL(vq);

  if (N <= 0 || s <= 0 || q < 0)
    rb_raise(rb_eArgError, "Random#zipf_mandelbrot: parameters must be N >0, s > 0, q >= 0");
  
  for (i=1, sum=0; i<=N; ++i)
    sum += 1.0/pow(i+q, s);

  u = rb_random_real(self);
  
  for (i=1; i<=N; ++i) {
    double p = 1.0/pow(i+q, s)/sum;
    if (u <= p)
      return INT2NUM(i);
    u -= p;
  }
  
  return INT2NUM(N);
}