diff --git a/miscprob.cc b/miscprob.cc
index 37396b1ff92e85273b6529e1ef55c9b1320afedb..830424d5dce2e18a3d165272c5af47eccbafe2b2 100644
--- a/miscprob.cc
+++ b/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;
diff --git a/miscprob.h b/miscprob.h
index 879186a3d57a67db7dfa6bca165d2be9b2f0edff..44099c9ef57837224f1a8000e58afb20a41cee53 100644
--- a/miscprob.h
+++ b/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);