1 /*****************************************************************************
2 * *
3 * UNURAN -- Universal Non-Uniform Random number generator *
4 * *
5 *****************************************************************************
6 * *
7 * Distribution objects of multivariate distributions having a *
8 * correlation matrix with constant off-diagonal elements *
9 * *
10 *****************************************************************************
11 * *
12 * Use correlation matrix of the form: *
13 * *
14 * +- -+ *
15 * | 1 r r ... r | *
16 * | r 1 r ... r | *
17 * | r r 1 ... r | *
18 * | ... | *
19 * | r r r ... 1 | *
20 * +- -+ *
21 * *
22 * in this case the inverse matrix is given as *
23 * *
24 * +- -+ *
25 * | a b b ... b b | *
26 * | b a b ... b b | *
27 * | b b a ... b b | *
28 * | ... | *
29 * | b b b ... a b | *
30 * | b b b ... b a | *
31 * +- -+ *
32 * *
33 * with *
34 * *
35 * a = (1+(d-2)*r) / (1+(d-2)*r-(d-1)*r^2) *
36 * b = - r / (1+(d-2)*r-(d-1)*r^2) *
37 * *
38 * and the determinant of the covariance matrix is *
39 * *
40 * det = (1-r)^(d-1) * (1+(d-1)*r) *
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 <unuran.h>
67 #include <config.h>
68 #include <unuran_config.h>
69 #include <distr/distr_struct.h>
70 #include <utils/matrix_source.h>
71 #include <utils/unur_fp_source.h>
72 #include <utils/unur_math_source.h>
73 #include <specfunct/unur_specfunct_source.h>
74 #include <float.h>
75 #include "testdistributions.h"
76
77 /*---------------------------------------------------------------------------*/
78
79 #define LOGNORMCONSTANT (distr->data.cvec.norm_constant)
80 #define NORMCONSTANT (distr->data.cvec.norm_constant)
81
82 /*---------------------------------------------------------------------------*/
83
84 /* Set correlation matrix with constant off-diagonal elements and its inverse */
85 static int _unur_vc_set_corrmatrix_constantrho( UNUR_DISTR *distr, int dim, double rho );
86
87 /*---------------------------------------------------------------------------*/
88
89 int
_unur_vc_set_corrmatrix_constantrho(UNUR_DISTR * distr,int dim,double rho)90 _unur_vc_set_corrmatrix_constantrho( UNUR_DISTR *distr, int dim, double rho )
91 /*----------------------------------------------------------------------*/
92 /* Create and set correlation matrix with constant off-diagonal */
93 /* elements and its inverse. */
94 /* */
95 /* parameters: */
96 /* distr ... pointer to distribution object */
97 /* dim ... dimension of random vector */
98 /* rho ... correlation between consecutive elements in process */
99 /* */
100 /* return: */
101 /* UNUR_SUCCESS ... success */
102 /* UNUR_FAILURE ... failure */
103 /*----------------------------------------------------------------------*/
104 {
105 #define idx(a,b) ((a)*dim+(b))
106
107 double *covar, *covar_inv;
108 double a, b, denominator;
109 int i,j;
110 int error = UNUR_FAILURE;
111
112 /* checking parameters */
113 if (_unur_iszero(rho)) /* nothing to do */
114 return UNUR_SUCCESS;
115
116 if ( rho<0. || rho>=1. || dim<1)
117 /* invalid arguments */
118 return UNUR_FAILURE;
119
120 /* entries of matrix */
121 denominator = 1.+(dim-2)*rho-(dim-1)*rho*rho;
122 a = (1.+(dim-2)*rho)/denominator;
123 b = -rho/denominator;
124
125 /* allocate memory */
126 covar = malloc( dim * dim * sizeof(double) );
127 covar_inv = malloc( dim * dim * sizeof(double) );
128 if (covar == NULL || covar_inv == NULL) {
129 if (covar) free (covar);
130 if (covar_inv) free (covar_inv);
131 return UNUR_FAILURE;
132 }
133
134 /* create correlation matrix */
135 for (i=0; i<dim; i++)
136 for (j=0; j<dim; j++)
137 covar[idx(i,j)] = (i==j)? 1.: rho;
138
139 /* create inverse of correlation matrix */
140 for (i=0; i<dim; i++)
141 for (j=0; j<dim; j++)
142 covar_inv[idx(i,j)] = (i==j) ? a : b;
143
144
145 /* set covariance matrix and its inverse */
146 if ( (unur_distr_cvec_set_covar( distr, covar ) == UNUR_SUCCESS) &&
147 (unur_distr_cvec_set_covar_inv( distr, covar_inv ) == UNUR_SUCCESS) )
148 error = UNUR_SUCCESS;
149
150 /* free working space */
151 free (covar); free (covar_inv);
152
153 /* return errorcode */
154 return error;
155
156 #undef idx
157 } /* end of _unur_vc_set_corrmatrix_constantrho() */
158
159 /*---------------------------------------------------------------------------*/
160
161 struct unur_distr *
unur_distr_multinormal_constantrho(int dim,const double * mean,double rho)162 unur_distr_multinormal_constantrho(int dim, const double *mean, double rho)
163 /*---------------------------------------------------------------------------*/
164 /* Multinormal distribution (corr-matrix : constant off-diagonal elements) */
165 /*---------------------------------------------------------------------------*/
166 {
167 struct unur_distr *distr;
168 double det_covar;
169
170 /* checking parameters */
171 if ( rho<0. || rho>=1. || dim<1)
172 return NULL;
173
174 /* get distribution object for multinormal distribution */
175 distr = unur_distr_multinormal_w_marginals( dim, mean, NULL );
176
177 /* set the name of distribution */
178 unur_distr_set_name(distr, "multinormal_constantrho");
179
180 /* set the correlation matrix and its inverse */
181 if (_unur_vc_set_corrmatrix_constantrho(distr, dim, rho) != UNUR_SUCCESS) {
182 /* error */
183 unur_distr_free(distr); distr = NULL; return NULL;
184 }
185
186 /* compute normalization constant */
187 det_covar = _unur_matrix_determinant(distr->dim, distr->data.cvec.covar);
188 LOGNORMCONSTANT = - ( distr->dim * log(2 * M_PI) + log(det_covar) ) / 2.;
189
190 /* return pointer to object */
191 return distr;
192
193 } /* end of unur_distr_multinormal_constantrho() */
194
195 /*---------------------------------------------------------------------------*/
196
197 struct unur_distr *
unur_distr_multicauchy_constantrho(int dim,const double * mean,double rho)198 unur_distr_multicauchy_constantrho(int dim, const double *mean, double rho)
199 /*---------------------------------------------------------------------------*/
200 /* Multinormal distribution (corr-matrix : constant off-diagonal elements) */
201 /*---------------------------------------------------------------------------*/
202 {
203 struct unur_distr *distr;
204 double det_covar;
205
206 /* checking parameters */
207 if ( rho<0. || rho>=1. || dim<1)
208 return NULL;
209
210 /* get distribution object for multicauchy distribution */
211 distr = unur_distr_multicauchy_w_marginals( dim, mean, NULL );
212
213 /* set the name of distribution */
214 unur_distr_set_name(distr, "multicauchy_constantrho");
215
216 /* set the correlation matrix and its inverse */
217 if (_unur_vc_set_corrmatrix_constantrho(distr, dim, rho) != UNUR_SUCCESS) {
218 /* error */
219 unur_distr_free(distr); distr = NULL; return NULL;
220 }
221
222 /* compute normalization constant */
223 det_covar = _unur_matrix_determinant(distr->dim, distr->data.cvec.covar);
224 LOGNORMCONSTANT = _unur_SF_ln_gamma((distr->dim+1)/2.)
225 - ( (distr->dim+1) * log(M_PI) + log(det_covar) ) / 2.;
226
227 /* return pointer to object */
228 return distr;
229
230 } /* end of unur_distr_multicauchy_constantrho() */
231
232 /*---------------------------------------------------------------------------*/
233
234 struct unur_distr *
unur_distr_multistudent_constantrho(int dim,double df,const double * mean,double rho)235 unur_distr_multistudent_constantrho(int dim, double df, const double *mean, double rho)
236 /*---------------------------------------------------------------------------*/
237 /* Multistudent distribution (corr-matrix : constant off-diagonal elements) */
238 /*---------------------------------------------------------------------------*/
239 {
240 struct unur_distr *distr;
241 double det_covar;
242
243 /* checking parameters */
244 if ( rho<0. || rho>=1. || dim<1)
245 return NULL;
246
247 /* get distribution object for multistudent distribution */
248 distr = unur_distr_multistudent_w_marginals( dim, df, mean, NULL );
249
250 /* set the name of distribution */
251 unur_distr_set_name(distr, "multistudent_constantrho");
252
253 /* set the correlation matrix and its inverse */
254 if (_unur_vc_set_corrmatrix_constantrho(distr, dim, rho) != UNUR_SUCCESS) {
255 /* error */
256 unur_distr_free(distr); distr = NULL; return NULL;
257 }
258
259 /* compute normalization constant */
260 det_covar = _unur_matrix_determinant(distr->dim, distr->data.cvec.covar);
261 LOGNORMCONSTANT = _unur_SF_ln_gamma((distr->dim+df)/2.) - _unur_SF_ln_gamma(df/2.)
262 - ( distr->dim * log(df*M_PI) + log(det_covar) ) / 2.;
263
264 /* return pointer to object */
265 return distr;
266
267 } /* end of unur_distr_multistudent_constantrho() */
268
269 /*---------------------------------------------------------------------------*/
270