*** empty log message ***

histogram.cc

@@ -54,31 +54,45 @@ namespace MISCMATHS {
 
       // smooth in i direction
       newhist=0;
+      ColumnVector kernel(3); 
+      // corresponds to Gaussian with sigma=0.8 voxels
+      //       kernel(1)=0.5;
+      //       kernel(2)=0.2283;      
+      //       kernel(3)=0.0219;
+      // corresponds to Gaussian with sigma=0.6 voxels
+      //       kernel(1)=0.6638;
+      //       kernel(2)=0.1655;      
+      //       kernel(3)=0.0026;
+
+      //gauss(0.5,5,1)
+      kernel(1)=0.7866;
+      kernel(2)=0.1065;      
+      kernel(3)=0.0003;
 
       for(int i=1; i<=bins; i++)
 	  {
 	    float val=0.5*histogram(i);
-	    float norm=0.5;
+	    float norm=kernel(1);
 
 	    if(i>1)
 	      {
-		val+=0.2283*(histogram(i-1));
-		norm+=0.2283;
+		val+=kernel(2)*(histogram(i-1));
+		norm+=kernel(2);
 	      }
 	    if(i>2)
 	      {
-		val+=0.0219*(histogram(i-2));
-		norm+=0.0219;		
+		val+=kernel(3)*(histogram(i-2));
+		norm+=kernel(3);		
 	      }
 	    if(i<bins)
 	      {
-		val+=0.2283*(histogram(i+1));
-		norm+=0.2283;
+		val+=kernel(2)*(histogram(i+1));
+		norm+=kernel(2);
 	      }
 	    if(i<bins-1)
 	      {
-		val+=0.0219*(histogram(i+2));
-		norm+=0.0219;		
+		val+=kernel(3)*(histogram(i+2));
+		norm+=kernel(3);		
 	      }
 	    val/=norm;
 

miscprob.cc

@@ -121,6 +121,7 @@ ReturnMatrix normpdf(const RowVector& vals, const float mu, const float var)
   return res;
 }
 
+
 ReturnMatrix normcdf(const RowVector& vals, const float mu, const float var)
 {
   RowVector res(vals);
@@ -176,7 +177,11 @@ ReturnMatrix gammapdf(const RowVector& vals, const float mu, const float var)
   {
     return std::exp(-0.5*(std::pow(val-mu,2)/var))*std::pow(2*M_PI*var,-0.5);
   }
-  
+  float lognormpdf(const float val, const float mu, const float var)
+  {
+    return -0.5*(std::pow(val-mu,2)/var+std::log(2*M_PI*var));
+  }
+    
 ReturnMatrix normpdf(const RowVector& vals, const RowVector& mu, const RowVector& var)
 {
   Matrix res(mu.Ncols(),vals.Ncols());
@@ -205,13 +210,27 @@ ReturnMatrix mvnrnd(const RowVector& mu, const SymmetricMatrix& covar, int nsamp
 
  float mvnpdf(const RowVector& vals, const RowVector& mu, const SymmetricMatrix& covar)
  {
-   return std::exp(-0.5*((vals-mu)*covar.i()*(vals-mu).t()).AsScalar())/(std::pow(covar.Determinant(),0.5)*std::pow(2*M_PI,vals.Ncols()/2.0));
+   if(vals.Ncols()==2)
+     return bvnpdf(vals,mu,covar);
+   else
+     return std::exp(-0.5*((vals-mu)*covar.i()*(vals-mu).t()).AsScalar())/(std::pow(covar.Determinant(),0.5)*std::pow(2*M_PI,vals.Ncols()/2.0));
+ }
+
+ float bvnpdf(const RowVector& vals, const RowVector& mu, const SymmetricMatrix& covar)
+ {
+   // bivariate normal pdf
+   double det=covar(1,1)*covar(2,2)-Sqr(covar(1,2));
+   float m1=vals(1)-mu(1);
+   float m2=vals(2)-mu(2);
+   float ss=(Sqr(m1)*covar(2,2)-2*m1*m2*covar(1,2)+Sqr(m2)*covar(1,1))/det;
+
+   return std::exp(-0.5*ss)/(std::pow(det,0.5)*std::pow(2*M_PI,vals.Ncols()/2.0));
  }
 
 // 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.
+//   // 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;}

miscprob.h

@@ -37,7 +37,10 @@ namespace MISCMATHS {
   
   float mvnpdf(const RowVector& vals, const RowVector& mu, const SymmetricMatrix& covar);
 
+  float bvnpdf(const RowVector& vals, const RowVector& mu, const SymmetricMatrix& covar);
+
   float normpdf(const float val, const float mu = 0, const float var = 1);
+  float lognormpdf(const float val, const float mu = 0, const float var = 1);
 
   ReturnMatrix normpdf(const RowVector& vals, const float mu = 0, const float var = 1);