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