1 /*****************************************************************************
2 * *
3 * UNURAN -- Universal Non-Uniform Random number generator *
4 * *
5 *****************************************************************************
6 * *
7 * FILE: c_pareto.c *
8 * *
9 * REFERENCES: *
10 * *
11 * [2] N.L. Johnson, S. Kotz and N. Balakrishnan *
12 * Continuous Univariate Distributions, *
13 * Volume 1, 2nd edition *
14 * John Wiley & Sons, Inc., New York, 1994 *
15 * *
16 *****************************************************************************
17 * *
18 * distr: Pareto distribution (of first kind) [2; ch.20, p.574] *
19 * *
20 * pdf: f(x) = x^(-(a+1)) *
21 * domain: x >= k *
22 * constant: a * k^a *
23 * *
24 * parameters: *
25 * 0: k > 0 ... location, shape *
26 * 1: a > 0 ... shape *
27 * *
28 *****************************************************************************
29 * *
30 * Copyright (c) 2000-2006 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/cont.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 static const char distr_name[] = "pareto";
62
63 /* parameters */
64 #define k params[0]
65 #define a params[1]
66
67 #define DISTR distr->data.cont
68 /* #define NORMCONSTANT (distr->data.cont.norm_constant) */
69
70 /* function prototypes */
71 static double _unur_pdf_pareto( double x, const UNUR_DISTR *distr );
72 static double _unur_dpdf_pareto( double x, const UNUR_DISTR *distr );
73 static double _unur_cdf_pareto( double x, const UNUR_DISTR *distr );
74 static double _unur_invcdf_pareto( double u, const UNUR_DISTR *distr );
75
76 static int _unur_upd_mode_pareto( UNUR_DISTR *distr );
77 static int _unur_upd_area_pareto( UNUR_DISTR *distr );
78 static int _unur_set_params_pareto( UNUR_DISTR *distr, const double *params, int n_params );
79
80 /*---------------------------------------------------------------------------*/
81
82 double
_unur_pdf_pareto(double x,const UNUR_DISTR * distr)83 _unur_pdf_pareto( double x, const UNUR_DISTR *distr )
84 {
85 register const double *params = DISTR.params;
86 if (x<k)
87 return 0.;
88 /* else */
89 return ((a/k) / pow(x/k, a + 1.) );
90
91 /* return ( (x<k) ? 0. : (a/k) / pow(x/k, a + 1.) ); */
92 } /* end of _unur_pdf_pareto() */
93
94 /*---------------------------------------------------------------------------*/
95
96 double
_unur_dpdf_pareto(double x,const UNUR_DISTR * distr)97 _unur_dpdf_pareto( double x, const UNUR_DISTR *distr )
98 {
99 register const double *params = DISTR.params;
100 return ( (x<k) ? 0. : a * (-a-1.) / (k * k) * pow(x/k,-a-2.) );
101 } /* end of _unur_dpdf_pareto() */
102
103 /*---------------------------------------------------------------------------*/
104
105 double
_unur_cdf_pareto(double x,const UNUR_DISTR * distr)106 _unur_cdf_pareto( double x, const UNUR_DISTR *distr )
107 {
108 register const double *params = DISTR.params;
109 return ( (x<k) ? 0. : (1. - pow(k/x,a)) );
110 } /* end of _unur_cdf_pareto() */
111
112 /*---------------------------------------------------------------------------*/
113
114 double
_unur_invcdf_pareto(double U,const UNUR_DISTR * distr)115 _unur_invcdf_pareto( double U, const UNUR_DISTR *distr )
116 {
117 register const double *params = DISTR.params;
118 double X;
119
120 X = pow(1-U, -1/a);
121 X *= k;
122 return X;
123 } /* end of _unur_invcdf_pareto() */
124
125 /*---------------------------------------------------------------------------*/
126
127 int
_unur_upd_mode_pareto(UNUR_DISTR * distr)128 _unur_upd_mode_pareto( UNUR_DISTR *distr )
129 {
130 DISTR.mode = DISTR.k;
131
132 /* mode must be in domain */
133 if (DISTR.mode < DISTR.domain[0])
134 DISTR.mode = DISTR.domain[0];
135 else if (DISTR.mode > DISTR.domain[1])
136 DISTR.mode = DISTR.domain[1];
137
138 return UNUR_SUCCESS;
139 } /* end of _unur_upd_mode_pareto() */
140
141 /*---------------------------------------------------------------------------*/
142
143 int
_unur_upd_area_pareto(UNUR_DISTR * distr)144 _unur_upd_area_pareto( UNUR_DISTR *distr )
145 {
146 /* normalization constant: none */
147
148 if (distr->set & UNUR_DISTR_SET_STDDOMAIN) {
149 DISTR.area = 1.;
150 return UNUR_SUCCESS;
151 }
152
153 /* else */
154 DISTR.area = ( _unur_cdf_pareto( DISTR.domain[1],distr)
155 - _unur_cdf_pareto( DISTR.domain[0],distr) );
156 return UNUR_SUCCESS;
157
158 } /* end of _unur_upd_area_pareto() */
159
160 /*---------------------------------------------------------------------------*/
161
162 int
_unur_set_params_pareto(UNUR_DISTR * distr,const double * params,int n_params)163 _unur_set_params_pareto( UNUR_DISTR *distr, const double *params, int n_params )
164 {
165 /* check number of parameters for distribution */
166 /* check new parameter for generator */
167 if (n_params < 2) {
168 _unur_error(distr_name,UNUR_ERR_DISTR_NPARAMS,"too few"); return UNUR_ERR_DISTR_NPARAMS; }
169 if (n_params > 2) {
170 _unur_warning(distr_name,UNUR_ERR_DISTR_NPARAMS,"too many");
171 n_params = 2; }
172 CHECK_NULL(params,UNUR_ERR_NULL);
173
174 /* check parameter k */
175 if (k <= 0.) {
176 _unur_error(distr_name,UNUR_ERR_DISTR_DOMAIN,"k <= 0");
177 return UNUR_ERR_DISTR_DOMAIN;
178 }
179
180 /* check parameter a */
181 if (a <= 0.) {
182 _unur_error(distr_name,UNUR_ERR_DISTR_DOMAIN,"a <= 0");
183 return UNUR_ERR_DISTR_DOMAIN;
184 }
185
186 /* copy parameters for standard form */
187 DISTR.k = k;
188 DISTR.a = a;
189
190 /* default parameters: none */
191 /* copy optional parameters: none */
192
193 /* store number of parameters */
194 DISTR.n_params = n_params;
195
196 /* set (standard) domain */
197 if (distr->set & UNUR_DISTR_SET_STDDOMAIN) {
198 DISTR.domain[0] = DISTR.k; /* left boundary */
199 DISTR.domain[1] = INFINITY; /* right boundary */
200 }
201
202 return UNUR_SUCCESS;
203 } /* end of _unur_set_params_pareto() */
204
205 /*---------------------------------------------------------------------------*/
206
207 struct unur_distr *
unur_distr_pareto(const double * params,int n_params)208 unur_distr_pareto( const double *params, int n_params )
209 {
210 register struct unur_distr *distr;
211
212 /* get new (empty) distribution object */
213 distr = unur_distr_cont_new();
214
215 /* set distribution id */
216 distr->id = UNUR_DISTR_PARETO;
217
218 /* name of distribution */
219 distr->name = distr_name;
220
221 /* how to get special generators */
222 /* DISTR.init = _unur_stdgen_pareto_init; */
223
224 /* functions */
225 DISTR.pdf = _unur_pdf_pareto; /* pointer to PDF */
226 DISTR.dpdf = _unur_dpdf_pareto; /* pointer to derivative of PDF */
227 DISTR.cdf = _unur_cdf_pareto; /* pointer to CDF */
228 DISTR.invcdf = _unur_invcdf_pareto; /* pointer to inverse CDF */
229
230 /* indicate which parameters are set */
231 distr->set = ( UNUR_DISTR_SET_DOMAIN |
232 UNUR_DISTR_SET_STDDOMAIN |
233 UNUR_DISTR_SET_MODE |
234 UNUR_DISTR_SET_PDFAREA );
235
236 /* set parameters for distribution */
237 if (_unur_set_params_pareto(distr,params,n_params)!=UNUR_SUCCESS) {
238 free(distr);
239 return NULL;
240 }
241
242 /* normalization constant: none */
243
244 /* mode and area below p.d.f. */
245 DISTR.mode = DISTR.k;
246 DISTR.area = 1.;
247
248 /* function for setting parameters and updating domain */
249 DISTR.set_params = _unur_set_params_pareto;
250
251 /* function for updating derived parameters */
252 DISTR.upd_mode = _unur_upd_mode_pareto; /* funct for computing mode */
253 DISTR.upd_area = _unur_upd_area_pareto; /* funct for computing area */
254
255 /* return pointer to object */
256 return distr;
257
258 } /* end of unur_distr_pareto() */
259
260 /*---------------------------------------------------------------------------*/
261 #undef k
262 #undef a
263 #undef DISTR
264 /*---------------------------------------------------------------------------*/
265