1 /*****************************************************************************
2 * *
3 * UNURAN -- Universal Non-Uniform Random number generator *
4 * *
5 *****************************************************************************
6 * *
7 * FILE: d_negativebinomial.c *
8 * *
9 * REFERENCES: *
10 * *
11 * [1] N.L. Johnson, S. Kotz and A.W. Kemp *
12 * Univariate Discrete Distributions, *
13 * 2nd edition *
14 * John Wiley & Sons, Inc., New York, 1992 *
15 * *
16 *****************************************************************************
17 * *
18 * distr: Negative Binomial distribution [1; ch.5.1, p.200] *
19 * *
20 * pmf: p(k) = (k+r-1 \choose r-1) * p^r * (1-p)^k *
21 * domain: 0 <= k < infinity *
22 * constant: 1 *
23 * *
24 * parameters: *
25 * 0: 0 < p < 1 *
26 * 1: r > 0 *
27 * *
28 *****************************************************************************
29 * *
30 * Copyright (c) 2000-2010 Wolfgang Hoermann and Josef Leydold *
31 * Department of Statistics and Mathematics, WU Wien, Austria *
32 * *
33 * This program is free software; you can redistribute it and/or modify *
34 * it under the terms of the GNU General Public License as published by *
35 * the Free Software Foundation; either version 2 of the License, or *
36 * (at your option) any later version. *
37 * *
38 * This program is distributed in the hope that it will be useful, *
39 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
40 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
41 * GNU General Public License for more details. *
42 * *
43 * You should have received a copy of the GNU General Public License *
44 * along with this program; if not, write to the *
45 * Free Software Foundation, Inc., *
46 * 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA *
47 * *
48 *****************************************************************************/
49
50 /*---------------------------------------------------------------------------*/
51
52 #include <unur_source.h>
53 #include <distr/distr_source.h>
54 #include <distr/discr.h>
55 #include <specfunct/unur_specfunct_source.h>
56 #include "unur_distributions.h"
57 #include "unur_distributions_source.h"
58 #include "unur_stddistr.h"
59
60 /*---------------------------------------------------------------------------*/
61
62 static const char distr_name[] = "negativebinomial";
63
64 /*---------------------------------------------------------------------------*/
65 /* parameters */
66 #define p params[0]
67 #define r params[1]
68
69 /*---------------------------------------------------------------------------*/
70
71 #define DISTR distr->data.discr
72 #define LOGNORMCONSTANT (distr->data.discr.norm_constant)
73
74 /*---------------------------------------------------------------------------*/
75 /* function prototypes */
76
77 static double _unur_pmf_negativebinomial( int k, const UNUR_DISTR *distr );
78 #ifdef _unur_SF_cdf_negativebinomial
79 static double _unur_cdf_negativebinomial( int k, const UNUR_DISTR *distr );
80 #endif
81 #ifdef _unur_SF_invcdf_negativebinomial
82 static int _unur_invcdf_negativebinomial( double u, const UNUR_DISTR *distr );
83 #endif
84
85 static int _unur_upd_mode_negativebinomial( UNUR_DISTR *distr );
86 static int _unur_upd_sum_negativebinomial( UNUR_DISTR *distr );
87 static int _unur_set_params_negativebinomial( UNUR_DISTR *distr, const double *params, int n_params );
88
89 /*---------------------------------------------------------------------------*/
90
91 double
_unur_pmf_negativebinomial(int k,const UNUR_DISTR * distr)92 _unur_pmf_negativebinomial(int k, const UNUR_DISTR *distr)
93 {
94 const double *params = DISTR.params;
95
96 if (k<0)
97 return 0.;
98
99 else
100 return exp( k*log(1-p)
101 + _unur_SF_ln_gamma(k+r) - _unur_SF_ln_gamma(k+1.) - LOGNORMCONSTANT ) ;
102
103 } /* end of _unur_pmf_negativebinomial() */
104
105 /*---------------------------------------------------------------------------*/
106 #ifdef _unur_SF_cdf_negativebinomial
107
108 double
_unur_cdf_negativebinomial(int k,const UNUR_DISTR * distr)109 _unur_cdf_negativebinomial(int k, const UNUR_DISTR *distr)
110 {
111 const double *params = DISTR.params;
112
113 if (k<0)
114 return 0.;
115
116 else
117 return _unur_SF_cdf_negativebinomial(k,r,p);
118
119 } /* end of _unur_cdf_negativebinomial() */
120
121 #endif
122 /*---------------------------------------------------------------------------*/
123 #ifdef _unur_SF_invcdf_negativebinomial
124
125 int
_unur_invcdf_negativebinomial(double u,const UNUR_DISTR * distr)126 _unur_invcdf_negativebinomial(double u, const UNUR_DISTR *distr)
127 {
128 const double *params = DISTR.params;
129 double x;
130
131 x = _unur_SF_invcdf_negativebinomial(u,r,p);
132 return ((x>=INT_MAX) ? INT_MAX : ((int) x));
133 } /* end of _unur_invcdf_negativebinomial() */
134
135 #endif
136 /*---------------------------------------------------------------------------*/
137
138 int
_unur_upd_mode_negativebinomial(UNUR_DISTR * distr)139 _unur_upd_mode_negativebinomial( UNUR_DISTR *distr )
140 {
141 if (DISTR.r > 1.) {
142 /* mode = floor( (r-1) * (1-p) / p ) */
143 /* (we add a guard against round-off errors */
144 DISTR.mode = (int) ((1.+UNUR_EPSILON) * (DISTR.r - 1.) * (1. - DISTR.p) / DISTR.p);
145 }
146 else { /* r <= 1. */
147 DISTR.mode = 0;
148 }
149
150 /* mode must be in domain */
151 if (DISTR.mode < DISTR.domain[0])
152 DISTR.mode = DISTR.domain[0];
153 else if (DISTR.mode > DISTR.domain[1])
154 DISTR.mode = DISTR.domain[1];
155
156 /* o.k. */
157 return UNUR_SUCCESS;
158 } /* end of _unur_upd_mode_negativebinomial() */
159
160 /*---------------------------------------------------------------------------*/
161
162 int
_unur_upd_sum_negativebinomial(UNUR_DISTR * distr)163 _unur_upd_sum_negativebinomial( UNUR_DISTR *distr )
164 {
165 /* log of normalization constant */
166 LOGNORMCONSTANT = - DISTR.r * log(DISTR.p) + _unur_SF_ln_gamma(DISTR.r);
167
168 if (distr->set & UNUR_DISTR_SET_STDDOMAIN) {
169 DISTR.sum = 1.;
170 return UNUR_SUCCESS;
171 }
172
173 #ifdef _unur_SF_cdf_negativebinomial
174 /* else */
175 DISTR.sum = ( _unur_cdf_negativebinomial( DISTR.domain[1],distr)
176 - _unur_cdf_negativebinomial( DISTR.domain[0]-1,distr) );
177 return UNUR_SUCCESS;
178 #else
179 return UNUR_ERR_DISTR_REQUIRED;
180 #endif
181
182 } /* end of _unur_upd_sum_negativebinomial() */
183
184 /*---------------------------------------------------------------------------*/
185
186 int
_unur_set_params_negativebinomial(UNUR_DISTR * distr,const double * params,int n_params)187 _unur_set_params_negativebinomial( UNUR_DISTR *distr, const double *params, int n_params )
188 {
189 /* check number of parameters for distribution */
190 if (n_params < 2) {
191 _unur_error(distr_name,UNUR_ERR_DISTR_NPARAMS,"too few"); return UNUR_ERR_DISTR_NPARAMS; }
192 if (n_params > 2) {
193 _unur_warning(distr_name,UNUR_ERR_DISTR_NPARAMS,"too many");
194 n_params = 2; }
195 CHECK_NULL(params,UNUR_ERR_NULL);
196
197 /* check parameters */
198 if (p <= 0. || p >= 1. || r <= 0.) {
199 _unur_error(distr_name,UNUR_ERR_DISTR_DOMAIN,"p <= 0 || p >= 1 || r <= 0");
200 return UNUR_ERR_DISTR_DOMAIN;
201 }
202
203 /* copy parameters for standard form */
204 DISTR.p = p;
205 DISTR.r = r;
206
207 /* default parameters: none */
208 /* copy optional parameters: none */
209
210 /* store number of parameters */
211 DISTR.n_params = n_params;
212
213 /* set (standard) domain: [0, infinity] */
214 if (distr->set & UNUR_DISTR_SET_STDDOMAIN) {
215 DISTR.domain[0] = 0; /* left boundary */
216 DISTR.domain[1] = INT_MAX; /* right boundary */
217 }
218
219 return UNUR_SUCCESS;
220 } /* end of _unur_set_params_negativebinomial() */
221
222 /*---------------------------------------------------------------------------*/
223
224 struct unur_distr *
unur_distr_negativebinomial(const double * params,int n_params)225 unur_distr_negativebinomial( const double *params, int n_params )
226 {
227 register struct unur_distr *distr;
228
229 /* get new (empty) distribution object */
230 distr = unur_distr_discr_new();
231
232 /* set distribution id */
233 distr->id = UNUR_DISTR_NEGATIVEBINOMIAL;
234
235 /* name of distribution */
236 distr->name = distr_name;
237
238 /* how to get special generators */
239 /* DISTR.init = _unur_stdgen_negativebinomial_init; */
240
241 /* functions */
242 DISTR.pmf = _unur_pmf_negativebinomial; /* pointer to PMF */
243 #ifdef _unur_SF_cdf_negativebinomial
244 DISTR.cdf = _unur_cdf_negativebinomial; /* pointer to CDF */
245 #endif
246 #ifdef _unur_SF_invcdf_negativebinomial
247 DISTR.invcdf = _unur_invcdf_negativebinomial; /* pointer to inverse CDF */
248 #endif
249
250 /* indicate which parameters are set */
251 distr->set = ( UNUR_DISTR_SET_DOMAIN |
252 UNUR_DISTR_SET_STDDOMAIN |
253 UNUR_DISTR_SET_PMFSUM |
254 UNUR_DISTR_SET_MODE );
255
256 /* set parameters for distribution */
257 if (_unur_set_params_negativebinomial(distr,params,n_params)!=UNUR_SUCCESS) {
258 free(distr);
259 return NULL;
260 }
261
262 /* log of normalization constant */
263 _unur_upd_sum_negativebinomial( distr );
264
265 /* mode and sum over PMF */
266 _unur_upd_mode_negativebinomial(distr);
267 DISTR.sum = 1.;
268
269 /* function for setting parameters and updating domain */
270 DISTR.set_params = _unur_set_params_negativebinomial;
271
272 /* function for updating derived parameters */
273 DISTR.upd_mode = _unur_upd_mode_negativebinomial; /* funct for computing mode */
274 DISTR.upd_sum = _unur_upd_sum_negativebinomial; /* funct for computing area */
275
276 /* return pointer to object */
277 return distr;
278
279 } /* end of unur_distr_negativebinomial() */
280
281 /*---------------------------------------------------------------------------*/
282 #undef p
283 #undef r
284 #undef DISTR
285 /*---------------------------------------------------------------------------*/
286