|
Packit |
67cb25 |
/* randist/wishart.c
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 2017 Timothée Flutre
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is free software; you can redistribute it and/or modify
|
|
Packit |
67cb25 |
* it under the terms of the GNU General Public License as published by
|
|
Packit |
67cb25 |
* the Free Software Foundation; either version 3 of the License, or (at
|
|
Packit |
67cb25 |
* your option) any later version.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is distributed in the hope that it will be useful, but
|
|
Packit |
67cb25 |
* WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
Packit |
67cb25 |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
Packit |
67cb25 |
* General Public License for more details.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* You should have received a copy of the GNU General Public License
|
|
Packit |
67cb25 |
* along with this program; if not, write to the Free Software
|
|
Packit |
67cb25 |
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <config.h>
|
|
Packit |
67cb25 |
#include <math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_rng.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_randist.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_vector.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_matrix.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_blas.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_linalg.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_sf_gamma.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Generate a random matrix from a Wishart distribution using the Bartlett
|
|
Packit |
67cb25 |
* decomposition, following Smith and Hocking, Journal of the Royal Statistical
|
|
Packit |
67cb25 |
* Society. Series C (Applied Statistics), Vol. 21, No. 3 (1972), pp. 341-345.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* df degrees of freedom
|
|
Packit |
67cb25 |
* L matrix resulting from the Cholesky decomposition of
|
|
Packit |
67cb25 |
* the scale matrix V = L L^T (dimension d x d)
|
|
Packit |
67cb25 |
* result output matrix (dimension d x d)
|
|
Packit |
67cb25 |
* work matrix used for intermediate computations (dimension d x d)
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_ran_wishart (const gsl_rng * r,
|
|
Packit |
67cb25 |
const double df,
|
|
Packit |
67cb25 |
const gsl_matrix * L,
|
|
Packit |
67cb25 |
gsl_matrix * result,
|
|
Packit |
67cb25 |
gsl_matrix * work)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (L->size1 != L->size2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR("L should be a square matrix", GSL_ENOTSQR);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (result->size1 != result->size2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR("result should be a square matrix", GSL_ENOTSQR);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (work->size1 != work->size2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR("work should be a square matrix", GSL_ENOTSQR);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (result->size1 != L->size1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR("incompatible dimensions of result matrix", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (work->size1 != L->size1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR("incompatible dimensions of work matrix", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (df <= L->size1 - 1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR("incompatible degrees of freedom", GSL_EDOM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* result: X = L A A^T L^T */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
size_t d = L->size1, i, j;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* insure the upper part of A is zero before filling its lower part */
|
|
Packit |
67cb25 |
gsl_matrix_set_zero(work);
|
|
Packit |
67cb25 |
for (i = 0; i < d; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix_set(work, i, i, sqrt(gsl_ran_chisq(r, df - i)));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (j = 0; j < i; ++j)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix_set(work, i, j, gsl_ran_ugaussian(r));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute L * A */
|
|
Packit |
67cb25 |
gsl_blas_dtrmm(CblasLeft, CblasLower, CblasNoTrans, CblasNonUnit, 1.0,
|
|
Packit |
67cb25 |
L, work);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute (L * A) * (L * A)^T */
|
|
Packit |
67cb25 |
gsl_blas_dsyrk(CblasUpper, CblasNoTrans, 1.0, work, 0.0, result);
|
|
Packit |
67cb25 |
for (i = 0; i < d; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
for (j = 0; j < i; ++j)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix_set(result, i, j, gsl_matrix_get(result, j, i));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute the log of the probability density function at a given quantile
|
|
Packit |
67cb25 |
* matrix for a Wishart distribution using the Cholesky decomposition of the
|
|
Packit |
67cb25 |
* scale matrix.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* X quantile matrix at which to evaluate the PDF (dimension d x d)
|
|
Packit |
67cb25 |
* L_X matrix resulting from the Cholesky decomposition of
|
|
Packit |
67cb25 |
* of the quantile matrix at which to evaluate the PDF
|
|
Packit |
67cb25 |
* X = L_X L_X^T (dimension d x d)
|
|
Packit |
67cb25 |
* df degrees of freedom
|
|
Packit |
67cb25 |
* L matrix resulting from the Cholesky decomposition of
|
|
Packit |
67cb25 |
* the scale matrix V = L L^T (dimension d x d)
|
|
Packit |
67cb25 |
* result output of the density (dimension 1)
|
|
Packit |
67cb25 |
* work matrix used for intermediate computations (dimension d x d)
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_ran_wishart_log_pdf (const gsl_matrix * X,
|
|
Packit |
67cb25 |
const gsl_matrix * L_X,
|
|
Packit |
67cb25 |
const double df,
|
|
Packit |
67cb25 |
const gsl_matrix * L,
|
|
Packit |
67cb25 |
double * result,
|
|
Packit |
67cb25 |
gsl_matrix * work)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (L->size1 != L->size2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR("L should be a square matrix", GSL_ENOTSQR);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (X->size1 != X->size2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR("X should be a square matrix", GSL_ENOTSQR);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (L_X->size1 != L_X->size2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR("L_X should be a square matrix", GSL_ENOTSQR);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (X->size1 != L->size1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR("incompatible dimensions of X matrix", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (L_X->size1 != L->size1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR("incompatible dimensions of L_X matrix", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (df <= L->size1 - 1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR("incompatible degrees of freedom", GSL_EDOM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t d = L->size1, i;
|
|
Packit |
67cb25 |
int status;
|
|
Packit |
67cb25 |
double log_mv_Ga, log_det_V, log_det_X, tr_Vinv_X;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute the log of the multivariate Gamma */
|
|
Packit |
67cb25 |
log_mv_Ga = d * (d-1) * 0.25 * log(M_PI);
|
|
Packit |
67cb25 |
for (i = 0; i < d; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
log_mv_Ga += gsl_sf_lngamma((df - i + 1) * 0.5);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute the log of the determinant of the scale matrix */
|
|
Packit |
67cb25 |
log_det_V = log(gsl_matrix_get(L, 0, 0));
|
|
Packit |
67cb25 |
for (i = 1; i < d; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
log_det_V += log(gsl_matrix_get(L, i, i));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
log_det_V = 2 * log_det_V;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute the log of the determinant of the quantile matrix */
|
|
Packit |
67cb25 |
log_det_X = log(gsl_matrix_get(L_X, 0, 0));
|
|
Packit |
67cb25 |
for (i = 1; i < d; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
log_det_X += log(gsl_matrix_get(L_X, i, i));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
log_det_X = 2 * log_det_X;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute the trace of V^(-1) X */
|
|
Packit |
67cb25 |
status = gsl_linalg_cholesky_solve_mat(L, X, work);
|
|
Packit |
67cb25 |
if (status)
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
tr_Vinv_X = gsl_matrix_get(work, 0, 0);
|
|
Packit |
67cb25 |
for (i = 1; i < d; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
tr_Vinv_X += gsl_matrix_get(work, i, i);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* add all to get the log of the PDF */
|
|
Packit |
67cb25 |
*result = - (0.5 * df * d) * log(2.0)
|
|
Packit |
67cb25 |
- (0.5 * df) * log_det_V
|
|
Packit |
67cb25 |
- log_mv_Ga
|
|
Packit |
67cb25 |
+ 0.5 * (df - d - 1) * log_det_X
|
|
Packit |
67cb25 |
- 0.5 * tr_Vinv_X;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_ran_wishart_pdf (const gsl_matrix * X,
|
|
Packit |
67cb25 |
const gsl_matrix * L_X,
|
|
Packit |
67cb25 |
const double df,
|
|
Packit |
67cb25 |
const gsl_matrix * L,
|
|
Packit |
67cb25 |
double * result,
|
|
Packit |
67cb25 |
gsl_matrix * work)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double logpdf;
|
|
Packit |
67cb25 |
int status = gsl_ran_wishart_log_pdf(X, L_X, df, L, &logpdf, work);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (status == GSL_SUCCESS)
|
|
Packit |
67cb25 |
*result = exp(logpdf);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
}
|