1 /*****************************************************************************
2 * *
3 * UNURAN -- Universal Non-Uniform Random number generator *
4 * *
5 *****************************************************************************
6 * *
7 * FILE: vc_multistudent.c *
8 * *
9 *****************************************************************************
10 * *
11 * distr: Multivariate Student distribution [5; ch.45, p.219] *
12 * *
13 * pdf: f(x) = 1 / ( 1 + (x-mu)^t . Sigma^-1 . (x-mu) / m)^(dim+m)/2 )*
14 * domain: Reals^(dim) *
15 * constant: Gamma((dim+m)/2) *
16 * / ( Gamma(m/2) (m*pi)^(dim/2) * sqrt(det(Sigma)) ) *
17 * *
18 * parameters: *
19 * 0: d.f. ... m (default : 1) *
20 * 1: mean ... mu (default : 0-vector) *
21 * 2: "covar" ... Sigma (default : identity matrix) *
22 * *
23 *****************************************************************************
24 * *
25 * standard form *
26 * *
27 * pdf: f(x) = 1 / ( 1 + x^t . x )^(dim+m)/2 *
28 * domain: Reals^(dim) *
29 * *
30 * parameters: *
31 * *
32 * d.f. = m *
33 * mean = (0,...,0) ... 0-vector *
34 * covar = identity matrix *
35 * *
36 *****************************************************************************
37 * *
38 * Copyright (c) 2000-2010 Wolfgang Hoermann and Josef Leydold *
39 * Department of Statistics and Mathematics, WU Wien, Austria *
40 * *
41 * This program is free software; you can redistribute it and/or modify *
42 * it under the terms of the GNU General Public License as published by *
43 * the Free Software Foundation; either version 2 of the License, or *
44 * (at your option) any later version. *
45 * *
46 * This program is distributed in the hope that it will be useful, *
47 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
48 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
49 * GNU General Public License for more details. *
50 * *
51 * You should have received a copy of the GNU General Public License *
52 * along with this program; if not, write to the *
53 * Free Software Foundation, Inc., *
54 * 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA *
55 * *
56 *****************************************************************************/
57
58 /*---------------------------------------------------------------------------*/
59
60 #include <unur_source.h>
61 #include <distr/distr_source.h>
62 #include <distr/distr.h>
63 #include <distr/cvec.h>
64 #include <specfunct/unur_specfunct_source.h>
65 #include <utils/matrix_source.h>
66 #include "unur_distributions.h"
67 #include "unur_distributions_source.h"
68 #include "unur_stddistr.h"
69
70 #ifdef USE_DEPRECATED_CODE
71 # include <distr/deprecated_distr.h>
72 #endif
73
74 /*---------------------------------------------------------------------------*/
75
76 static const char distr_name[] = "multistudent";
77
78 /*---------------------------------------------------------------------------*/
79 /* parameters */
80 #define nu params[0]
81
82 /*---------------------------------------------------------------------------*/
83
84 #define DISTR distr->data.cvec
85 #define LOGNORMCONSTANT (distr->data.cvec.norm_constant)
86
87 /*---------------------------------------------------------------------------*/
88 /* function prototypes */
89
90 static double _unur_pdf_multistudent( const double *x, UNUR_DISTR *distr );
91 static double _unur_logpdf_multistudent( const double *x, UNUR_DISTR *distr );
92 static int _unur_dlogpdf_multistudent( double *result, const double *x, UNUR_DISTR *distr );
93 static double _unur_pdlogpdf_multistudent( const double *x, int coord, UNUR_DISTR *distr );
94
95 static int _unur_set_params_multistudent( UNUR_DISTR *distr, double df );
96 static int _unur_upd_mode_multistudent( UNUR_DISTR *distr );
97 static int _unur_upd_volume_multistudent( UNUR_DISTR *distr );
98
99 /*---------------------------------------------------------------------------*/
100
101 double
_unur_pdf_multistudent(const double * x,UNUR_DISTR * distr)102 _unur_pdf_multistudent( const double *x, UNUR_DISTR *distr )
103 {
104 return exp( _unur_logpdf_multistudent( x, distr ) );
105 } /* end of _unur_pdf_multistudent() */
106
107 /*---------------------------------------------------------------------------*/
108
109 double
_unur_logpdf_multistudent(const double * x,UNUR_DISTR * distr)110 _unur_logpdf_multistudent( const double *x, UNUR_DISTR *distr )
111 {
112 #define idx(a,b) ((a)*dim+(b))
113
114 int i,j, dim;
115 double *mean;
116 const double *covar_inv;
117
118 double xx; /* argument used in the evaluation of exp(-xx/2) */
119 double cx; /* element of multiplication of covariance matrix and x */
120
121 dim = distr->dim;
122
123 if (DISTR.mean == NULL) {
124 if (DISTR.covar != NULL) {
125 _unur_warning(distr->name,UNUR_ERR_SHOULD_NOT_HAPPEN,"");
126 }
127 /* standard form */
128 xx=0.;
129 for (i=0; i<dim; i++) { xx += x[i]*x[i]; }
130 return ( - (dim+DISTR.nu)/2. * log(1+xx/DISTR.nu) + LOGNORMCONSTANT);
131 }
132
133 /* mean vector */
134 mean = DISTR.mean;
135
136 /* get inverse of covariance matrix */
137 covar_inv = unur_distr_cvec_get_covar_inv(distr);
138 if (covar_inv==NULL)
139 /* inverse of covariance matrix not available */
140 return INFINITY;
141
142 xx=0.; /* resetting exponential function argument */
143 for (i=0; i<dim; i++) {
144 cx=0.;
145 /* multiplication of inverse covariance matrix and (x-mean) */
146 for (j=0; j<dim; j++) {
147 cx += covar_inv[idx(i,j)] * (x[j]-mean[j]);
148 }
149 xx += (x[i]-mean[i])*cx;
150 }
151
152 return (- (dim+DISTR.nu)/2. * log(1+xx/DISTR.nu) + LOGNORMCONSTANT);
153
154 #undef idx
155 } /* end of _unur_logpdf_multistudent() */
156
157 /*---------------------------------------------------------------------------*/
158
159 int
_unur_dlogpdf_multistudent(double * result,const double * x,UNUR_DISTR * distr)160 _unur_dlogpdf_multistudent( double *result, const double *x, UNUR_DISTR *distr )
161 {
162 #define idx(a,b) ((a)*dim+(b))
163
164 int i, j, dim;
165 double xx, cx;
166 double *mean;
167 const double *covar_inv;
168
169 dim = distr->dim;
170 mean = DISTR.mean;
171
172 /* get inverse of covariance matrix */
173 covar_inv = unur_distr_cvec_get_covar_inv(distr);
174 if (covar_inv==NULL)
175 /* inverse of covariance matrix not available */
176 return UNUR_FAILURE;
177
178 /* calculating (x-mean)' Sigma^-1 (x-mean) */
179 xx=0.;
180 for (i=0; i<dim; i++) {
181 cx=0.;
182 /* multiplication of inverse covariance matrix and (x-mean) */
183 for (j=0; j<dim; j++) {
184 cx += covar_inv[idx(i,j)] * (x[j]-mean[j]);
185 }
186 xx += (x[i]-mean[i])*cx;
187 }
188
189 /* calculation of the dlogpdf components */
190 for (i=0; i<dim; i++) {
191 result[i] = 0.;
192
193 for (j=0; j<dim; j++)
194 result[i] -= (x[j]-mean[j]) * (covar_inv[idx(i,j)]+covar_inv[idx(j,i)]);
195
196 result[i] *= .5*(dim+DISTR.nu)/(DISTR.nu+xx);
197 }
198
199 return UNUR_SUCCESS;
200
201 #undef idx
202 } /* end of _unur_dlogpdf_multistudent() */
203
204 /*---------------------------------------------------------------------------*/
205
206 double
_unur_pdlogpdf_multistudent(const double * x,int coord,UNUR_DISTR * distr)207 _unur_pdlogpdf_multistudent( const double *x, int coord, UNUR_DISTR *distr )
208 {
209 #define idx(a,b) ((a)*dim+(b))
210
211 int i,j;
212 double xx, cx;
213 const double *covar_inv;
214 double result;
215
216 int dim = distr->dim;
217 double *mean = DISTR.mean;
218
219 /* check arguments */
220 if (coord < 0 || coord >= dim) {
221 _unur_warning(distr->name,UNUR_ERR_DISTR_DOMAIN,"invalid coordinate");
222 return INFINITY;
223 }
224
225 /* get inverse of covariance matrix */
226 covar_inv = unur_distr_cvec_get_covar_inv(distr);
227 if (covar_inv==NULL)
228 /* inverse of covariance matrix not available */
229 return INFINITY;
230
231 /* calculating (x-mean)' Sigma^-1 (x-mean) */
232 xx=0.;
233 for (i=0; i<dim; i++) {
234 cx=0.;
235 /* multiplication of inverse covariance matrix and (x-mean) */
236 for (j=0; j<dim; j++) {
237 cx += covar_inv[idx(i,j)] * (x[j]-mean[j]);
238 }
239 xx += (x[i]-mean[i])*cx;
240 }
241
242 /* calculation of pdlogpdf of given component */
243 result = 0.;
244 for (j=0; j<dim; j++)
245 result -= (x[j]-mean[j]) * (covar_inv[idx(coord,j)]+covar_inv[idx(j,coord)]);
246 result *= .5*(dim+DISTR.nu)/(DISTR.nu+xx);
247
248 return result;
249
250 #undef idx
251 } /* end of _unur_pdlogpdf_multistudent() */
252
253 /*---------------------------------------------------------------------------*/
254
255 int
_unur_set_params_multistudent(UNUR_DISTR * distr,double df)256 _unur_set_params_multistudent( UNUR_DISTR *distr, double df )
257 {
258 /* check parameter m (degrees of freedom) */
259 if ( df <= 0. ) {
260 _unur_error(distr_name,UNUR_ERR_DISTR_DOMAIN,"nu <= 0");
261 return UNUR_ERR_DISTR_DOMAIN;
262 }
263
264 /* copy parameters for standard form */
265 DISTR.nu = df;
266
267 /* store number of parameters */
268 DISTR.n_params = 1;
269
270 return UNUR_SUCCESS;
271 } /* end of _unur_set_params_multistudent() */
272
273 /*---------------------------------------------------------------------------*/
274
275 int
_unur_upd_mode_multistudent(UNUR_DISTR * distr)276 _unur_upd_mode_multistudent( UNUR_DISTR *distr )
277 {
278 /* TODO: checking if mode is inside domain */
279
280 if (DISTR.mode == NULL)
281 DISTR.mode = _unur_xmalloc( distr->dim * sizeof(double) );
282 memcpy( DISTR.mode, DISTR.mean, distr->dim * sizeof(double) );
283
284 return UNUR_SUCCESS;
285 } /* end of _unur_upd_mode_multistudent() */
286
287 /*---------------------------------------------------------------------------*/
288
289 int
_unur_upd_volume_multistudent(UNUR_DISTR * distr)290 _unur_upd_volume_multistudent( UNUR_DISTR *distr )
291 {
292 /* TODO: checking for modified domain */
293
294 double det_covar;
295
296 /* log of normalization constant */
297 det_covar = (DISTR.covar == NULL)
298 ? 1. : _unur_matrix_determinant(distr->dim, DISTR.covar);
299 LOGNORMCONSTANT = _unur_SF_ln_gamma((distr->dim+DISTR.nu)/2.)
300 - _unur_SF_ln_gamma(DISTR.nu/2.)
301 - ( distr->dim * log(DISTR.nu*M_PI) + log(det_covar) ) / 2.;
302
303 return UNUR_SUCCESS;
304 } /* end of _unur_upd_volume_multistudent() */
305
306 /*****************************************************************************/
307 /** **/
308 /** Make distribution object **/
309 /** **/
310 /*****************************************************************************/
311
312 /*---------------------------------------------------------------------------*/
313
314 struct unur_distr *
unur_distr_multistudent(int dim,double df,const double * mean,const double * covar)315 unur_distr_multistudent( int dim, double df, const double *mean, const double *covar )
316 /*
317 @code{UNUR_DISTR *unur_distr_multistudent(int dim, const double nu, const double *mean, const double *covar)}
318 creates a distribution object for the multivariate Student t-distribution with
319 @var{dim} components and @var{nu} degrees of freedom.
320 @var{mean} is an array of size @var{dim}.
321 A NULL pointer for @var{mean} is interpreted as the zero
322 vector (0,@dots{},0).
323 @var{covar} is an array of size @var{dim}x@var{dim} and holds the
324 covariance matrix, where the rows of the matrix are stored
325 consecutively in this array. The NULL pointer can be used
326 instead the identity matrix.
327 If @var{covar} is not a valid covariance matrix (i.e., not positive
328 definite) then no distribution object is created and NULL is returned.
329
330 For standard form of the distribution use the null vector for @var{mean} and
331 the identity matrix for @var{covar}.
332 */
333 {
334 struct unur_distr *distr;
335 double det_covar; /* determinant of covariance matrix */
336
337 /* get new (empty) distribution object */
338 distr = unur_distr_cvec_new(dim);
339
340 /* check new parameter for generator */
341 if (distr == NULL) {
342 /* error: dim < 1 */
343 return NULL;
344 }
345
346 /* set distribution id */
347 distr->id = UNUR_DISTR_MSTUDENT;
348
349 /* name of distribution */
350 distr->name = distr_name;
351
352 /* how to get special generators */
353 DISTR.init = NULL;
354
355 /* copy (and check) parameters */
356 if ( (_unur_set_params_multistudent(distr, df)!=UNUR_SUCCESS) ||
357 (unur_distr_cvec_set_mean(distr,mean)!=UNUR_SUCCESS) ||
358 (unur_distr_cvec_set_covar(distr,covar)!=UNUR_SUCCESS) ) {
359 unur_distr_free( distr );
360 return NULL;
361 }
362
363 /* functions */
364 DISTR.pdf = _unur_pdf_multistudent; /* pointer to PDF */
365 DISTR.logpdf = _unur_logpdf_multistudent; /* pointer to logPDF */
366 DISTR.dpdf = _unur_distr_cvec_eval_dpdf_from_dlogpdf; /* pointer to gradient of PDF */
367 DISTR.dlogpdf = _unur_dlogpdf_multistudent; /* pointer to gradient of logPDF */
368 DISTR.pdpdf = _unur_distr_cvec_eval_pdpdf_from_pdlogpdf; /* pointer to part. deriv. of PDF */
369 DISTR.pdlogpdf = _unur_pdlogpdf_multistudent; /* pointer to partial derivative of logPDF */
370
371 #ifdef USE_DEPRECATED_CODE
372 /* set standardized marginal distributions */
373 {
374 struct unur_distr *stdmarginal = unur_distr_student(&df,1);
375 unur_distr_cvec_set_stdmarginals(distr,stdmarginal);
376 unur_distr_free(stdmarginal);
377 }
378 #endif
379
380 /* domain */
381
382 /* log of normalization constant */
383 /* constant: Gamma((dim+1)/2) / ( pi^((dim+1)/2) * sqrt(det(Sigma)) ) */
384
385 /* constant: Gamma((dim+df)/2) */
386 /* / ( Gamma(df/2) (df*pi)^(dim/2) * sqrt(det(Sigma)) ) */
387 det_covar = (DISTR.covar == NULL) ? 1. : _unur_matrix_determinant(dim, DISTR.covar);
388 LOGNORMCONSTANT = _unur_SF_ln_gamma((distr->dim+df)/2.) - _unur_SF_ln_gamma(df/2.)
389 - ( distr->dim * log(df*M_PI) + log(det_covar) ) / 2.;
390
391 /* mode */
392 DISTR.mode = _unur_xmalloc( distr->dim * sizeof(double) );
393 memcpy( DISTR.mode, DISTR.mean, distr->dim * sizeof(double) );
394
395 /* volume below p.d.f. */
396 DISTR.volume = 1.;
397
398 /* function for updating derived parameters */
399 DISTR.upd_mode = _unur_upd_mode_multistudent; /* funct for computing mode */
400 DISTR.upd_volume = _unur_upd_volume_multistudent; /* funct for computing volume */
401
402 /* indicate which parameters are set (additional to mean and covariance) */
403 distr->set |= ( UNUR_DISTR_SET_STDDOMAIN |
404 UNUR_DISTR_SET_PDFVOLUME |
405 UNUR_DISTR_SET_MODE );
406
407 /* return pointer to object */
408 return distr;
409
410 } /* end of unur_distr_multistudent() */
411
412 /*---------------------------------------------------------------------------*/
413 #undef DISTR
414 /*---------------------------------------------------------------------------*/
415