1 /*****************************************************************************
2 * *
3 * UNURAN -- Universal Non-Uniform Random number generator *
4 * *
5 *****************************************************************************
6 * *
7 * FILE: c_extremeII.c *
8 * *
9 * REFERENCES: *
10 * *
11 * [3] N.L. Johnson, S. Kotz and N. Balakrishnan *
12 * Continuous Univariate Distributions, *
13 * Volume 2, 2nd edition *
14 * John Wiley & Sons, Inc., New York, 1995 *
15 * *
16 *****************************************************************************
17 * *
18 * distr: Extreme value type II distribution [3; ch.22, p.2] *
19 * (also Frechet-type distribution) *
20 * *
21 * *
22 * Type II (also Frechet-type distribution) *
23 * *
24 * cdf: F(x) = exp( -((x-zeta)/theta)^(-k) ) *
25 * pdf: f(x) = exp( -((x-zeta)/theta)^(-k)) * ((x-zeta)/theta)^(-k-1) *
26 * domain: zeta < x <infinity *
27 * constant: k / theta *
28 * *
29 * parameters: 3 *
30 * 0: k > 0 ... shape *
31 * 1: zeta (0) ... location *
32 * 2: theta > 0 (1) ... scale *
33 * *
34 *...........................................................................*
35 * *
36 * Standard Form: *
37 * *
38 * cdf: F(x) = exp( -x^(-k) ) *
39 * pdf: f(x) = exp( -x^(-k) ) * x^(-k-1) *
40 * domain: 0 < x <infinity *
41 * constant: k *
42 * *
43 * parameters: 1 *
44 * 0: k ... shape *
45 * *
46 * 1: zeta = 0 *
47 * 2: theta = 1 *
48 * *
49 *****************************************************************************
50 * *
51 * Copyright (c) 2000-2006 Wolfgang Hoermann and Josef Leydold *
52 * Department of Statistics and Mathematics, WU Wien, Austria *
53 * *
54 * This program is free software; you can redistribute it and/or modify *
55 * it under the terms of the GNU General Public License as published by *
56 * the Free Software Foundation; either version 2 of the License, or *
57 * (at your option) any later version. *
58 * *
59 * This program is distributed in the hope that it will be useful, *
60 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
61 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
62 * GNU General Public License for more details. *
63 * *
64 * You should have received a copy of the GNU General Public License *
65 * along with this program; if not, write to the *
66 * Free Software Foundation, Inc., *
67 * 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA *
68 * *
69 *****************************************************************************/
70
71 /*---------------------------------------------------------------------------*/
72
73 #include <unur_source.h>
74 #include <distr/distr_source.h>
75 #include <distr/cont.h>
76 #include <specfunct/unur_specfunct_source.h>
77 #include "unur_distributions.h"
78 #include "unur_distributions_source.h"
79 #include "unur_stddistr.h"
80
81 /*---------------------------------------------------------------------------*/
82
83 static const char distr_name[] = "extremeII";
84
85 /* parameters */
86 #define k params[0] /* shape */
87 #define zeta params[1] /* location */
88 #define theta params[2] /* scale */
89
90 #define DISTR distr->data.cont
91 #define LOGNORMCONSTANT (distr->data.cont.norm_constant)
92
93 /* function prototypes */
94 static double _unur_pdf_extremeII( double x, const UNUR_DISTR *distr );
95 static double _unur_dpdf_extremeII( double x, const UNUR_DISTR *distr );
96 static double _unur_cdf_extremeII( double x, const UNUR_DISTR *distr );
97 static double _unur_invcdf_extremeII( double u, const UNUR_DISTR *distr );
98
99 static int _unur_upd_mode_extremeII( UNUR_DISTR *distr );
100 static int _unur_upd_area_extremeII( UNUR_DISTR *distr );
101 static int _unur_set_params_extremeII( UNUR_DISTR *distr, const double *params, int n_params );
102
103 /*---------------------------------------------------------------------------*/
104
105 double
_unur_pdf_extremeII(double x,const UNUR_DISTR * distr)106 _unur_pdf_extremeII( double x, const UNUR_DISTR *distr )
107 {
108 register double xk;
109 register const double *params = DISTR.params;
110
111 if (DISTR.n_params > 1)
112 /* standardize */
113 x = (x - zeta) / theta;
114
115 /* standard form */
116
117 if (x<=0.)
118 return 0.;
119
120 xk = pow( x, -k - 1.);
121 return ( exp( -xk * x - LOGNORMCONSTANT) * xk * k );
122
123 } /* end of _unur_pdf_extremeII() */
124
125 /*---------------------------------------------------------------------------*/
126
127 double
_unur_dpdf_extremeII(double x,const UNUR_DISTR * distr)128 _unur_dpdf_extremeII( double x, const UNUR_DISTR *distr )
129 {
130 register double factor = 1.;
131 register double xk;
132 register const double *params = DISTR.params;
133
134 if (DISTR.n_params > 1) {
135 /* standardize */
136 factor = 1. / (theta * theta);
137 x = (x - zeta) / theta;
138 }
139
140 /* standard form */
141
142 if (x<=0.)
143 return 0.;
144
145 xk = pow(x, k);
146 return (- factor * exp(-1./xk) * k * (xk + k*(xk - 1.)) / pow(x,2. + 2.*k)) ;
147
148 } /* end of unur_dpdf_extremeII() */
149
150 /*---------------------------------------------------------------------------*/
151
152 double
_unur_cdf_extremeII(double x,const UNUR_DISTR * distr)153 _unur_cdf_extremeII( double x, const UNUR_DISTR *distr )
154 {
155 register const double *params = DISTR.params;
156
157 if (DISTR.n_params > 1)
158 /* standardize */
159 x = (x - zeta) / theta;
160
161 /* standard form */
162
163 if (x<=0.)
164 return 0.;
165
166 return ( exp( -pow( x, -k ) ) );
167
168 } /* end of _unur_cdf_extremeII() */
169
170 /*---------------------------------------------------------------------------*/
171
172 double
_unur_invcdf_extremeII(double U,const UNUR_DISTR * distr)173 _unur_invcdf_extremeII( double U, const UNUR_DISTR *distr )
174 {
175 register const double *params = DISTR.params;
176 double X;
177
178 X = exp( -log( -log(U) )/k );
179 return ((DISTR.n_params==1) ? X : zeta + theta * X );
180 } /* end of _unur_invcdf_extremeII() */
181
182 /*---------------------------------------------------------------------------*/
183
184 int
_unur_upd_mode_extremeII(UNUR_DISTR * distr)185 _unur_upd_mode_extremeII( UNUR_DISTR *distr )
186 {
187 DISTR.mode = DISTR.zeta + pow( DISTR.k/(DISTR.k+1.), 1/DISTR.k ) * DISTR.theta;
188
189 /* mode must be in domain */
190 if (DISTR.mode < DISTR.domain[0])
191 DISTR.mode = DISTR.domain[0];
192 else if (DISTR.mode > DISTR.domain[1])
193 DISTR.mode = DISTR.domain[1];
194
195 return UNUR_SUCCESS;
196 } /* end of _unur_upd_mode_extremeII() */
197
198 /*---------------------------------------------------------------------------*/
199
200 int
_unur_upd_area_extremeII(UNUR_DISTR * distr)201 _unur_upd_area_extremeII( UNUR_DISTR *distr )
202 {
203 /* log of normalization constant */
204 LOGNORMCONSTANT = log(DISTR.theta);
205
206
207 if (distr->set & UNUR_DISTR_SET_STDDOMAIN) {
208 DISTR.area = 1.;
209 return UNUR_SUCCESS;
210 }
211
212 /* else */
213 DISTR.area = ( _unur_cdf_extremeII( DISTR.domain[1],distr)
214 - _unur_cdf_extremeII( DISTR.domain[0],distr) );
215 return UNUR_SUCCESS;
216
217 } /* end of _unur_upd_area_extremeII() */
218
219 /*---------------------------------------------------------------------------*/
220
221 int
_unur_set_params_extremeII(UNUR_DISTR * distr,const double * params,int n_params)222 _unur_set_params_extremeII( UNUR_DISTR *distr, const double *params, int n_params )
223 {
224
225 /* check number of parameters for distribution */
226 if (n_params < 1) {
227 _unur_error(distr_name,UNUR_ERR_DISTR_NPARAMS,"too few"); return UNUR_ERR_DISTR_NPARAMS; }
228 if (n_params > 3) {
229 _unur_warning(distr_name,UNUR_ERR_DISTR_NPARAMS,"too many");
230 n_params = 3; }
231 CHECK_NULL(params,UNUR_ERR_NULL);
232
233 /* check parameter k */
234 if (k <= 0.) {
235 _unur_error(distr_name,UNUR_ERR_DISTR_DOMAIN,"k <= 0");
236 return UNUR_ERR_DISTR_DOMAIN;
237 }
238
239 /* check parameter theta */
240 if (n_params > 2 && theta <= 0.) {
241 _unur_error(distr_name,UNUR_ERR_DISTR_DOMAIN,"theta <= 0");
242 return UNUR_ERR_DISTR_DOMAIN;
243 }
244
245 /* copy parameters for standard form */
246 DISTR.k = k;
247
248 /* default parameters */
249 DISTR.zeta = 0.;
250 DISTR.theta = 1.;
251
252 /* copy optional parameters */
253 switch (n_params) {
254 case 3:
255 DISTR.theta = theta;
256 case 2:
257 DISTR.zeta = zeta;
258 n_params = 3; /* number of parameters for non-standard form */
259 default:
260 break;
261 }
262
263 /* store number of parameters */
264 DISTR.n_params = n_params;
265
266 /* set (standard) domain */
267 if (distr->set & UNUR_DISTR_SET_STDDOMAIN) {
268 DISTR.domain[0] = DISTR.zeta; /* left boundary */
269 DISTR.domain[1] = INFINITY; /* right boundary */
270 }
271
272 return UNUR_SUCCESS;
273 } /* end of _unur_set_params_extremeII() */
274
275 /*---------------------------------------------------------------------------*/
276
277 struct unur_distr *
unur_distr_extremeII(const double * params,int n_params)278 unur_distr_extremeII( const double *params, int n_params )
279 {
280 register struct unur_distr *distr;
281
282 /* get new (empty) distribution object */
283 distr = unur_distr_cont_new();
284
285 /* set distribution id */
286 distr->id = UNUR_DISTR_EXTREME_II;
287
288 /* name of distribution */
289 distr->name = distr_name;
290
291 /* how to get special generators */
292 /* DISTR.init = _unur_stdgen_extremeII_init; */
293
294 /* functions */
295 DISTR.pdf = _unur_pdf_extremeII; /* pointer to PDF */
296 DISTR.dpdf = _unur_dpdf_extremeII; /* pointer to derivative of PDF */
297 DISTR.cdf = _unur_cdf_extremeII; /* pointer to CDF */
298 DISTR.invcdf = _unur_invcdf_extremeII; /* pointer to inverse CDF */
299
300 /* indicate which parameters are set */
301 distr->set = ( UNUR_DISTR_SET_DOMAIN |
302 UNUR_DISTR_SET_STDDOMAIN |
303 UNUR_DISTR_SET_MODE |
304 UNUR_DISTR_SET_PDFAREA );
305
306 /* set parameters for distribution */
307 if (_unur_set_params_extremeII(distr,params,n_params)!=UNUR_SUCCESS) {
308 free(distr);
309 return NULL;
310 }
311
312 /* log of normalization constant */
313 LOGNORMCONSTANT = log(DISTR.theta);
314
315 /* mode and area below p.d.f. */
316 DISTR.mode = DISTR.zeta + pow( DISTR.k/(DISTR.k+1.), 1/DISTR.k ) * DISTR.theta;
317 DISTR.area = 1.;
318
319 /* function for setting parameters and updating domain */
320 DISTR.set_params = _unur_set_params_extremeII;
321
322 /* function for updating derived parameters */
323 DISTR.upd_mode = _unur_upd_mode_extremeII; /* funct for computing mode */
324 DISTR.upd_area = _unur_upd_area_extremeII; /* funct for computing area */
325
326 /* return pointer to object */
327 return distr;
328
329 } /* end of unur_distr_extremeII() */
330
331 /*---------------------------------------------------------------------------*/
332 #undef c
333 #undef alpha
334 #undef zeta
335 #undef DISTR
336 /*---------------------------------------------------------------------------*/
337