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