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