added sampling from a distribution given a histogram

miscprob.cc

@@ -36,6 +36,26 @@ 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;

miscprob.h

@@ -23,6 +23,8 @@ namespace MISCMATHS {
  
   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);