1 /*****************************************************************************
2 * *
3 * UNURAN -- Universal Non-Uniform Random number generator *
4 * *
5 *****************************************************************************
6 * *
7 * FILE: c_powerexponential_gen.c *
8 * *
9 * Special generators for Power-exponential distribution *
10 * *
11 *****************************************************************************
12 * *
13 * Copyright (c) 2000-2010 Wolfgang Hoermann and Josef Leydold *
14 * Department of Statistics and Mathematics, WU Wien, Austria *
15 * *
16 * This program is free software; you can redistribute it and/or modify *
17 * it under the terms of the GNU General Public License as published by *
18 * the Free Software Foundation; either version 2 of the License, or *
19 * (at your option) any later version. *
20 * *
21 * This program is distributed in the hope that it will be useful, *
22 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
23 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
24 * GNU General Public License for more details. *
25 * *
26 * You should have received a copy of the GNU General Public License *
27 * along with this program; if not, write to the *
28 * Free Software Foundation, Inc., *
29 * 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA *
30 * *
31 *****************************************************************************/
32
33 /*---------------------------------------------------------------------------*/
34
35 #include <unur_source.h>
36 #include <methods/cstd.h>
37 #include <methods/cstd_struct.h>
38 #include <specfunct/unur_specfunct_source.h>
39 #include "unur_distributions_source.h"
40
41 /*---------------------------------------------------------------------------*/
42 /* init routines for special generators */
43
44 inline static int powerexponential_epd_init( struct unur_gen *gen );
45
46 /*---------------------------------------------------------------------------*/
47 /* abbreviations */
48
49 #define PAR ((struct unur_cstd_par*)par->datap) /* data for parameter object */
50 #define GEN ((struct unur_cstd_gen*)gen->datap) /* data for generator object */
51 #define DISTR gen->distr->data.cont /* data for distribution in generator object */
52
53 #define MAX_gen_params 2 /* maximal number of parameters for generator */
54
55 #define uniform() _unur_call_urng(gen->urng) /* call for uniform prng */
56
57 /*---------------------------------------------------------------------------*/
58 #define tau (DISTR.params[0]) /* shape */
59 /*---------------------------------------------------------------------------*/
60
61 /*****************************************************************************/
62 /** **/
63 /** Inititialize **/
64 /** **/
65 /*****************************************************************************/
66
67 /*---------------------------------------------------------------------------*/
68
69 int
_unur_stdgen_powerexponential_init(struct unur_par * par,struct unur_gen * gen)70 _unur_stdgen_powerexponential_init( struct unur_par *par, struct unur_gen *gen )
71 /*----------------------------------------------------------------------*/
72 /* initialize special generator for power-exponential distribution */
73 /* if gen == NULL then only check existance of variant. */
74 /* */
75 /* parameters: */
76 /* par ... pointer to parameter for building generator object */
77 /* gen ... pointer to generator object */
78 /* */
79 /* return: */
80 /* UNUR_SUCCESS ... on success */
81 /* error code ... on error */
82 /*----------------------------------------------------------------------*/
83 {
84 /* one of par and gen must not be the NULL pointer */
85 switch ((par) ? par->variant : gen->variant) {
86
87 case 0: /* DEFAULT */
88 case 1: /* Transformed density rejection */
89 { /* check parameters of distribution */
90 double d_tau = (par) ? par->distr->data.cont.params[0] : tau;
91 if (d_tau < 1.) {
92 _unur_error(NULL,UNUR_ERR_GEN_CONDITION,"");
93 return UNUR_ERR_GEN_CONDITION;
94 }
95 }
96
97 /* tau >= 1 !!!! */
98 _unur_cstd_set_sampling_routine(gen, _unur_stdgen_sample_powerexponential_epd );
99 return powerexponential_epd_init( gen );
100
101 default: /* no such generator */
102 return UNUR_FAILURE;
103 }
104
105 } /* end of _unur_stdgen_powerexponential_init() */
106
107 /*---------------------------------------------------------------------------*/
108
109 /*****************************************************************************/
110 /** **/
111 /** Special generators **/
112 /** **/
113 /*****************************************************************************/
114
115 /*---------------------------------------------------------------------------*/
116
117 /*****************************************************************************
118 * *
119 * Powerexponential Distribution: Transformed density rejection *
120 * *
121 *****************************************************************************
122 * *
123 * FUNCTION: - samples a random number from the Power-exponential *
124 * distribution with parameter tau >= 1 by using the *
125 * non-universal rejection method for logconcave densities. *
126 * *
127 * REFERENCE: - L. Devroye (1986): Non-Uniform Random Variate Generation, *
128 * Springer Verlag, New York. *
129 * *
130 * Implemented by K. Lehner, 1990 *
131 * Revised by F. Niederl, August 1992 *
132 *****************************************************************************
133 * WinRand (c) 1995 Ernst Stadlober, Institut fuer Statistitk, TU Graz *
134 *****************************************************************************/
135
136 /*---------------------------------------------------------------------------*/
137 #define s GEN->gen_param[0]
138 #define sm1 GEN->gen_param[1]
139 /*---------------------------------------------------------------------------*/
140
141 inline static int
powerexponential_epd_init(struct unur_gen * gen)142 powerexponential_epd_init( struct unur_gen *gen )
143 {
144 /* check arguments */
145 CHECK_NULL(gen,UNUR_ERR_NULL);
146 COOKIE_CHECK(gen,CK_CSTD_GEN,UNUR_ERR_COOKIE);
147
148 if (GEN->gen_param == NULL) {
149 GEN->n_gen_param = MAX_gen_params;
150 GEN->gen_param = _unur_xmalloc(GEN->n_gen_param * sizeof(double));
151 }
152
153 /* -X- setup code -X- */
154 s = 1. / tau;
155 sm1 = 1. - s;
156 /* -X- end of setup code -X- */
157
158 return UNUR_SUCCESS;
159
160 } /* end of powerexponential_epd_init() */
161
162 double
_unur_stdgen_sample_powerexponential_epd(struct unur_gen * gen)163 _unur_stdgen_sample_powerexponential_epd( struct unur_gen *gen )
164 {
165 /* -X- generator code -X- */
166 double U,u1,V,X,y;
167
168 /* check arguments */
169 CHECK_NULL(gen,INFINITY);
170 COOKIE_CHECK(gen,CK_CSTD_GEN,INFINITY);
171
172 do {
173 U = 2. * uniform() - 1.; /* U(-1.0/1.0) */
174 u1 = fabs(U); /* u1=|u| */
175 V = uniform(); /* U(0/1) */
176
177 if (u1 <= sm1)
178 /* Uniform hat-function for x <= (1-1/tau) */
179 X = u1;
180 else {
181 /* Exponential hat-function for x > (1-1/tau) */
182 y = tau * (1. - u1); /* U(0/1) */
183 X = sm1 - s * log(y);
184 V *= y;
185 }
186 } while (log(V) > -exp(log(X)*tau)); /* Acceptance/Rejection */
187
188 /* Random sign */
189 if (U > 0.)
190 X = -X;
191
192 /* -X- end of generator code -X- */
193
194 return X;
195
196 } /* end of _unur_stdgen_sample_powerexponential_epd() */
197
198 /*---------------------------------------------------------------------------*/
199 #undef s
200 #undef sm1
201 /*---------------------------------------------------------------------------*/
202