Cybernetics

/***************************************************************************************
 *  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;
}

# 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