add a rmultinom() function to Eidos

Author: bhaller

EidosScribe/EidosHelpFunctions.rtf

@@ -2187,6 +2187,32 @@ Returns a vector of
 \f2 \
 \pard\pardeftab397\li720\fi-446\ri720\sb180\sa60\partightenfactor0
 
+\f1\fs18 \cf2 (integer)rmultinom(integer$\'a0n, integer$\'a0size, numeric\'a0prob)\
+\pard\pardeftab397\li547\ri720\sb60\sa60\partightenfactor0
+
+\f3\fs20 \cf2 Returns a matrix of 
+\f1\fs18 n
+\f3\fs20  
+\f0\b random vectors from the specified multinomial distribution
+\f3\b0 .  For each random vector drawn from the multinomial distribution, 
+\f1\fs18 size
+\f3\fs20  objects will be put into 
+\f7\i k
+\f3\i0  boxes, where 
+\f1\fs18 prob
+\f3\fs20  is a vector of probabilities of length 
+\f7\i k
+\f3\i0  specifying the probability for each box.  If 
+\f1\fs18 prob
+\f3\fs20  is not normalized to sum to 
+\f1\fs18 1.0
+\f3\fs20 , its entries are treated as weights and normalized appropriately.  The draws are returned as a matrix with 
+\f7\i k
+\f3\i0  rows (one row per box) and 
+\f1\fs18 n
+\f3\fs20  columns (one column per drawn random vector).\
+\pard\pardeftab397\li720\fi-446\ri720\sb180\sa60\partightenfactor0
+
 \f1\fs18 \cf2 \expnd0\expndtw0\kerning0
 (float)rmvnorm(integer$\'a0n, numeric\'a0mu, numeric\'a0sigma)\
 \pard\pardeftab720\li547\ri720\sb60\sa60\partightenfactor0

QtSLiM/help/EidosHelpFunctions.html

@@ -194,6 +194,8 @@
 <p class="p5">Returns a vector of <span class="s2">n</span> <b>random draws from a Laplace distribution</b> with shape parameter <span class="s2">b</span>.<span class="Apple-converted-space">  </span>The <span class="s2">b</span> parameter may either be a singleton, specifying a single value to be used for all of the draws, or it may be a vector of length <span class="s2">n</span>, specifying a value for each draw.</p>
 <p class="p2">(float)rlnorm(integer$<span class="s1"> </span>n, [numeric meanlog = 0], [numeric sdlog = 1])</p>
 <p class="p3">Returns a vector of <span class="s2">n</span> <b>random draws from a lognormal distribution</b> with mean <span class="s2">meanlog</span> and standard deviation <span class="s2">sdlog</span>, specified on the log scale.<span class="Apple-converted-space">  </span>The <span class="s2">meanlog</span> and <span class="s2">sdlog</span> parameters may either be singletons, specifying a single value to be used for all of the draws, or they may be vectors of length <span class="s2">n</span>, specifying a value for each draw.</p>
+<p class="p4">(integer)rmultinom(integer$ n, integer$ size, numeric prob)</p>
+<p class="p5">Returns a matrix of <span class="s2">n</span> <b>random vectors from the specified multinomial distribution</b>.<span class="Apple-converted-space">  </span>For each random vector drawn from the multinomial distribution, <span class="s2">size</span> objects will be put into <i>k</i> boxes, where <span class="s2">prob</span> is a vector of probabilities of length <i>k</i> specifying the probability for each box.<span class="Apple-converted-space">  </span>If <span class="s2">prob</span> is not normalized to sum to <span class="s2">1.0</span>, its entries are treated as weights and normalized appropriately.<span class="Apple-converted-space">  </span>The draws are returned as a matrix with <i>k</i> rows (one row per box) and <span class="s2">n</span> columns (one column per drawn random vector).</p>
 <p class="p4"><span class="s5">(float)rmvnorm(integer$ n, numeric mu, numeric sigma)</span></p>
 <p class="p5"><span class="s5">Returns a matrix of </span><span class="s9">n</span><span class="s5"> <b>random draws from a <i>k</i>-dimensional multivariate normal distribution</b> with a length <i>k</i> mean vector </span><span class="s9">mu</span><span class="s5"> and a <i>k</i> × <i>k</i> variance-covariance matrix </span><span class="s9">sigma</span><span class="s5">.<span class="Apple-converted-space">  </span>The </span><span class="s9">mu</span><span class="s5"> and </span><span class="s9">sigma</span><span class="s5"> parameters are used for all </span><span class="s9">n</span><span class="s5"> draws.<span class="Apple-converted-space">  </span>The draws are returned as a matrix with </span><span class="s9">n</span><span class="s5"> rows (one row per draw) and <i>k</i> columns (one column per dimension).<span class="Apple-converted-space">  </span>The number of dimensions <i>k</i> must be at least two; for <i>k</i>=1, use </span><span class="s9">rnorm()</span><span class="s5">.</span></p>
 <p class="p5"><span class="s5">Cholesky decomposition of the variance-covariance matrix </span><span class="s9">sigma</span><span class="s5"> is involved as an internal step, and this requires that </span><span class="s9">sigma</span><span class="s5"> be positive-definite; if it is not, an error will result.<span class="Apple-converted-space">  </span>When more than one draw is needed, it is much more efficient to call </span><span class="s9">rmvnorm()</span><span class="s5"> once to generate all of the draws, since the Cholesky decomposition of </span><span class="s9">sigma</span><span class="s5"> can then be done just once.</span></p>

VERSIONS

@@ -6,6 +6,7 @@ Note that not every commit will be logged here; that is what the Github commit h
 
 
 development head (in the master branch):
+	add a rmultinom() function to draw from a multinomial distribution, matching rmultinom() in R
 
 
 version 5.2 (Eidos version 4.2):

eidos/eidos_functions.cpp

@@ -185,6 +185,7 @@ const std::vector<EidosFunctionSignature_CSP> &EidosInterpreter::BuiltInFunction
 		signatures->emplace_back((EidosFunctionSignature *)(new EidosFunctionSignature("rgeom",				Eidos_ExecuteFunction_rgeom,		kEidosValueMaskInt))->AddInt_S(gEidosStr_n)->AddFloat("p"));
 		signatures->emplace_back((EidosFunctionSignature *)(new EidosFunctionSignature("rlaplace",			Eidos_ExecuteFunction_rlaplace,		kEidosValueMaskFloat))->AddInt_S(gEidosStr_n)->AddNumeric_O("b", gStaticEidosValue_Integer1));
 		signatures->emplace_back((EidosFunctionSignature *)(new EidosFunctionSignature("rlnorm",			Eidos_ExecuteFunction_rlnorm,		kEidosValueMaskFloat))->AddInt_S(gEidosStr_n)->AddNumeric_O("meanlog", gStaticEidosValue_Integer0)->AddNumeric_O("sdlog", gStaticEidosValue_Integer1));
+		signatures->emplace_back((EidosFunctionSignature *)(new EidosFunctionSignature("rmultinom",			Eidos_ExecuteFunction_rmultinom,	kEidosValueMaskInt))->AddInt_S(gEidosStr_n)->AddInt_S(gEidosStr_size)->AddNumeric("prob"));
 		signatures->emplace_back((EidosFunctionSignature *)(new EidosFunctionSignature("rmvnorm",			Eidos_ExecuteFunction_rmvnorm,		kEidosValueMaskFloat))->AddInt_S(gEidosStr_n)->AddNumeric("mu")->AddNumeric("sigma"));
 		signatures->emplace_back((EidosFunctionSignature *)(new EidosFunctionSignature("rnbinom",			Eidos_ExecuteFunction_rnbinom,		kEidosValueMaskInt))->AddInt_S(gEidosStr_n)->AddNumeric("size")->AddFloat("prob"));
 		signatures->emplace_back((EidosFunctionSignature *)(new EidosFunctionSignature("rnorm",				Eidos_ExecuteFunction_rnorm,		kEidosValueMaskFloat))->AddInt_S(gEidosStr_n)->AddNumeric_O("mean", gStaticEidosValue_Integer0)->AddNumeric_O("sd", gStaticEidosValue_Integer1));

eidos/eidos_functions.h

@@ -122,6 +122,7 @@ EidosValue_SP Eidos_ExecuteFunction_rgamma(const std::vector<EidosValue_SP> &p_a
 EidosValue_SP Eidos_ExecuteFunction_rgeom(const std::vector<EidosValue_SP> &p_arguments, EidosInterpreter &p_interpreter);
 EidosValue_SP Eidos_ExecuteFunction_rlaplace(const std::vector<EidosValue_SP> &p_arguments, EidosInterpreter &p_interpreter);
 EidosValue_SP Eidos_ExecuteFunction_rlnorm(const std::vector<EidosValue_SP> &p_arguments, EidosInterpreter &p_interpreter);
+EidosValue_SP Eidos_ExecuteFunction_rmultinom(const std::vector<EidosValue_SP> &p_arguments, EidosInterpreter &p_interpreter);
 EidosValue_SP Eidos_ExecuteFunction_rmvnorm(const std::vector<EidosValue_SP> &p_arguments, EidosInterpreter &p_interpreter);
 EidosValue_SP Eidos_ExecuteFunction_rnbinom(const std::vector<EidosValue_SP> &p_arguments, EidosInterpreter &p_interpreter);
 EidosValue_SP Eidos_ExecuteFunction_rnorm(const std::vector<EidosValue_SP> &p_arguments, EidosInterpreter &p_interpreter);

eidos/eidos_functions_distributions.cpp

@@ -1656,6 +1656,108 @@ EidosValue_SP Eidos_ExecuteFunction_rlnorm(const std::vector<EidosValue_SP> &p_a
 	return result_SP;
 }
 
+// (integer)rmultinom(integer$ n, integer$ size, numeric prob)
+EidosValue_SP Eidos_ExecuteFunction_rmultinom(const std::vector<EidosValue_SP> &p_arguments, __attribute__((unused)) EidosInterpreter &p_interpreter)
+{
+	EidosValue_SP result_SP(nullptr);
+	
+	EidosValue *arg_n = p_arguments[0].get();
+	EidosValue *arg_size = p_arguments[1].get();
+	EidosValue *arg_prob = p_arguments[2].get();
+	
+	// get and bounds-check the number of draws we will do
+	int64_t num_draws = arg_n->IntAtIndex_NOCAST(0, nullptr);
+	
+	if (num_draws < 1)
+		EIDOS_TERMINATION << "ERROR (Eidos_ExecuteFunction_rmultinom): function rmultinom() requires n to be greater than or equal to 1 (" << num_draws << " supplied)." << EidosTerminate(nullptr);
+	
+	// get and check the size of each multinomial draw
+	// the upper limit is because the GSL requires size to be cast down to unsigned int
+	int64_t size = arg_size->IntAtIndex_NOCAST(0, nullptr);
+	
+	if (size < 0)
+		EIDOS_TERMINATION << "ERROR (Eidos_ExecuteFunction_rmultinom): function rmultinom() requires size to be greater than or equal to 0." << EidosTerminate(nullptr);
+	if (size > 1000000000)
+		EIDOS_TERMINATION << "ERROR (Eidos_ExecuteFunction_rmultinom): function rmultinom() requires size to be less than or equal to 1000000000." << EidosTerminate(nullptr);
+	
+	// get and check the number of "boxes", i.e., the number of probabilities supplied
+	// the upper limit is arbitrary, but it seems wise to have a limit
+	int k = arg_prob->Count();
+	
+	if (k <= 0)
+		EIDOS_TERMINATION << "ERROR (Eidos_ExecuteFunction_rmultinom): function rmultinom() requires prob to be of length 1 or greater." << EidosTerminate(nullptr);
+	if (k > 1000000)
+		EIDOS_TERMINATION << "ERROR (Eidos_ExecuteFunction_rmultinom): function rmultinom() requires prob to be of length 1000000 or less." << EidosTerminate(nullptr);
+	
+	// fetch the probability vector into a local static buffer, checking as we go; note the gsl normalizes prob_vec for us
+	double prob_sum = 0.0;
+	static std::vector<double> prob_vec;
+	prob_vec.reserve(k);
+	prob_vec.resize(0);
+	
+	for (int prob_index = 0; prob_index < k; ++prob_index)
+	{
+		double prob = arg_prob->FloatAtIndex_CAST(prob_index, nullptr);
+		
+		if (prob < 0.0)
+			EIDOS_TERMINATION << "ERROR (Eidos_ExecuteFunction_rmultinom): function rmultinom() requires all probabilities in prob to be >= 0.0." << EidosTerminate(nullptr);
+		if (!std::isfinite(prob))
+			EIDOS_TERMINATION << "ERROR (Eidos_ExecuteFunction_rmultinom): function rmultinom() requires all probabilities in prob to be finite (not INF or NAN)." << EidosTerminate(nullptr);
+		
+		prob_sum += prob;
+		prob_vec.push_back(prob);
+	}
+	
+	if (prob_sum == 0.0)
+		EIDOS_TERMINATION << "ERROR (Eidos_ExecuteFunction_rmultinom): function rmultinom() requires the sum of prob to be greater than zero; there must be at least one non-zero probability." << EidosTerminate(nullptr);
+	if (!std::isfinite(prob_sum))
+		EIDOS_TERMINATION << "ERROR (Eidos_ExecuteFunction_rmultinom): function rmultinom() could not compute the multinomial draws because the sum of prob overflowed to infinity." << EidosTerminate(nullptr);
+	
+	if (k == 1)
+	{
+		// If k == 1, every object goes into the same single box, so the result is trivial
+		EidosValue_Int *int_result = (new (gEidosValuePool->AllocateChunk()) EidosValue_Int())->resize_no_initialize((int)num_draws);
+		result_SP = EidosValue_SP(int_result);
+		
+		for (int draw_index = 0; draw_index < num_draws; ++draw_index)
+			int_result->set_int_no_check(size, draw_index);
+	}
+	else if (size == 0)
+	{
+		// if size == 0, zero objects are drawn in each multinomial trial, so the result is trivial
+		EidosValue_Int *int_result = (new (gEidosValuePool->AllocateChunk()) EidosValue_Int())->reserve((int)num_draws * k);
+		result_SP = EidosValue_SP(int_result);
+		
+		for (int draw_index = 0; draw_index < num_draws; ++draw_index)
+			for (int k_index = 0; k_index < k; ++k_index)
+				int_result->push_int_no_check(0);
+	}
+	else
+	{
+		// the main case, where we need to call gsl_ran_multinomial() to do the work
+		// we put the drawn values into a local static buffer, and then copy them into the Eidos result
+		static std::vector<unsigned int> draw_vec;
+		draw_vec.resize(k);
+		
+		gsl_rng *rng_gsl = EIDOS_GSL_RNG(omp_get_thread_num());
+		EidosValue_Int *int_result = (new (gEidosValuePool->AllocateChunk()) EidosValue_Int())->reserve((int)num_draws * k);
+		result_SP = EidosValue_SP(int_result);
+		
+		for (int draw_index = 0; draw_index < num_draws; ++draw_index)
+		{
+			gsl_ran_multinomial(rng_gsl, k, (unsigned int)size, prob_vec.data(), draw_vec.data());
+			
+			for (int k_index = 0; k_index < k; ++k_index)
+				int_result->push_int_no_check(draw_vec[k_index]);
+		}
+	}
+	
+	int64_t dim_buffer[2] = {k, num_draws};
+	result_SP->SetDimensions(2, dim_buffer);
+	
+	return result_SP;
+}
+
 // (float)rmvnorm(integer$ n, numeric mu, numeric sigma)
 EidosValue_SP Eidos_ExecuteFunction_rmvnorm(const std::vector<EidosValue_SP> &p_arguments, __attribute__((unused)) EidosInterpreter &p_interpreter)
 {

eidos/eidos_test_builtins.h

@@ -1160,6 +1160,16 @@ if (abs(m - 5) > 0.02) stop('Mismatch in expectation vs. realization of rlnorm()
 
 // ***********************************************************************************************
 
+setSeed(asInteger(clock() * 100000));
+x = rmultinom(10, 1000000000, c(0.1, 0.5, 0.1, 5.0, 0.25));
+if (!identical(colSums(x), rep(1000000000, 10))) stop('ERROR (rmultinom): (internal error) colSums incorrect');
+r = rowSums(x);
+norm = c(0.1, 0.5, 0.1, 5.0, 0.25) / sum(c(0.1, 0.5, 0.1, 5.0, 0.25));
+expected = 1000000000 * norm * 10;
+if (sum(abs(r - expected)) > 300000) stop('Mismatch in expectation vs. realization of rmultinom() - could be random chance (but very unlikely), rerun test');
+
+// ***********************************************************************************************
+
 setSeed(asInteger(clock() * 100000));
 x = rmvnorm(100000, c(-1, 5), matrix(c(1, 0.2, 0.2, 2), nrow=2));
 m1 = mean(x[,0]);

eidos/eidos_test_functions_statistics.cpp

@@ -961,6 +961,17 @@ void _RunFunctionDistributionTests(void)
 	EidosAssertScriptSuccess("rlnorm(1, NAN, 100);", gStaticEidosValue_FloatNAN);
 	EidosAssertScriptSuccess("rlnorm(1, 1, NAN);", gStaticEidosValue_FloatNAN);
 
+	// rmultinom()
+	EidosAssertScriptRaise("rmultinom(0, 10, c(0.1, 1.0, 10.0));", 0, "requires n to be");
+	EidosAssertScriptRaise("rmultinom(5, -1, c(0.1, 1.0, 10.0));", 0, "requires size to be");
+	EidosAssertScriptRaise("rmultinom(5, 10, float(0));", 0, "requires prob to be");
+	EidosAssertScriptRaise("rmultinom(5, 10, c(-0.1, 1.0, 10.0));", 0, "requires all probabilities in prob to be");
+	EidosAssertScriptRaise("rmultinom(5, 10, c(INF, 1.0, 10.0));", 0, "requires all probabilities in prob to be");
+	EidosAssertScriptRaise("rmultinom(5, 10, c(NAN, 1.0, 10.0));", 0, "requires all probabilities in prob to be");
+	EidosAssertScriptRaise("rmultinom(5, 10, c(0.0, 0.0, 0.0));", 0, "requires the sum of prob");
+	EidosAssertScriptSuccess_L("x = rmultinom(5, 10, c(0.1)); identical(x, matrix(rep(10, 5*1), ncol=5));", true);
+	EidosAssertScriptSuccess_L("x = rmultinom(5, 0, c(0.1, 1.0, 10.0)); identical(x, matrix(rep(0, 5*3), ncol=5));", true);
+	
 	// rmvnorm()
 	EidosAssertScriptRaise("rmvnorm(0, c(0,2), matrix(c(10,3), nrow=2));", 0, "requires n to be");
 	EidosAssertScriptRaise("rmvnorm(5, matrix(c(0,0)), matrix(c(10,3,3,2), nrow=2));", 0, "plain vector of length k");