From c3ba8588852d6dd5677104a0d371693ba0174507 Mon Sep 17 00:00:00 2001
From: Saad Jbabdi <saad@fmrib.ox.ac.uk>
Date: Mon, 14 May 2007 17:52:24 +0000
Subject: [PATCH] added wishart distribution for infinite gmm
---
miscprob.cc | 71 +++++++++++++++++++++++++++++++++++++++++++++++++++++
miscprob.h | 7 ++++++
2 files changed, 78 insertions(+)
diff --git a/miscprob.cc b/miscprob.cc
index 00b08a3..276de51 100644
--- a/miscprob.cc
+++ b/miscprob.cc
@@ -66,7 +66,33 @@ ReturnMatrix normpdf(const RowVector& vals, const float mu, const float var)
res.Release();
return res;
}
+float normpdf(const ColumnVector& vals, const ColumnVector& mu, const SymmetricMatrix& sigma)
+{
+ float res;
+ LogAndSign ld=(2*M_PI*sigma).LogDeterminant();
+
+ res = std::exp(-0.5*(
+ ((vals-mu).t()*sigma.i()*(vals-mu)).AsScalar()
+ +ld.LogValue()
+ ));
+ return res;
+}
+ReturnMatrix normpdf(const Matrix& vals, const ColumnVector& mu, const SymmetricMatrix& sigma)
+{
+ RowVector res(vals);
+ LogAndSign ld = (2*M_PI*sigma).LogDeterminant();
+ Matrix isigma = sigma.i();
+
+ for (int mc=1; mc<=res.Ncols(); mc++){
+ res(mc) = std::exp( -0.5*(
+ ((vals.Column(mc)-mu).t()*isigma*(vals.Column(mc)-mu)).AsScalar()
+ +ld.LogValue()
+ ));
+ }
+ res.Release();
+ return res;
+}
ReturnMatrix normcdf(const RowVector& vals, const float mu, const float var)
{
RowVector res(vals);
@@ -145,6 +171,51 @@ ReturnMatrix mvnrnd(const RowVector& mu, const SymmetricMatrix& covar, int nsamp
return mvn.next(nsamp);
}
+// Saad: Wishart and inverseWishart Random Generator
+ReturnMatrix wishrnd(const SymmetricMatrix& sigma,const int dof){
+ // compute cholesky factor for sigma
+ LowerTriangularMatrix L = Cholesky(sigma);
+
+ // for small degrees of freedom, use the definition
+ int n = sigma.Nrows();
+ Matrix X;
+ if(dof <= 81+n ){
+ X.ReSize(dof,n);
+ X = normrnd(dof,n) * L.t();
+ }
+ // otherwise, use Smith & Hocking procedure
+ else{
+ X.ReSize(n,n);
+ Matrix A(n,n);
+ for(int i=1;i<=n;i++){
+ Gamma G((dof-i+1)/2);
+ G.Set(rand()/float(RAND_MAX));
+ for(int j=1;j<=n;j++){
+ if (i>j) { A(i,j) = 0; }
+ else if(i<j) { A(i,j) = normrnd(1,1).AsScalar(); }
+ else { A(i,j) = std::sqrt(2*G.Next()); }
+ }
+ }
+ X = A * L.t();
+ }
+
+ SymmetricMatrix res(n);
+ res << X.t() * X;
+
+ res.Release();
+ return res;
+}
+ReturnMatrix iwishrnd(const SymmetricMatrix& sigma,const int dof){
+ // assumes inv-Wishart(sigma.i(),dof)
+
+ SymmetricMatrix res;
+ res = wishrnd(sigma,dof);
+ res = res.i();
+
+ res.Release();
+ return res;
+}
+
ReturnMatrix perms(const int n){
if(n<=1){
Matrix P(1,1);
diff --git a/miscprob.h b/miscprob.h
index f0fb80b..127de4f 100644
--- a/miscprob.h
+++ b/miscprob.h
@@ -17,6 +17,7 @@
#include "stdlib.h"
using namespace NEWMAT;
+using namespace NEWRAN;
namespace MISCMATHS {
@@ -35,6 +36,10 @@ namespace MISCMATHS {
ReturnMatrix normpdf(const RowVector& vals, const RowVector& mus,
const RowVector& vars);
+ float normpdf(const ColumnVector& val,const ColumnVector& mu,const SymmetricMatrix& sigma);
+
+ ReturnMatrix normpdf(const Matrix& val,const ColumnVector& mu,const SymmetricMatrix& sigma);
+
ReturnMatrix normcdf(const RowVector& vals, const float mu = 0, const float var = 1);
ReturnMatrix gammapdf(const RowVector& vals, const float mu = 0, const float var = 1);
@@ -44,6 +49,8 @@ namespace MISCMATHS {
// returns n! * n matrix of all possible permutations
ReturnMatrix perms(const int n);
+ ReturnMatrix wishrnd(const SymmetricMatrix&,const int);
+ ReturnMatrix iwishrnd(const SymmetricMatrix&,const int);
class Mvnormrandm
{
--
GitLab