fork(9) download
copy 
  1. // File: alias_method.cpp
  2. // Author: Jaewon Jung (jj@notsharingmy.info)
  3. //
  4. // This is a C++ port of Keith's Java implementation here:
  5. //  http://w...content-available-to-author-only...z.com/interesting/code/?dir=alias-method
  6. //
  7. // An implementation of the alias method implemented using Vose's algorithm.
  8. // The alias method allows for efficient sampling of random values from a
  9. // discrete probability distribution (i.e. rolling a loaded die) in O(1) time
  10. // each after O(n) preprocessing time.
  11. //
  12. // For a complete writeup on the alias method, including the intuition and
  13. // important proofs, please see the article "Darts, Dice, and Coins: Smpling
  14. // from a Discrete Distribution" at
  15. //
  16. //                 http://w...content-available-to-author-only...z.com/darts-dice-coins/
  17. //
  18. #include <random>
  19. #include <vector>
  20. #include <algorithm>
  21. #include <cstdio>
  22. #include <ctime>
  23.  
  24. // Some utility class and function for wrapping the C rand() library function 
  25. struct rand_engine
  26. {
  27. 	int operator() ()
  28. 	{
  29. 		return rand();
  30. 	}
  31. };
  32.  
  33. template <class T>
  34. double gen_uniform_normalized_double(T& random)
  35. {
  36. 	std::uniform_real_distribution<double> dist(0.0, 1.0);
  37. 	return dist(random);
  38. }
  39.  
  40. template <>
  41. double gen_uniform_normalized_double<rand_engine>(rand_engine& random)
  42. {
  43. 	return static_cast<double>(rand()) / static_cast<double>(RAND_MAX);
  44. }
  45.  
  46. template <class T>
  47. class alias_method
  48. {
  49. public:
  50. 	/**
  51. 	 * Constructs a new alias_method to sample from a discrete distribution and
  52. 	 * hand back outcomes based on the probability distribution.
  53. 	 * <p>
  54. 	 * Given as input a list of probabilities corresponding to outcomes 0, 1,
  55. 	 * ..., n - 1, along with the random number generator that should be used
  56. 	 * as the underlying generator, this constructor creates the probability 
  57. 	 * and alias tables needed to efficiently sample from this distribution.
  58. 	 *
  59. 	 * @param probabilities The list of probabilities.
  60. 	 * @param random The random number generator
  61. 	 */
  62. 	alias_method(const std::vector<double>& probability, T& random)
  63. 		: probability_(probability), random_(random)
  64. 	{
  65.     // Allocate space for the alias table.
  66. 		alias_.resize(probability_.size());
  67.  
  68. 		// Compute the average probability and cache it for later use.
  69. 		const double average = 1.0 / probability_.size();
  70.  
  71. 		// two stacks to act as worklists as we populate the tables
  72. 		std::vector<int> small, large;
  73.  
  74. 		// Populate the stacks with the input probabilities.
  75. 		for(size_t i=0; i<probability_.size(); ++i)
  76. 		{
  77. 			// If the probability is below the average probability, then we add
  78. 			// it to the small list; otherwise we add it to the large list.
  79. 			if (probability_[i] >= average)
  80. 				large.push_back(i);
  81. 			else
  82. 				small.push_back(i);
  83. 		}
  84.  
  85. 		// As a note: in the mathematical specification of the algorithm, we
  86. 		// will always exhaust the small list before the big list.  However,
  87. 		// due to floating point inaccuracies, this is not necessarily true.
  88. 		// Consequently, this inner loop (which tries to pair small and large
  89. 		// elements) will have to check that both lists aren't empty.
  90. 		while(!small.empty() && !large.empty()) 
  91. 		{
  92. 			// Get the index of the small and the large probabilities.
  93. 			int less = small.back(); small.pop_back();
  94. 			int more = large.back(); large.pop_back();
  95.  
  96. 			alias_[less] = more;
  97.  
  98. 			// Decrease the probability of the larger one by the appropriate
  99. 			// amount.
  100. 			probability_[more] = (probability_[more] + probability_[less]) - average;
  101.  
  102. 			// If the new probability is less than the average, add it into the
  103. 			// small list; otherwise add it to the large list.
  104. 			if (probability_[more] >= average)
  105. 				large.push_back(more);
  106. 			else
  107. 				small.push_back(more);
  108. 		}
  109.  
  110. 		// At this point, everything is in one list, which means that the
  111. 		// remaining probabilities should all be 1/n.  Based on this, set them
  112. 		// appropriately.  Due to numerical issues, we can't be sure which
  113. 		// stack will hold the entries, so we empty both.
  114. 		while(!small.empty())
  115. 		{
  116. 			probability_[small.back()] = average;
  117. 			small.pop_back();
  118. 		}
  119. 		while(!large.empty())
  120. 		{
  121. 			probability_[large.back()] = average;
  122. 			large.pop_back();
  123. 		}
  124.  
  125. 		// These probabilities have not yet been scaled up to be such that
  126. 		// 1/n is given weight 1.0.  We do this here instead.
  127. 		int n = static_cast<int>(probability_.size());
  128. 		std::transform(probability_.cbegin(), probability_.cend(), probability_.begin(), 
  129. 			[n](double p){ return p * n; });
  130. 	}
  131.  
  132. 	/**
  133. 	 * Samples a value from the underlying distribution.
  134. 	 *
  135. 	 * @return A random value sampled from the underlying distribution.
  136. 	 */
  137. 	int next()
  138. 	{
  139. 		// Generate a fair die roll to determine which column to inspect.
  140. 		int column = random_() % probability_.size();
  141.  
  142. 		// Generate a biased coin toss to determine which option to pick.
  143. 		bool coinToss = gen_uniform_normalized_double(random_) < probability_[column];
  144.  
  145. 		// Based on the outcome, return either the column or its alias.
  146. 		return coinToss? column : alias_[column];
  147. 	}
  148.  
  149. private:
  150.   // The probability and alias tables.
  151. 	std::vector<int> alias_;
  152. 	std::vector<double> probability_;
  153.   // The random number generator used to sample from the distribution.
  154. 	T& random_;
  155. };
  156.  
  157. template <class T>
  158. alias_method<T> make_alias_method(const std::vector<double>& probabilities, T& random)
  159. {
  160. 	return alias_method<T>(probabilities, random);
  161. }
  162.  
  163. int main(int argc, char* argv[])
  164. {
  165. 	std::vector<double> probabilities;
  166. 	probabilities.push_back(1.0/20.0);
  167. 	probabilities.push_back(7.0/20.0);
  168. 	probabilities.push_back(3.0/20.0);
  169. 	probabilities.push_back(5.0/20.0);
  170. 	probabilities.push_back(4.0/20.0);
  171.  
  172. 	//rand_engine rnd;
  173. 	//srand(clock());
  174. 	std::mt19937 rnd;
  175. 	rnd.seed(clock());
  176. 	auto am = make_alias_method(probabilities, rnd);
  177.  
  178. 	const int tries = 1000000;
  179. 	int count[5] = { 0, 0, 0, 0, 0 };
  180. 	for(int i=0; i<tries; ++i)
  181. 	{
  182. 		++count[am.next()];
  183. 	}
  184.  
  185. 	for(int i=0; i<5; ++i)
  186. 		printf("%f\n", 20.0*static_cast<double>(count[i])/static_cast<double>(tries));
  187.  
  188. 	return 0;
  189. }
  190.  
Success #stdin #stdout 0.08s 2964KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
0.993480
6.989840
3.000700
5.017620
3.998360
https://ideone.com/gzFzI
language:
C++14 (gcc 8.3)
created:
12 years ago
visibility:
secret

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!