Source code by Ralph Silva
from
here
/***************************************************************************************
* Multivariate Normal density function and random number generator
* Multivariate Student t density function and random number generator
* Wishart random number generator
* Using GSL -> www.gnu.org/software/gsl
*
* Copyright (C) 2006 Ralph dos Santos Silva
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
* AUTHOR
* Ralph dos Santos Silva, [EMAIL PROTECTED]
* March, 2006
***************************************************************************************/
#include <math.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
/*****************************************************************************************************************/
/*****************************************************************************************************************/
int rmvnorm(const gsl_rng *r, const int n, const gsl_vector *mean, const gsl_matrix *var, gsl_vector *result){
/* multivariate normal distribution random number generator */
/*
* n dimension of the random vetor
* mean vector of means of size n
* var variance matrix of dimension n x n
* result output variable with a sigle random vector normal distribution generation
*/
int k;
gsl_matrix *work = gsl_matrix_alloc(n,n);
gsl_matrix_memcpy(work,var);
gsl_linalg_cholesky_decomp(work);
for(k=0; k<n; k++)
gsl_vector_set( result, k, gsl_ran_ugaussian(r) );
gsl_blas_dtrmv(CblasLower, CblasNoTrans, CblasNonUnit, work, result);
gsl_vector_add(result,mean);
gsl_matrix_free(work);
return 0;
}
/*****************************************************************************************************************/
/*****************************************************************************************************************/
double dmvnorm(const int n, const gsl_vector *x, const gsl_vector *mean, const gsl_matrix *var){
/* multivariate normal density function */
/*
* n dimension of the random vetor
* mean vector of means of size n
* var variance matrix of dimension n x n
*/
int s;
double ax,ay;
gsl_vector *ym, *xm;
gsl_matrix *work = gsl_matrix_alloc(n,n),
*winv = gsl_matrix_alloc(n,n);
gsl_permutation *p = gsl_permutation_alloc(n);
gsl_matrix_memcpy( work, var );
gsl_linalg_LU_decomp( work, p, &s );
gsl_linalg_LU_invert( work, p, winv );
ax = gsl_linalg_LU_det( work, s );
gsl_matrix_free( work );
gsl_permutation_free( p );
xm = gsl_vector_alloc(n);
gsl_vector_memcpy( xm, x);
gsl_vector_sub( xm, mean );
ym = gsl_vector_alloc(n);
gsl_blas_dsymv(CblasUpper,1.0,winv,xm,0.0,ym);
gsl_matrix_free( winv );
gsl_blas_ddot( xm, ym, &ay);
gsl_vector_free(xm);
gsl_vector_free(ym);
ay = exp(-0.5*ay)/sqrt( pow((2*M_PI),n)*ax );
return ay;
}
/*****************************************************************************************************************/
/*****************************************************************************************************************/
int rmvt(const gsl_rng *r, const int n, const gsl_vector *location, const gsl_matrix *scale, const int dof, gsl_vector *result){
/* multivariate Student t distribution random number generator */
/*
* n dimension of the random vetor
* location vector of locations of size n
* scale scale matrix of dimension n x n
* dof degrees of freedom
* result output variable with a single random vector normal distribution generation
*/
int k;
gsl_matrix *work = gsl_matrix_alloc(n,n);
double ax = 0.5*dof;
ax = gsl_ran_gamma(r,ax,(1/ax)); /* gamma distribution */
gsl_matrix_memcpy(work,scale);
gsl_matrix_scale(work,(1/ax)); /* scaling the matrix */
gsl_linalg_cholesky_decomp(work);
for(k=0; k<n; k++)
gsl_vector_set( result, k, gsl_ran_ugaussian(r) );
gsl_blas_dtrmv(CblasLower, CblasNoTrans, CblasNonUnit, work, result);
gsl_vector_add(result, location);
gsl_matrix_free(work);
return 0;
}
/*****************************************************************************************************************/
/*****************************************************************************************************************/
double dmvt(const int n, const gsl_vector *x, const gsl_vector *location, const gsl_matrix *scale, const int dof){
/* multivariate Student t density function */
/*
* n dimension of the random vetor
* location vector of locations of size n
* scale scale matrix of dimension n x n
* dof degrees of freedom
*/
int s;
double ax,ay,az=0.5*(dof + n);
gsl_vector *ym, *xm;
gsl_matrix *work = gsl_matrix_alloc(n,n),
*winv = gsl_matrix_alloc(n,n);
gsl_permutation *p = gsl_permutation_alloc(n);
gsl_matrix_memcpy( work, scale );
gsl_linalg_LU_decomp( work, p, &s );
gsl_linalg_LU_invert( work, p, winv );
ax = gsl_linalg_LU_det( work, s );
gsl_matrix_free( work );
gsl_permutation_free( p );
xm = gsl_vector_alloc(n);
gsl_vector_memcpy( xm, x);
gsl_vector_sub( xm, location );
ym = gsl_vector_alloc(n);
gsl_blas_dsymv(CblasUpper,1.0,winv,xm,0.0,ym);
gsl_matrix_free( winv );
gsl_blas_ddot( xm, ym, &ay);
gsl_vector_free(xm);
gsl_vector_free(ym);
ay = pow((1+ay/dof),-az)*gsl_sf_gamma(az)/(gsl_sf_gamma(0.5*dof)*sqrt( pow((dof*M_PI),n)*ax ));
return ay;
}
/*****************************************************************************************************************/
/*****************************************************************************************************************/
int rwishart(const gsl_rng *r, const int n, const int dof, const gsl_matrix *scale, gsl_matrix *result){
/* Wishart distribution random number generator */
/*
* n gives the dimension of the random matrix
* dof degrees of freedom
* scale scale matrix of dimension n x n
* result output variable with a single random matrix Wishart distribution generation
*/
int k,l;
gsl_matrix *work = gsl_matrix_calloc(n,n);
for(k=0; k<n; k++){
gsl_matrix_set( work, k, k, sqrt( gsl_ran_chisq( r, (dof-k) ) ) );
for(l=0; l<k; l++){
gsl_matrix_set( work, k, l, gsl_ran_ugaussian(r) );
}
}
gsl_matrix_memcpy(result,scale);
gsl_linalg_cholesky_decomp(result);
gsl_blas_dtrmm(CblasLeft,CblasLower,CblasNoTrans,CblasNonUnit,1.0,result,work);
gsl_blas_dsyrk(CblasUpper,CblasNoTrans,1.0,work,0.0,result);
return 0;
}
/***************************************************************************************
* Multivariate Normal density function and random number generator
* Multivariate Student t density function and random number generator
* Wishart random number generator
* Using GSL -> www.gnu.org/software/gsl
*
* Copyright (C) 2006 Ralph dos Santos Silva
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
* AUTHOR
* Ralph dos Santos Silva, [EMAIL PROTECTED]
* March, 2006
***************************************************************************************/
int rmvnorm(const gsl_rng *r, const int n, const gsl_vector *mean, const gsl_matrix *var, gsl_vector *result);
double dmvnorm(const int n, const gsl_vector *x, const gsl_vector *mean, const gsl_matrix *var);
int rmvt(const gsl_rng *r, const int n, const gsl_vector *location, const gsl_matrix *scale, const int dof, gsl_vector *result);
double dmvt(const int n, const gsl_vector *x, const gsl_vector *locationn, const gsl_matrix *scale, const int dof);
int rwishart(const gsl_rng *r, const int n, const int dof, const gsl_matrix *scale, gsl_matrix *result);
/***************************************************************************************
* A testing program for multivariate distributions
* Using GSL -> www.gnu.org/software/gsl
* Copyright (C) 2006 Ralph dos Santos Silva
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
* AUTHOR
* Ralph dos Santos Silva, [EMAIL PROTECTED]
* March, 2006
***************************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <stddef.h>
#include <time.h>
#include <math.h>
/* GSL - GNU Scientific Library */
/* #define GSL_CHECK_RANGE_OFF */
#include <gsl/gsl_math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_cdf.h>
/* ----------------------------- */
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
/* ----------------------------- */
#include "rmv.h"
int main(){
FILE *f;
int k;
unsigned long int initime;
double result;
gsl_vector *x = gsl_vector_calloc(3),
*mean = gsl_vector_calloc(3);
gsl_matrix *m = gsl_matrix_alloc(3,3),
*rm = gsl_matrix_alloc(3,3);
gsl_rng *r;
time(&initime);
r = gsl_rng_alloc(gsl_rng_mt19937);
gsl_rng_set(r, initime);
printf("The SEED is %li.\n\n",initime);
gsl_matrix_set(m,0,0, 1.0);
gsl_matrix_set(m,0,1, 0.2);
gsl_matrix_set(m,0,2,-0.9);
gsl_matrix_set(m,1,0, 0.2);
gsl_matrix_set(m,1,1, 1.0);
gsl_matrix_set(m,1,2, 0.1);
gsl_matrix_set(m,2,0,-0.9);
gsl_matrix_set(m,2,1, 0.1);
gsl_matrix_set(m,2,2, 1.0);
gsl_vector_set(x,0,0.1);
gsl_vector_set(x,1,0.1);
gsl_vector_set(x,2,0.1);
result = dmvnorm(3,x,mean,m);
fprintf(stdout,"norm = %g\n", result);
gsl_vector_set(x,0,0.1);
gsl_vector_set(x,1,0.1);
gsl_vector_set(x,2,0.1);
result = dmvt(3,x,mean,m,5);
fprintf(stdout,"t = %g\n", result);
f = fopen("test-n.txt","w");
for(k=0; k<10; k++){
rmvnorm(r,3,mean,m,x);
fprintf(f, "%g %g %g\n",gsl_vector_get(x,0),gsl_vector_get(x,1),gsl_vector_get(x,2));
}
fclose(f);
f = fopen("test-t.txt","w");
for(k=0; k<10; k++){
rmvt(r,3,mean,m,5,x);
fprintf(f, "%g %g %g\n",gsl_vector_get(x,0),gsl_vector_get(x,1),gsl_vector_get(x,2));
}
fclose(f);
f = fopen("test-w.txt","w");
for(k=0; k<1000; k++){
rwishart(r,3,5,m,rm);
fprintf(f, "%g %g %g\n",gsl_matrix_get(rm,0,0),gsl_matrix_get(rm,0,1),gsl_matrix_get(rm,0,2));
fprintf(f, "%g %g %g\n",gsl_matrix_get(rm,1,0),gsl_matrix_get(rm,1,1),gsl_matrix_get(rm,1,2));
fprintf(f, "%g %g %g\n",gsl_matrix_get(rm,2,0),gsl_matrix_get(rm,2,1),gsl_matrix_get(rm,2,2));
}
fclose(f);
return 0;
}
Makefile
# Copyright (C) 2006 Ralph dos Santos Silva
#LIBDIR = /home/rsilva1/local/libsun
LIBDIR = /usr/lib
INCLUDEDIR = /usr/include
LIBGSL = -lgsl
LIBCBLAS = -lcblas
LIBATLAS = -latlas
LIBMATH = -lm
CFLAGS = -O2 -static -g -Wall -W -DHAVE_INLINE -I$(INCLUDEDIR)
CC = gcc
.c.o:
$(CC) $(CFLAGS) -c $<
all: test.o rmv.o
$(CC) -o test test.o rmv.o -L$(LIBDIR) $(LIBGSL) $(LIBCBLAS) $(LIBATLAS) $(LIBMATH)
clean:
rm -f *.o test
