*** empty log message ***

Makefile

@@ -4,8 +4,8 @@ include ${FSLCONFDIR}/default.mk
 
 PROJNAME = miscmaths
 
-USRINCFLAGS = -I${INC_NEWMAT} -I${INC_BOOST} -I${INC_NEWRAN} -I${INC_PROB} -I${INC_ZLIB}
-USRLDFLAGS = -L${LIB_NEWMAT} -L${LIB_NEWRAN} -L${LIB_PROB} -L${LIB_ZLIB}
+USRINCFLAGS = -I${INC_NEWMAT} -I${INC_BOOST} -I${INC_PROB} -I${INC_ZLIB}
+USRLDFLAGS = -L${LIB_NEWMAT} -L${LIB_PROB} -L${LIB_ZLIB}
 
 OBJS = miscmaths.o optimise.o miscprob.o kernel.o histogram.o base2z.o t2z.o f2z.o volume.o volumeseries.o minimize.o cspline.o sparse_matrix.o sparsefn.o rungekutta.o nonlin.o bfmatrix.o
 #OBJS = miscmaths.o optimise.o miscprob.o kernel.o histogram.o base2z.o t2z.o f2z.o volume.o volumeseries.o minimize.o cspline.o
@@ -22,6 +22,3 @@ quick:${OBJS} quick.o
 
 libmiscmaths.a: ${OBJS}
 	${AR} -r libmiscmaths.a ${OBJS}
-
-
-

miscprob.cc

@@ -12,11 +12,65 @@
 #include "stdlib.h"
 #include "newmatio.h"
 #include <iostream>
+// #include "gam.h"
 
 using namespace NEWMAT;
 
 namespace MISCMATHS {
 
+//   ReturnMatrix betarnd(const int dim1, const int dim2, const float a, const float b)
+//   {
+//     // Devroye, L. (1986) Non-Uniform Random Variate Generation, Springer-Verlag.
+
+//     int tdim = dim2;
+
+//     if(tdim<0){tdim=dim1;}
+
+    
+//     Matrix g1=gammarnd(dim1, tdim, a, 1);
+//     Matrix g2=gammarnd(dim1, tdim, b, 1);
+
+//     Matrix res(dim1,tdim);
+//     for (int mc=1; mc<=res.Ncols(); mc++) {
+//       for (int mr=1; mr<=res.Nrows(); mr++) {
+// 	res(mr,mc)=g1(mr,mc)/(g1(mr,mc)+g2(mr,mc));
+//       }
+//     }
+    
+//     res.Release();
+//     return res;
+//   }
+
+ReturnMatrix betapdf(const RowVector& vals, const float a, const float b)
+{
+
+  RowVector res(vals);
+
+  if(a<0 || b<0)
+    {
+      throw Exception("Negative a or b in call to Miscprob::betapdf");
+    }
+
+  for (int mc=1; mc<=res.Ncols(); mc++)
+    {
+      float x=vals(mc);
+      if(x<0)
+	{
+	  res(mc)=0;
+	}
+      else
+	{
+	  float logkerna=(a-1)*std::log(x);
+	  float logkernb=(b-1)*std::log(1-x);
+	  float betaln_ab=lgam(a)+lgam(b)-lgam(a+b);
+	  res(mc)=std::exp(logkerna+logkernb-betaln_ab);	 
+	}
+    }
+  
+  res.Release();
+  return res;
+}
+
 ReturnMatrix unifrnd(const int dim1, const int dim2, const float start, const float end)
 {
   int tdim = dim2;
@@ -36,26 +90,6 @@ ReturnMatrix unifrnd(const int dim1, const int dim2, const float start, const fl
   return res;
 }
 
-  // SJ. Generates a sample from a distribution given the histogram
-int distribrnd(const ColumnVector& histo){
-  int res=1;
-  ColumnVector cumsum(histo.Nrows());
-  float sum=0.0;
-  for(int k=1;k<=histo.Nrows();k++){
-    sum += histo(k);
-    cumsum(k) = sum;
-  }
-  float U=rand()/float(RAND_MAX);
-  U *= sum;
-  for(int k=1;k<=histo.Nrows();k++){
-    if(U<cumsum(k)){
-      res = k;
-      break;
-    }
-  }
-  return res;
-}
-
 ReturnMatrix normrnd(const int dim1, const int dim2, const float mu, const float sigma)
 {
   int tdim = dim2;
@@ -86,33 +120,7 @@ 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);
@@ -191,53 +199,26 @@ 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 gammarnd(const int dim1, const int dim2, 
+// 			const float a, const float b)
+// {
+//   // Marsaglia, G. and Tsang, W.W. (2000) "A Simple Method for Generating Gamma Variables", ACM Trans. Math. Soft. 26(3):363-372.
+
+//   int tdim = dim2;
+//   if(tdim<0){tdim=dim1;}
+//   Matrix res(dim1,tdim);
+
+//   Gam& gam=Gam::getInstance();
+//   gam.setParams(a,b);
+
+//   for (int mc=1; mc<=res.Ncols(); mc++) {
+//     for (int mr=1; mr<=res.Nrows(); mr++) {
+//       res(mr,mc)=gam.rnd();
+//     }
+//   }
+//   res.Release();
+//   return res;
+// }
 
 ReturnMatrix perms(const int n){
   if(n<=1){

miscprob.h

@@ -13,18 +13,21 @@
 #define __miscprob_h
 
 #include "miscmaths.h"
-#include "libprob/libprob.h"
+#include "libprob.h"
 #include "stdlib.h"
 
 using namespace NEWMAT;
-// using namespace NEWRAN;
 
 namespace MISCMATHS {
- 
+
+//   ReturnMatrix betarnd(const int dim1, const int dim2, 
+// 		       const float a, const float b); 
+
+  ReturnMatrix betapdf(const RowVector& vals, 
+		       const float a, const float b); 
+
   ReturnMatrix unifrnd(const int dim1 = 1, const int dim2 = -1, 
 		       const float start = 0, const float end = 1);
-
-  int distribrnd(const ColumnVector& histo);
   
   ReturnMatrix normrnd(const int dim1 = 1, const int dim2 = -1, 
 		       const float mu = 0, const float sigma = 1);
@@ -38,21 +41,18 @@ 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);
 
   ReturnMatrix gammacdf(const RowVector& vals, const float mu = 0, const float var = 1);
 
+//   ReturnMatrix gammarnd(const int dim1, const int dim2, 
+// 			const float a, const float b);
+
   // 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
     {
@@ -76,6 +76,30 @@ namespace MISCMATHS {
 	  ret.Release();
 	  return ret;
 	}
+
+      ReturnMatrix next(const RowVector& pmu, int nsamp = 1)  
+	{
+	  mu=pmu;
+
+	  Matrix ret = ones(nsamp, 1)*mu + normrnd(nsamp,mu.Ncols())*covarw;
+	  ret.Release();
+	  return ret;
+	}
+
+      void setcovar(const SymmetricMatrix& pcovar)
+	{
+	  covar=pcovar;
+
+	  mu.ReSize(covar.Nrows());
+	  mu=0;
+
+	  Matrix eig_vec;
+	  DiagonalMatrix eig_val;
+	  EigenValues(covar,eig_val,eig_vec);
+
+	  covarw = sqrt(eig_val)*eig_vec.t();
+	}
+
     private:      
 
       RowVector mu;