added wishart distribution for infinite gmm

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);

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
     {