From fb5779b24b2a3302fe1b6adb946d300d1a407674 Mon Sep 17 00:00:00 2001
From: Saad Jbabdi <saad@fmrib.ox.ac.uk>
Date: Tue, 3 Jul 2007 08:53:53 +0000
Subject: [PATCH] added sampling from a distribution given a histogram

---
 miscprob.cc | 20 ++++++++++++++++++++
 miscprob.h  |  2 ++
 2 files changed, 22 insertions(+)

diff --git a/miscprob.cc b/miscprob.cc
index 37396b1..830424d 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 879186a..44099c9 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);
-- 
GitLab