1 /*****************************************************************************
2 * *
3 * UNURAN -- Universal Non-Uniform Random number generator *
4 * *
5 *****************************************************************************
6 * *
7 * FILE: c_student_gen.c *
8 * *
9 * Special generators for Student's t 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 student_trouo_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 6 /* maximal number of parameters for generator */
54
55 #define uniform() _unur_call_urng(gen->urng) /* call for uniform prng */
56
57 #define nu (DISTR.params[0]) /* shape (degrees of freedom) */
58
59 /*---------------------------------------------------------------------------*/
60
61 /*****************************************************************************/
62 /** **/
63 /** Inititialize **/
64 /** **/
65 /*****************************************************************************/
66
67 /*---------------------------------------------------------------------------*/
68
69 int
_unur_stdgen_student_init(struct unur_par * par,struct unur_gen * gen)70 _unur_stdgen_student_init( struct unur_par *par, struct unur_gen *gen )
71 /*----------------------------------------------------------------------*/
72 /* initialize special generator for Student's 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: /* Polar Method */
89 _unur_cstd_set_sampling_routine(gen, _unur_stdgen_sample_student_tpol );
90 return UNUR_SUCCESS;
91
92 case 2: /* Ratio of Uniforms */
93 if (par!=NULL && par->distr->data.cont.params[0] < 1.) { /* nu < 1 */
94 _unur_error(NULL,UNUR_ERR_GEN_CONDITION,"");
95 return UNUR_ERR_GEN_CONDITION;
96 }
97 /* nu >= 1 !!!! */
98 _unur_cstd_set_sampling_routine(gen, _unur_stdgen_sample_student_trouo );
99 return student_trouo_init( gen );
100
101 default: /* no such generator */
102 return UNUR_FAILURE;
103 }
104
105 } /* end of _unur_stdgen_chi_init() */
106
107 /*---------------------------------------------------------------------------*/
108
109 /*****************************************************************************/
110 /** **/
111 /** Special generators **/
112 /** **/
113 /*****************************************************************************/
114
115 /*---------------------------------------------------------------------------*/
116
117 /*****************************************************************************
118 * *
119 * Student's t Distribution: Polar Method *
120 * *
121 *****************************************************************************
122 * *
123 * FUNCTION: - samples a random number from Student's t distribution with *
124 * parameters nu > 0. *
125 * *
126 * REFERENCE: - R.W. Bailey (1994): Polar generation of random variates *
127 * with the t-distribution, *
128 * Mathematics of Computation 62, 779-781. *
129 * *
130 * Implemented by F. Niederl, 1994 *
131 *****************************************************************************
132 * *
133 * The polar method of Box/Muller for generating Normal variates is adapted *
134 * to the Student-t distribution. The two generated variates are not *
135 * independent and the expected no. of uniforms per variate is 2.5464. *
136 * *
137 *****************************************************************************
138 * WinRand (c) 1995 Ernst Stadlober, Institut fuer Statistitk, TU Graz *
139 *****************************************************************************/
140
141 double
_unur_stdgen_sample_student_tpol(struct unur_gen * gen)142 _unur_stdgen_sample_student_tpol( struct unur_gen *gen )
143 {
144 /* -X- generator code -X- */
145 double u,v,w;
146
147 /* check arguments */
148 CHECK_NULL(gen,INFINITY);
149 COOKIE_CHECK(gen,CK_CSTD_GEN,INFINITY);
150
151 do {
152 u = 2. * uniform() - 1.;
153 v = 2. * uniform() - 1.;
154 w = u * u + v * v;
155 } while (w > 1.);
156
157 return(u * sqrt( nu * ( exp(- 2. / nu * log(w)) - 1.) / w));
158 /* -X- end of generator code -X- */
159 } /* end of _unur_stdgen_sample_student_tpol() */
160
161 /*---------------------------------------------------------------------------*/
162
163 /*****************************************************************************
164 * *
165 * Student's t Distribution: Ratio of Uniforms *
166 * *
167 *****************************************************************************
168 * *
169 * FUNCTION: - samples a random number from Student's t distribution with *
170 * parameters nu >= 1. *
171 * *
172 * REFERENCE: - A.J. Kinderman, J.F. Monahan (1980): *
173 * New methods for generating Student's t and gamma variables, *
174 * Computing 25, 369-377. *
175 * *
176 * Implemented by R. Kremer, 1990 *
177 *****************************************************************************
178 * WinRand (c) 1995 Ernst Stadlober, Institut fuer Statistitk, TU Graz *
179 *****************************************************************************/
180
181 /*---------------------------------------------------------------------------*/
182 #define c (GEN->gen_param[0])
183 #define e (GEN->gen_param[1])
184 #define p (GEN->gen_param[2])
185 #define q (GEN->gen_param[3])
186 #define r (GEN->gen_param[4])
187 #define vm (GEN->gen_param[5])
188 /*---------------------------------------------------------------------------*/
189
190 inline static int
student_trouo_init(struct unur_gen * gen)191 student_trouo_init( struct unur_gen *gen )
192 {
193 /* check arguments */
194 CHECK_NULL(gen,UNUR_ERR_NULL);
195 COOKIE_CHECK(gen,CK_CSTD_GEN,UNUR_ERR_COOKIE);
196
197 if (GEN->gen_param == NULL) {
198 GEN->n_gen_param = MAX_gen_params;
199 GEN->gen_param = _unur_xmalloc(GEN->n_gen_param * sizeof(double));
200 }
201
202 /* -X- setup code -X- */
203 if (nu < 1.) {
204 _unur_error(NULL,UNUR_ERR_GEN_CONDITION,"");
205 return UNUR_ERR_GEN_CONDITION;
206 }
207
208 r = 1. / nu;
209 p = 1. / (1. + r);
210 q = -0.25 * (nu + 1.);
211 c = 4. * pow(p, q);
212 e = 16. / c;
213 vm = (nu>1.0) ? sqrt(p+p) * pow( (1.-r)*p, 0.25*(nu-1.) ) : 1.;
214 /* -X- end of setup code -X- */
215
216 return UNUR_SUCCESS;
217
218 } /* end of student_trouo_init() */
219
220 double
_unur_stdgen_sample_student_trouo(struct unur_gen * gen)221 _unur_stdgen_sample_student_trouo( struct unur_gen *gen )
222 {
223 /* -X- generator code -X- */
224 double tru,u,v;
225
226 /* check arguments */
227 CHECK_NULL(gen,INFINITY);
228 COOKIE_CHECK(gen,CK_CSTD_GEN,INFINITY);
229
230 while (1) {
231
232 /* step 1 */
233 u = uniform();
234
235 /* step 2 */
236 v = uniform();
237 v = vm * (v + v - 1.);
238 tru = v / u;
239
240 /* step 3 */
241 if ( c * u <= 5. - tru * tru)
242 break;
243 if (nu >= 3.)
244 if (u * (tru * tru + 3.) >= e)
245 continue; /* goto 1 */ /* step 4 */
246 if ( u <= pow(1. + tru * tru * r, q))
247 break;
248 }
249
250 return tru;
251
252 } /* end of _unur_stdgen_sample_student_trouo() */
253
254 /*---------------------------------------------------------------------------*/
255 #undef c
256 #undef e
257 #undef p
258 #undef q
259 #undef r
260 #undef vm
261 /*---------------------------------------------------------------------------*/
262