1 /*****************************************************************************
2 * *
3 * UNURAN -- Universal Non-Uniform Random number generator *
4 * *
5 *****************************************************************************
6 * *
7 * FILE: d_geometric.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: Geometric distribution [1; ch.5.2, p.201] *
19 * *
20 * pmf: p(k) = p * (1-p)^k *
21 * domain: 0 <= k < infinity *
22 * constant: 1 *
23 * *
24 * parameters: *
25 * 0: 0 < p < 1 ... shape *
26 * *
27 *****************************************************************************
28 * *
29 * Copyright (c) 2000-2006 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[] = "geometric";
62
63 /* parameters */
64 #define p params[0]
65
66 #define DISTR distr->data.discr
67 /* #define NORMCONSTANT (distr->data.discr.norm_constant) */
68
69 /*---------------------------------------------------------------------------*/
70
71 /* function prototypes */
72 static double _unur_pmf_geometric( int k, const UNUR_DISTR *distr );
73 static double _unur_cdf_geometric( int k, const UNUR_DISTR *distr );
74 static int _unur_invcdf_geometric( double u, const UNUR_DISTR *distr );
75
76 static int _unur_upd_mode_geometric( UNUR_DISTR *distr );
77 static int _unur_upd_sum_geometric( UNUR_DISTR *distr );
78 static int _unur_set_params_geometric( UNUR_DISTR *distr, const double *params, int n_params );
79
80 /*---------------------------------------------------------------------------*/
81
82 double
_unur_pmf_geometric(int k,const UNUR_DISTR * distr)83 _unur_pmf_geometric(int k, const UNUR_DISTR *distr)
84 {
85 return ((k<0) ? 0. : DISTR.p * pow( 1. - DISTR.p, (double)k ));
86 } /* end of _unur_pmf_geometric() */
87
88 /*---------------------------------------------------------------------------*/
89
90 double
_unur_cdf_geometric(int k,const UNUR_DISTR * distr)91 _unur_cdf_geometric(int k, const UNUR_DISTR *distr)
92 {
93 return ((k<0) ? 0. : (1. - pow(1. - DISTR.p, k+1.)) );
94 } /* end of _unur_cdf_geometric() */
95
96 /*---------------------------------------------------------------------------*/
97
98 int
_unur_invcdf_geometric(double u,const UNUR_DISTR * distr)99 _unur_invcdf_geometric(double u, const UNUR_DISTR *distr)
100 {
101 double x;
102
103 if (_unur_isone(DISTR.p))
104 return 0;
105
106 /* else */
107 x = ceil(log1p(-u) / log1p(-DISTR.p) - 1.);
108
109 return ((x>=INT_MAX) ? INT_MAX : ((int) x));
110 } /* end of _unur_invcdf_geometric() */
111
112 /*---------------------------------------------------------------------------*/
113
114 int
_unur_upd_mode_geometric(UNUR_DISTR * distr)115 _unur_upd_mode_geometric( UNUR_DISTR *distr )
116 {
117 DISTR.mode = 0;
118
119 /* mode must be in domain */
120 if (DISTR.mode < DISTR.domain[0] || DISTR.mode > DISTR.domain[1])
121 DISTR.mode = (DISTR.domain[0]<0) ? 0 : DISTR.domain[0];
122
123 /* o.k. */
124 return UNUR_SUCCESS;
125 } /* end of _unur_upd_mode_geometric() */
126
127 /*---------------------------------------------------------------------------*/
128
129 int
_unur_upd_sum_geometric(UNUR_DISTR * distr)130 _unur_upd_sum_geometric( UNUR_DISTR *distr )
131 {
132 /* normalization constant: none */
133
134 if (distr->set & UNUR_DISTR_SET_STDDOMAIN) {
135 DISTR.sum = 1.;
136 return UNUR_SUCCESS;
137 }
138
139 /* else */
140 DISTR.sum = ( _unur_cdf_geometric( DISTR.domain[1],distr)
141 - _unur_cdf_geometric( DISTR.domain[0]-1,distr) );
142 return UNUR_SUCCESS;
143
144 } /* end of _unur_upd_sum_geometric() */
145
146 /*---------------------------------------------------------------------------*/
147
148 int
_unur_set_params_geometric(UNUR_DISTR * distr,const double * params,int n_params)149 _unur_set_params_geometric( UNUR_DISTR *distr, const double *params, int n_params )
150 {
151 /* check number of parameters for distribution */
152 if (n_params < 1) {
153 _unur_error(distr_name,UNUR_ERR_DISTR_NPARAMS,"too few"); return UNUR_ERR_DISTR_NPARAMS; }
154 if (n_params > 1) {
155 _unur_warning(distr_name,UNUR_ERR_DISTR_NPARAMS,"too many");
156 n_params = 1; }
157 CHECK_NULL(params,UNUR_ERR_NULL);
158
159 /* check parameter p */
160 if (p <= 0. || p >= 1.) {
161 _unur_error(distr_name,UNUR_ERR_DISTR_DOMAIN,"p <= 0 || p >= 1");
162 return UNUR_ERR_DISTR_DOMAIN;
163 }
164
165 /* copy parameters for standard form */
166 DISTR.p = p;
167
168 /* default parameters: none */
169 /* copy optional parameters: none */
170
171 /* store number of parameters */
172 DISTR.n_params = n_params;
173
174 /* set (standard) domain: [0,inifinity] */
175 if (distr->set & UNUR_DISTR_SET_STDDOMAIN) {
176 DISTR.domain[0] = 0; /* left boundary */
177 DISTR.domain[1] = INT_MAX; /* right boundary */
178 }
179
180 return UNUR_SUCCESS;
181 } /* end of _unur_set_params_geometric() */
182
183 /*---------------------------------------------------------------------------*/
184
185 struct unur_distr *
unur_distr_geometric(const double * params,int n_params)186 unur_distr_geometric( const double *params, int n_params )
187 {
188 register struct unur_distr *distr;
189
190 /* get new (empty) distribution object */
191 distr = unur_distr_discr_new();
192
193 /* set distribution id */
194 distr->id = UNUR_DISTR_GEOMETRIC;
195
196 /* name of distribution */
197 distr->name = distr_name;
198
199 /* how to get special generators */
200 /* DISTR.init = _unur_stdgen_geometric_init; */
201
202 /* functions */
203 DISTR.pmf = _unur_pmf_geometric; /* pointer to PMF */
204 DISTR.cdf = _unur_cdf_geometric; /* pointer to CDF */
205 DISTR.invcdf = _unur_invcdf_geometric; /* pointer to invsere of CDF */
206
207 /* indicate which parameters are set */
208 distr->set = ( UNUR_DISTR_SET_DOMAIN |
209 UNUR_DISTR_SET_STDDOMAIN |
210 UNUR_DISTR_SET_MODE |
211 UNUR_DISTR_SET_PMFSUM );
212
213 /* set parameters for distribution */
214 if (_unur_set_params_geometric(distr,params,n_params)!=UNUR_SUCCESS) {
215 free(distr);
216 return NULL;
217 }
218
219 /* log of normalization constant: none */
220
221 /* mode and sum over PMF */
222 DISTR.mode = 0;
223 DISTR.sum = 1.;
224
225 /* function for setting parameters and updating domain */
226 DISTR.set_params = _unur_set_params_geometric;
227
228 /* function for updating derived parameters */
229 DISTR.upd_mode = _unur_upd_mode_geometric; /* funct for computing mode */
230 DISTR.upd_sum = _unur_upd_sum_geometric; /* funct for computing sum */
231
232 /* return pointer to object */
233 return distr;
234
235 } /* end of unur_distr_geometric() */
236
237 /*---------------------------------------------------------------------------*/
238 #undef p
239 #undef DISTR
240 /*---------------------------------------------------------------------------*/
241