/*
* Copyright (c) 2002-2006 Samit Basu
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*/
#include "Core.hpp"
#include "Exception.hpp"
#include "Malloc.hpp"
#include <math.h>
#include <stdio.h>
#include "RanLib.hpp"
static bool initialized = false;
//!
//@Module SEED Seed the Random Number Generator
//@@Section RANDOM
//@@Usage
//Seeds the random number generator using the given integer seeds.
//Changing the seed allows you to choose which pseudo-random
//sequence is generated. The seed takes two @|uint32| values:
//@[
// seed(s,t)
//@]
//where @|s| and @|t| are the seed values. Note that due to limitations
//in @|ranlib|, the values of @|s,t| must be between @|0 <= s,t <= 2^30|.
//@@Example
//Here's an example of how the seed value can be used to reproduce
//a specific random number sequence.
//@<
//seed(32,41);
//rand(1,5)
//seed(32,41);
//rand(1,5)
//@>
//!
ArrayVector SeedFunction(int nargout, const ArrayVector& arg) {
if (arg.size() != 2)
throw Exception("Seed function requires a two integer arguments");
Array tmp1(arg[0]);
uint32 seedval1;
Array tmp2(arg[1]);
uint32 seedval2;
seedval1 = tmp1.getContentsAsIntegerScalar();
seedval2 = tmp2.getContentsAsIntegerScalar();
setall(seedval1,seedval2);
init_genrand(seedval1);
initialized = true;
return ArrayVector();
}
//!
//@Module RANDBETA Beta Deviate Random Number Generator
//@@Section RANDOM
//@@Usage
//Creates an array of beta random deviates based on the supplied
//two parameters. The general syntax for @|randbeta| is
//@[
// y = randbeta(alpha, beta)
//@]
//where @|alpha| and @|beta| are the two parameters of the
//random deviate. There are three forms for calling @|randbeta|.
//The first uses two vectors @|alpha| and @|beta| of the same
//size, in which case the output @|y| is the same size as both
//inputs, and each deviate uses the corresponding values of @|alpha|
//and @|beta| from the arguments. In the other forms, either
//@|alpha| or @|beta| are scalars.
//@@Function Internals
//The probability density function (PDF) of a beta random variable
//is
//\[
//f(x) = x^(a-1) * (1-x)^(b-1) / B(a,b)
//\]
//for @|x| between 0 and 1. The function @|B(a,b)| is defined so
//that the integral of @|f(x)| is 1.
//@@Example
//Here is a plot of the PDF of a beta random variable with @|a=3|,
//@|b=7|.
//@<
//a = 3; b = 7;
//x = (0:100)/100; t = x.^(a-1).*(1-x).^(b-1);
//t = t/(sum(t)*.01);
//plot(x,t);
//mprint betapdf
//@>
//which is plotted as
//@figure betapdf
//If we generate a few random deviates with these values,
//we see they are distributed around the peak of roughly
//@|0.25|.
//@<
//randbeta(3*ones(1,5),7*ones(1,5))
//@>
//!
ArrayVector RandBetaFunction(int nargout, const ArrayVector& arg) {
if (arg.size() != 2)
throw Exception("randbeta requires two parameter arguments");
Array arg1(arg[0]);
Array arg2(arg[1]);
if (arg1.isEmpty() || arg2.isEmpty())
return singleArrayVector(Array::emptyConstructor());
// Check the logic to see if one or both are scalar values
if (!(arg1.isScalar() || arg2.isScalar() || (arg1.dimensions().equals(arg2.dimensions()))))
throw Exception("randbeta requires either one of the two arguments to be a scalar, or both arguments to be the same size");
int arg1_advance;
int arg2_advance;
arg1_advance = (arg1.isScalar()) ? 0 : 1;
arg2_advance = (arg2.isScalar()) ? 0 : 1;
// Output dimension is the larger of the two
Dimensions outDims;
if (arg1.getLength() > arg2.getLength()) {
outDims = arg1.dimensions();
} else {
outDims = arg2.dimensions();
}
arg1.promoteType(FM_FLOAT);
arg2.promoteType(FM_FLOAT);
float *dp;
dp = (float *) Malloc(sizeof(float)*outDims.getElementCount());
float *p1;
p1 = (float*) arg1.getDataPointer();
float *p2;
p2 = (float*) arg2.getDataPointer();
int i;
for (i=0;i<outDims.getElementCount();i++)
dp[i] = genbet(p1[i*arg1_advance],p2[i*arg2_advance]);
return singleArrayVector(Array(FM_FLOAT,outDims,dp));
}
//!
//@Module RANDI Uniformly Distributed Integer
//@@Section RANDOM
//@@Usage
//Generates an array of uniformly distributed integers between
//the two supplied limits. The general syntax for @|randi| is
//@[
// y = randi(low,high)
//@]
//where @|low| and @|high| are arrays of integers. Scalars
//can be used for one of the arguments. The output @|y| is
//a uniformly distributed pseudo-random number between @|low|
//and @|high| (inclusive).
//@@Example
//Here is an example of a set of random integers between
//zero and 5:
//@<
//randi(zeros(1,6),5*ones(1,6))
//@>
//!
ArrayVector RandIFunction(int nargout, const ArrayVector& arg) {
if (arg.size() != 2)
throw Exception("randi requires two parameter arguments");
Array arg1(arg[0]);
Array arg2(arg[1]);
if (arg1.isEmpty() || arg2.isEmpty())
return singleArrayVector(Array::emptyConstructor());
// Check the logic to see if one or both are scalar values
if (!(arg1.isScalar() || arg2.isScalar() || (arg1.dimensions().equals(arg2.dimensions()))))
throw Exception("randbeta requires either one of the two arguments to be a scalar, or both arguments to be the same size");
int arg1_advance;
int arg2_advance;
arg1_advance = (arg1.isScalar()) ? 0 : 1;
arg2_advance = (arg2.isScalar()) ? 0 : 1;
// Output dimension is the larger of the two
Dimensions outDims;
if (arg1.getLength() > arg2.getLength()) {
outDims = arg1.dimensions();
} else {
outDims = arg2.dimensions();
}
arg1.promoteType(FM_INT32);
arg2.promoteType(FM_INT32);
int32 *dp;
dp = (int32 *) Malloc(sizeof(int32)*outDims.getElementCount());
int32 *p1;
p1 = (int32*) arg1.getDataPointer();
int32 *p2;
p2 = (int32*) arg2.getDataPointer();
int i;
for (i=0;i<outDims.getElementCount();i++)
dp[i] = ignuin(p1[i*arg1_advance],p2[i*arg2_advance]);
return singleArrayVector(Array(FM_INT32,outDims,dp));
}
//!
//@Module RANDCHI Generate Chi-Square Random Variable
//@@Section RANDOM
//@@Usage
//Generates a vector of chi-square random variables with the
//given number of degrees of freedom. The general syntax for
//its use is
//@[
// y = randchi(n)
//@]
//where @|n| is an array containing the degrees of freedom for
//each generated random variable.
//@@Function Internals
//A chi-square random variable is essentially distributed as
//the squared Euclidean norm of a vector of standard Gaussian random
//variables. The number of degrees of freedom is generally the
//number of elements in the vector. In general, the PDF of
//a chi-square random variable is
//\[
// f(x) = \frac{x^{r/2-1}e^{-x/2}}{\Gamma(r/2)2^{r/2}}
//\]
//@@Example
//First, a plot of the PDF for a family of chi-square random variables
//@<
//f = zeros(7,100);
//x = (1:100)/10;
//for n=1:7;t=x.^(n/2-1).*exp(-x/2);f(n,:)=10*t/sum(t);end
//plot(x,f');
//mprint chipdf
//@>
//The PDF is below:
//@figure chipdf
//Here is an example of using @|randchi| and @|randn| to compute
//some chi-square random variables with four degrees of freedom.
//@<
//randchi(4*ones(1,6))
//sum(randn(4,6).^2)
//@>
//!
ArrayVector RandChiFunction(int nargout, const ArrayVector& arg) {
if (arg.size() != 1)
throw Exception("randchi requires exactly one parameter (the vector of degrees of freedom)");
Array arg1(arg[0]);
if (arg1.isEmpty())
return singleArrayVector(Array::emptyConstructor());
arg1.promoteType(FM_FLOAT);
// Output dimension is the larger of the two
Dimensions outDims(arg1.dimensions());
float *dp;
dp = (float *) Malloc(sizeof(float)*outDims.getElementCount());
float *p1;
p1 = (float*) arg1.getDataPointer();
int i;
for (i=0;i<outDims.getElementCount();i++) {
if (p1[i] <= 0)
throw Exception("argument to randchi must be positive");
dp[i] = genchi(p1[i]);
}
return singleArrayVector(Array(FM_FLOAT,outDims,dp));
}
//!
//@Module RANDEXP Generate Exponential Random Variable
//@@Section RANDOM
//@@Usage
//Generates a vector of exponential random variables with
//the specified parameter. The general syntax for its use is
//@[
// y = randexp(lambda)
//@]
//where @|lambda| is a vector containing the parameters for
//the generated random variables.
//@@Function Internals
//The exponential random variable is usually associated with
//the waiting time between events in a Poisson random process.
//The PDF of an exponential random variable is:
//\[
// f(x) = \lambda e^{-\lambda x}
//\]
//@@Example
//Here is an example of using the @|randexp| function to generate
//some exponentially distributed random variables
//@<
//randexp(ones(1,6))
//@>
//!
ArrayVector RandExpFunction(int nargout, const ArrayVector& arg) {
if (arg.size() != 1)
throw Exception("randexp requires exactly one parameter (the vector of means)");
Array arg1(arg[0]);
if (arg1.isEmpty())
return singleArrayVector(Array::emptyConstructor());
arg1.promoteType(FM_FLOAT);
// Output dimension is the larger of the two
Dimensions outDims(arg1.dimensions());
float *dp;
dp = (float *) Malloc(sizeof(float)*outDims.getElementCount());
float *p1;
p1 = (float*) arg1.getDataPointer();
int i;
for (i=0;i<outDims.getElementCount();i++)
dp[i] = genexp(p1[i]);
return singleArrayVector(Array(FM_FLOAT,outDims,dp));
}
//!
//@Module RANDP Generate Poisson Random Variable
//@@Section RANDOM
//@@Usage
//Generates a vector Poisson random variables with the given
//parameters. The general syntax for its use is
//@[
// y = randp(nu),
//@]
//where @|nu| is an array containing the rate parameters
//for the generated random variables.
//@@Function Internals
//A Poisson random variable is generally defined by taking the
//limit of a binomial distribution as the sample size becomes
//large, with the expected number of successes being fixed (so
//that the probability of success decreases as @|1/N|).
//The Poisson distribution is given by
//\[
// P_{\nu}(n) = \frac{\nu^n e^{-nu}}{n!}.
//\]
//@@Example
//Here is an exmaple of using @|randp| to generate some Poisson
//random variables, and also using @|randbin| to do the same
//using @|N=1000| trials to approximate the Poisson result.
//@<
//randp(33*ones(1,10))
//randbin(1000*ones(1,10),33/1000*ones(1,10))
//@>
//!
ArrayVector RandPoissonFunction(int nargout, const ArrayVector& arg) {
if (arg.size() != 1)
throw Exception("randp requires exactly one parameter (the vector of means)");
Array arg1(arg[0]);
if (arg1.isEmpty())
return singleArrayVector(Array::emptyConstructor());
arg1.promoteType(FM_FLOAT);
// Output dimension is the larger of the two
Dimensions outDims(arg1.dimensions());
int32 *dp;
dp = (int32 *) Malloc(sizeof(int32)*outDims.getElementCount());
float *p1;
p1 = (float*) arg1.getDataPointer();
int i;
for (i=0;i<outDims.getElementCount();i++)
dp[i] = ignpoi(p1[i]);
return singleArrayVector(Array(FM_INT32,outDims,dp));
}
//!
//@Module RANDBIN Generate Binomial Random Variables
//@@Section RANDOM
//@@Usage
//Generates random variables with a binomial distribution.
//The general syntax for its use is
//@[
// y = randbin(N,p)
//@]
//where @|N| is a vector representing the number of Bernoulli
//trials, and @|p| is the success probability associated with each
//trial.
//@@Function Internals
//A Binomial random variable describes the number of successful
//outcomes from @|N| Bernoulli trials, with the probability of
//success in each trial being @|p|. The probability distribution
//is
//\[
// P(n) = \frac{N!}{n!(N-n)!}p^n(1-p)^{N-n}
//\]
//@@Example
//Here we generate @|10| binomial random variables, corresponding
//to @|N=100| trials, each with probability @|p=0.1|, using
//both @|randbin| and then again using @|rand| (to simulate the trials):
//@<
//randbin(100,.1*ones(1,10))
//sum(rand(100,10)<0.1)
//@>
//!
ArrayVector RandBinFunction(int nargout, const ArrayVector& arg) {
if (arg.size() != 2)
throw Exception("randbin requires two parameter arguments");
Array arg1(arg[0]);
Array arg2(arg[1]);
if (arg1.isEmpty() || arg2.isEmpty())
return singleArrayVector(Array::emptyConstructor());
// Check the logic to see if one or both are scalar values
if (!(arg1.isScalar() || arg2.isScalar() || (arg1.dimensions().equals(arg2.dimensions()))))
throw Exception("randbin requires either one of the two arguments to be a scalar, or both arguments to be the same size");
int arg1_advance;
int arg2_advance;
arg1_advance = (arg1.isScalar()) ? 0 : 1;
arg2_advance = (arg2.isScalar()) ? 0 : 1;
// Output dimension is the larger of the two
Dimensions outDims;
if (arg1.getLength() > arg2.getLength()) {
outDims = arg1.dimensions();
} else {
outDims = arg2.dimensions();
}
arg1.promoteType(FM_UINT32);
arg2.promoteType(FM_FLOAT);
uint32 *dp;
dp = (uint32 *) Malloc(sizeof(uint32)*outDims.getElementCount());
uint32 *p1;
p1 = (uint32*) arg1.getDataPointer();
float *p2;
p2 = (float*) arg2.getDataPointer();
int i;
for (i=0;i<outDims.getElementCount();i++) {
dp[i] = ignbin(p1[i*arg1_advance],p2[i*arg2_advance]);
}
return singleArrayVector(Array(FM_UINT32,outDims,dp));
}
//!
//@Module RANDNBIN Generate Negative Binomial Random Variables
//@@Section RANDOM
//@@Usage
//Generates random variables with a negative binomial distribution.
//The general syntax for its use is
//@[
// y = randnbin(r,p)
//@]
//where @|r| is a vector of integers representing the number of
//successes, and @|p| is the probability of success.
//@@Function Internals
//A negative binomial random variable describes the number of failures
//@|x| that occur in @|x+r| bernoulli trials, with a success on the
//@|x+r| trial. The pdf is given by
//\[
// P_{r,p}(x)=\left(\begin{matrix}x+r-1\\r-1\end{matrix}\right)p^r(1-p)^x.
//\]
//@@Example
//Here we generate some negative binomial random variables:
//@<
//randnbin(3*ones(1,4),.01)
//randnbin(6*ones(1,4),.01)
//@>
//!
ArrayVector RandNBinFunction(int nargout, const ArrayVector& arg) {
if (arg.size() != 2)
throw Exception("randnbin requires two parameter arguments");
Array arg1(arg[0]);
Array arg2(arg[1]);
if (arg1.isEmpty() || arg2.isEmpty())
return singleArrayVector(Array::emptyConstructor());
// Check the logic to see if one or both are scalar values
if (!(arg1.isScalar() || arg2.isScalar() || (arg1.dimensions().equals(arg2.dimensions()))))
throw Exception("randnbin requires either one of the two arguments to be a scalar, or both arguments to be the same size");
int arg1_advance;
int arg2_advance;
arg1_advance = (arg1.isScalar()) ? 0 : 1;
arg2_advance = (arg2.isScalar()) ? 0 : 1;
// Output dimension is the larger of the two
Dimensions outDims;
if (arg1.getLength() > arg2.getLength()) {
outDims = arg1.dimensions();
} else {
outDims = arg2.dimensions();
}
arg1.promoteType(FM_UINT32);
arg2.promoteType(FM_FLOAT);
uint32 *dp;
dp = (uint32 *) Malloc(sizeof(uint32)*outDims.getElementCount());
uint32 *p1;
p1 = (uint32*) arg1.getDataPointer();
float *p2;
p2 = (float*) arg2.getDataPointer();
int i;
for (i=0;i<outDims.getElementCount();i++) {
dp[i] = ignnbn(p1[i*arg1_advance],p2[i*arg2_advance]);
}
return singleArrayVector(Array(FM_UINT32,outDims,dp));
}
//!
//@Module RANDF Generate F-Distributed Random Variable
//@@Section RANDOM
//@@Usage
//Generates random variables with an F-distribution. The general
//syntax for its use is
//@[
// y = randf(n,m)
//@]
//where @|n| and @|m| are vectors of the number of degrees of freedom
//in the numerator and denominator of the chi-square random variables
//whose ratio defines the statistic.
//@@Function Internals
//The statistic @|F_{n,m}| is defined as the ratio of two chi-square
//random variables:
//\[
// F_{n,m} = \frac{\chi_n^2/n}{\chi_m^2/m}
//\]
//The PDF is given by
//\[
// f_{n,m} = \frac{m^{m/2}n^{n/2}x^{n/2-1}}{(m+nx)^{(n+m)/2}B(n/2,m/2)},
//\]
//where @|B(a,b)| is the beta function.
//@@Example
//Here we use @|randf| to generate some F-distributed random variables,
//and then again using the @|randchi| function:
//@<
//randf(5*ones(1,9),7)
//randchi(5*ones(1,9))./randchi(7*ones(1,9))
//@>
//!
ArrayVector RandFFunction(int nargout, const ArrayVector& arg) {
if (arg.size() != 2)
throw Exception("randf requires two parameter arguments");
Array arg1(arg[0]);
Array arg2(arg[1]);
if (arg1.isEmpty() || arg2.isEmpty())
return singleArrayVector(Array::emptyConstructor());
// Check the logic to see if one or both are scalar values
if (!(arg1.isScalar() || arg2.isScalar() || (arg1.dimensions().equals(arg2.dimensions()))))
throw Exception("randf requires either one of the two arguments to be a scalar, or both arguments to be the same size");
int arg1_advance;
int arg2_advance;
arg1_advance = (arg1.isScalar()) ? 0 : 1;
arg2_advance = (arg2.isScalar()) ? 0 : 1;
// Output dimension is the larger of the two
Dimensions outDims;
if (arg1.getLength() > arg2.getLength()) {
outDims = arg1.dimensions();
} else {
outDims = arg2.dimensions();
}
arg1.promoteType(FM_FLOAT);
arg2.promoteType(FM_FLOAT);
float *dp;
dp = (float *) Malloc(sizeof(float)*outDims.getElementCount());
float *p1;
p1 = (float*) arg1.getDataPointer();
float *p2;
p2 = (float*) arg2.getDataPointer();
int i;
for (i=0;i<outDims.getElementCount();i++) {
if ((p1[i*arg1_advance] <= 0) || (p2[i*arg2_advance]) <= 0)
throw Exception("randf requires positive arguments");
dp[i] = genf(p1[i*arg1_advance],p2[i*arg2_advance]);
}
return singleArrayVector(Array(FM_FLOAT,outDims,dp));
}
//!
//@Module RANDGAMMA Generate Gamma-Distributed Random Variable
//@@Section RANDOM
//@@Usage
//Generates random variables with a gamma distribution. The general
//syntax for its use is
//@[
// y = randgamma(a,r),
//@]
//where @|a| and @|r| are vectors describing the parameters of the
//gamma distribution. Roughly speaking, if @|a| is the mean time between
//changes of a Poisson random process, and we wait for the @|r| change,
//the resulting wait time is Gamma distributed with parameters @|a|
//and @|r|.
//@@Function Internals
//The Gamma distribution arises in Poisson random processes. It represents
//the waiting time to the occurance of the @|r|-th event in a process with
//mean time @|a| between events. The probability distribution of a Gamma
//random variable is
//\[
// P(x) = \frac{a^r x^{r-1} e^{-ax}}{\Gamma(r)}.
//\]
//Note also that for integer values of @|r| that a Gamma random variable
//is effectively the sum of @|r| exponential random variables with parameter
//@|a|.
//@@Example
//Here we use the @|randgamma| function to generate Gamma-distributed
//random variables, and then generate them again using the @|randexp|
//function.
//@<
//randgamma(1,15*ones(1,9))
//sum(randexp(ones(15,9)))
//@>
//!
ArrayVector RandGammaFunction(int nargout, const ArrayVector& arg) {
if (arg.size() != 2)
throw Exception("randgamma requires two parameter arguments");
Array arg1(arg[0]);
Array arg2(arg[1]);
if (arg1.isEmpty() || arg2.isEmpty())
return singleArrayVector(Array::emptyConstructor());
// Check the logic to see if one or both are scalar values
if (!(arg1.isScalar() || arg2.isScalar() || (arg1.dimensions().equals(arg2.dimensions()))))
throw Exception("randgamma requires either one of the two arguments to be a scalar, or both arguments to be the same size");
int arg1_advance;
int arg2_advance;
arg1_advance = (arg1.isScalar()) ? 0 : 1;
arg2_advance = (arg2.isScalar()) ? 0 : 1;
// Output dimension is the larger of the two
Dimensions outDims;
if (arg1.getLength() > arg2.getLength()) {
outDims = arg1.dimensions();
} else {
outDims = arg2.dimensions();
}
arg1.promoteType(FM_FLOAT);
arg2.promoteType(FM_FLOAT);
float *dp;
dp = (float *) Malloc(sizeof(float)*outDims.getElementCount());
float *p1;
p1 = (float*) arg1.getDataPointer();
float *p2;
p2 = (float*) arg2.getDataPointer();
int i;
for (i=0;i<outDims.getElementCount();i++) {
dp[i] = gengam(p1[i*arg1_advance],p2[i*arg2_advance]);
}
return singleArrayVector(Array(FM_FLOAT,outDims,dp));
}
//!
//@Module RANDMULTI Generate Multinomial-distributed Random Variables
//@@Section RANDOM
//@@Usage
//This function generates samples from a multinomial distribution
//given the probability of each outcome. The general syntax for
//its use is
//@[
// y = randmulti(N,pvec)
//@]
//where @|N| is the number of experiments to perform, and @|pvec|
//is the vector of probabilities describing the distribution of
//outcomes.
//@@Function Internals
//A multinomial distribution describes the number of times each
//of @|m| possible outcomes occurs out of @|N| trials, where each
//outcome has a probability @|p_i|. More generally, suppose that
//the probability of a Bernoulli random variable @|X_i| is @|p_i|,
//and that
//\[
// \sum_{i=1}^{m} p_i = 1.
//\]
//Then the probability that @|X_i| occurs @|x_i| times is
//\[
// P_N(x_1,x_2,\ldots,x_n) = \frac{N!}{x_1!\cdots x_n!} p_1^{x_1}\cdots p_n^{x_n}.
//\]
//@@Example
//Suppose an experiment has three possible outcomes, say heads,
//tails and edge, with probabilities @|0.4999|, @|0.4999| and
//@|0.0002|, respectively. Then if we perform ten thousand coin
//flips we get
//@<
//randmulti(10000,[0.4999,0.4999,0.0002])
//@>
//!
ArrayVector RandMultiFunction(int nargout, const ArrayVector& arg) {
if (arg.size() != 2)
throw Exception("randmulti requires two parameter arguments");
Array arg1(arg[0]);
Array arg2(arg[1]);
if (arg1.isEmpty() || arg2.isEmpty())
return singleArrayVector(Array::emptyConstructor());
int N = arg1.getContentsAsIntegerScalar();
if (N<0)
throw Exception("number of events to generate for randmulti must be a nonnegative integer");
arg2.promoteType(FM_FLOAT);
// Verify the correctness of the probability argument
float *dp;
dp = (float*) arg2.getReadWriteDataPointer();
float Psum = 0.0;
int i;
for (i=0;i<arg2.getLength();i++) {
if ((dp[i] < 0) || (dp[i] > 1))
throw Exception("probabiliy vector argument to randmulti must have all elements between 0 and 1");
Psum += dp[i];
}
for (i=0;i<arg2.getLength();i++) {
dp[i] /= Psum;
}
Dimensions outDims;
outDims = arg2.dimensions();
int32 *ip = (int32*) Malloc(sizeof(int32)*arg2.getLength());
genmul(N,dp,arg2.getLength(),(long int*) ip);
return singleArrayVector(Array(FM_INT32,outDims,ip));
}
//!
//@Module RANDNCHI Generate Noncentral Chi-Square Random Variable
//@@Section RANDOM
//@@Usage
//Generates a vector of non-central chi-square random variables
//with the given number of degrees of freedom and the given
//non-centrality parameters. The general syntax for its use is
//@[
// y = randnchi(n,mu)
//@]
//where @|n| is an array containing the degrees of freedom for
//each generated random variable (with each element of @|n| >= 1),
//and @|mu| is the non-centrality shift (must be positive).
//@@Function Internals
//A non-central chi-square random variable is the sum of a chisquare
//deviate with @|n-1| degrees of freedom plus the square of a normal
//deviate with mean @|mu| and standard deviation 1.
//@@Examples
//Here is an example of a non-central chi-square random variable:
//@<
//randnchi(5*ones(1,9),0.3)
//@>
//!
ArrayVector RandNChiFunction(int nargout, const ArrayVector& arg) {
if (arg.size() != 2)
throw Exception("randnchi requires two parameter arguments");
Array arg1(arg[0]);
Array arg2(arg[1]);
if (arg1.isEmpty() || arg2.isEmpty())
return singleArrayVector(Array::emptyConstructor());
// Check the logic to see if one or both are scalar values
if (!(arg1.isScalar() || arg2.isScalar() || (arg1.dimensions().equals(arg2.dimensions()))))
throw Exception("randnchi requires either one of the two arguments to be a scalar, or both arguments to be the same size");
int arg1_advance;
int arg2_advance;
arg1_advance = (arg1.isScalar()) ? 0 : 1;
arg2_advance = (arg2.isScalar()) ? 0 : 1;
// Output dimension is the larger of the two
Dimensions outDims;
if (arg1.getLength() > arg2.getLength()) {
outDims = arg1.dimensions();
} else {
outDims = arg2.dimensions();
}
arg1.promoteType(FM_FLOAT);
arg2.promoteType(FM_FLOAT);
float *dp;
dp = (float *) Malloc(sizeof(float)*outDims.getElementCount());
float *p1;
p1 = (float*) arg1.getDataPointer();
float *p2;
p2 = (float*) arg2.getDataPointer();
int i;
for (i=0;i<outDims.getElementCount();i++) {
if (p1[i*arg1_advance] <= 1.0)
throw Exception("degrees of freedom argument must be > 1.0");
if (p2[i*arg2_advance] < 0.0)
throw Exception("noncentrality parameter must be positive");
dp[i] = gennch(p1[i*arg1_advance],p2[i*arg2_advance]);
}
return singleArrayVector(Array(FM_FLOAT,outDims,dp));
}
//!
//@Module RANDNF Generate Noncentral F-Distribution Random Variable
//@@Section RANDOM
//@@Usage
//Generates a vector of non-central F-distributed random variables
//with the specified parameters. The general syntax for its use is
//@[
// y = randnf(n,m,c)
//@]
//where @|n| is the number of degrees of freedom in the numerator,
//and @|m| is the number of degrees of freedom in the denominator.
//The vector @|c| determines the non-centrality shift of the numerator.
//@@Function Internals
//A non-central F-distributed random variable is the ratio of a
//non-central chi-square random variable and a central chi-square random
//variable, i.e.,
//\[
// F_{n,m,c} = \frac{\chi_{n,c}^2/n}{\chi_m^2/m}.
//\]
//@@Example
//Here we use the @|randf| to generate some non-central F-distributed
//random variables:
//@<
//randnf(5*ones(1,9),7,1.34)
//@>
//!
ArrayVector RandNFFunction(int nargout, const ArrayVector& arg) {
if (arg.size() != 3)
throw Exception("randnf requires three parameter arguments");
Array arg1(arg[0]);
Array arg2(arg[1]);
Array arg3(arg[2]);
if (arg1.isEmpty() || arg2.isEmpty() || arg3.isEmpty())
return singleArrayVector(Array::emptyConstructor());
int arg1_advance;
int arg2_advance;
int arg3_advance;
arg1_advance = (arg1.isScalar()) ? 0 : 1;
arg2_advance = (arg2.isScalar()) ? 0 : 1;
arg3_advance = (arg3.isScalar()) ? 0 : 1;
if (arg1_advance && arg2_advance && (!arg1.dimensions().equals(arg2.dimensions())))
throw Exception("vector arguments to randnf must be the same size");
if (arg1_advance && arg3_advance && (!arg1.dimensions().equals(arg3.dimensions())))
throw Exception("vector arguments to randnf must be the same size");
if (arg2_advance && arg3_advance && (!arg2.dimensions().equals(arg3.dimensions())))
throw Exception("vector arguments to randnf must be the same size");
// Output dimension is the larger of the two
Dimensions outDims;
if ((arg1.getLength() > arg2.getLength()) && (arg1.getLength() > arg3.getLength()))
outDims = arg1.dimensions();
else if ((arg2.getLength() > arg1.getLength()) && (arg2.getLength() > arg3.getLength()))
outDims = arg2.dimensions();
else if ((arg3.getLength() > arg1.getLength()) && (arg3.getLength() > arg2.getLength()))
outDims = arg3.dimensions();
arg1.promoteType(FM_FLOAT);
arg2.promoteType(FM_FLOAT);
arg3.promoteType(FM_FLOAT);
float *dp;
dp = (float *) Malloc(sizeof(float)*outDims.getElementCount());
float *p1;
p1 = (float*) arg1.getDataPointer();
float *p2;
p2 = (float*) arg2.getDataPointer();
float *p3;
p3 = (float*) arg3.getDataPointer();
int i;
for (i=0;i<outDims.getElementCount();i++) {
if (p1[i*arg1_advance] <= 1.0)
throw Exception("numerator degrees of freedom argument must be > 1.0");
if (p2[i*arg2_advance] <= 0.0)
throw Exception("denominator degrees of freedom must be positive");
if (p3[i*arg2_advance] < 0.0)
throw Exception("noncentrality parameter must be non-negative");
dp[i] = gennf(p1[i*arg1_advance],p2[i*arg2_advance],p3[i*arg3_advance]);
}
return singleArrayVector(Array(FM_FLOAT,outDims,dp));
}
void InitializeRandGen() {
unsigned long init[4]={0x923, 0x234, 0x405, 0x456}, length=4;
init_by_array(init, length);
initialized = true;
}
ArrayVector RandStateControl(const ArrayVector& arg) {
string key = arg[0].getContentsAsStringUpper();
if (!(key=="STATE"))
throw Exception("expecting string 'state' as first argument");
if (arg.size() == 1) {
uint32 *mt = (uint32*) Malloc(sizeof(uint32)*625);
GetRandStateVect(mt);
return singleArrayVector(Array(FM_UINT32,Dimensions(625,1),mt));
} else {
Array statevec(arg[1]);
if ((statevec.isScalar()) && (statevec.getContentsAsIntegerScalar() == 0)) {
InitializeRandGen();
return ArrayVector();
} else {
statevec.promoteType(FM_UINT32);
if (statevec.getLength() != 625)
throw Exception("illegal state vector - must be of length 625");
SetRandStateVect((uint32*)statevec.getDataPointer());
return ArrayVector();
}
}
}
//!
//@Module RANDN Gaussian (Normal) Random Number Generator
//@@Section RANDOM
//@@Usage
//Creates an array of pseudo-random numbers of the specified size.
//The numbers are normally distributed with zero mean and a unit
//standard deviation (i.e., @|mu = 0, sigma = 1|).
// Two seperate syntaxes are possible. The first syntax specifies the array
//dimensions as a sequence of scalar dimensions:
//@[
// y = randn(d1,d2,...,dn).
//@]
//The resulting array has the given dimensions, and is filled with
//random numbers. The type of @|y| is @|double|, a 64-bit floating
//point array. To get arrays of other types, use the typecast
//functions.
//
//The second syntax specifies the array dimensions as a vector,
//where each element in the vector specifies a dimension length:
//@[
// y = randn([d1,d2,...,dn]).
//@]
//This syntax is more convenient for calling @|randn| using a
//variable for the argument.
//
//Finally, @|randn| supports two additional forms that allow
//you to manipulate the state of the random number generator.
//The first retrieves the state
//@[
// y = randn('state')
//@]
//which is a 625 length integer vector. The second form sets
//the state
//@[
// randn('state',y)
//@]
//or alternately, you can reset the random number generator with
//@[
// randn('state',0)
//@]
//@@Function Internals
//Recall that the
//probability density function (PDF) of a normal random variable is
//\[
//f(x) = \frac{1}{\sqrt{2\pi \sigma^2}} e^{\frac{-(x-\mu)^2}{2\sigma^2}}.
//\]
//The Gaussian random numbers are generated from pairs of uniform random numbers using a transformation technique.
//@@Example
//The following example demonstrates an example of using the first form of the @|randn| function.
//@<
//randn(2,2,2)
//@>
//The second example demonstrates the second form of the @|randn| function.
//@<
//randn([2,2,2])
//@>
//In the next example, we create a large array of 10000 normally distributed pseudo-random numbers. We then shift the mean to 10, and the variance to 5. We then numerically calculate the mean and variance using @|mean| and @|var|, respectively.
//@<
//x = 10+sqrt(5)*randn(1,10000);
//mean(x)
//var(x)
//@>
//Now, we use the state manipulation functions of @|randn| to exactly reproduce
//a random sequence. Note that unlike using @|seed|, we can exactly control where
//the random number generator starts by saving the state.
//@<
//randn('state',0) % restores us to startup conditions
//a = randn(1,3) % random sequence 1
//b = randn('state'); % capture the state vector
//c = randn(1,3) % random sequence 2
//randn('state',b); % restart the random generator so...
//c = randn(1,3) % we get random sequence 2 again
//@>
//!
ArrayVector RandnFunction(int nargout, const ArrayVector& arg) {
if (!initialized)
InitializeRandGen();
int i;
Array t;
Dimensions dims;
int32 *dp;
if (arg.size() == 0)
dims.makeScalar();
else {
if (arg[0].isString())
return RandStateControl(arg);
// Case 1 - all of the entries are scalar
bool allScalars;
allScalars = true;
for (i=0;i<arg.size();i++)
allScalars &= arg[i].isScalar();
if (allScalars) {
t = arg[0];
if (arg.size() == 1) {
// If all scalars and only one argument - we want a square zero matrix
dims.set(0,t.getContentsAsIntegerScalar());
dims.set(1,dims.get(0));
} else {
// If all scalars and and multiple arguments, we count dimensions
for (i=0;i<arg.size();i++) {
t = arg[i];
dims.set(i,t.getContentsAsIntegerScalar());
}
}
} else {
if (arg.size() > 1)
throw Exception("Arguments to randn function must be either all scalars or a single vector");
t = arg[0];
t.promoteType(FM_UINT32);
dp = (int*) t.getDataPointer();
for (i=0;i<t.getLength();i++)
dims.set(i,dp[i]);
}
bool allPositive;
allPositive = true;
for (i=0;i<dims.getLength();i++)
allPositive &= (dims.get(i) >= 0);
if (!allPositive)
throw Exception("Randn function requires positive arguments");
}
double *qp;
qp = (double*) Malloc(sizeof(double)*dims.getElementCount());
int len;
int j;
len = dims.getElementCount();
for (j=0;j<len;j+=2) {
double x1, x2, w, y1, y2;
do {
x1 = 2.0*genrand_res53()-1.0;
x2 = 2.0*genrand_res53()-1.0;
w = x1 * x1 + x2 * x2;
} while ( w >= 1.0 );
w = sqrt( (-2.0 * log( w ) ) / w );
y1 = x1 * w;
y2 = x2 * w;
qp[j] = y1;
if (j<len-1)
qp[j+1] = y2;
}
return singleArrayVector(Array(FM_DOUBLE,dims,qp));
}
//!
//@Module RAND Uniform Random Number Generator
//@@Section RANDOM
//@@Usage
//Creates an array of pseudo-random numbers of the specified size.
//The numbers are uniformly distributed on @|[0,1)|.
//Two seperate syntaxes are possible. The first syntax specifies the array
//dimensions as a sequence of scalar dimensions:
//@[
// y = rand(d1,d2,...,dn).
//@]
//The resulting array has the given dimensions, and is filled with
//random numbers. The type of @|y| is @|double|, a 64-bit floating
//point array. To get arrays of other types, use the typecast
//functions.
//
//The second syntax specifies the array dimensions as a vector,
//where each element in the vector specifies a dimension length:
//@[
// y = rand([d1,d2,...,dn]).
//@]
//This syntax is more convenient for calling @|rand| using a
//variable for the argument.
//
//Finally, @|rand| supports two additional forms that allow
//you to manipulate the state of the random number generator.
//The first retrieves the state
//@[
// y = rand('state')
//@]
//which is a 625 length integer vector. The second form sets
//the state
//@[
// rand('state',y)
//@]
//or alternately, you can reset the random number generator with
//@[
// rand('state',0)
//@]
//@@Example
//The following example demonstrates an example of using the first form of the @|rand| function.
//@<
//rand(2,2,2)
//@>
//The second example demonstrates the second form of the @|rand| function.
//@<
//rand([2,2,2])
//@>
//The third example computes the mean and variance of a large number of uniform random numbers. Recall that the mean should be @|1/2|, and the variance should be @|1/12 ~ 0.083|.
//@<
//x = rand(1,10000);
//mean(x)
//var(x)
//@>
//Now, we use the state manipulation functions of @|rand| to exactly reproduce
//a random sequence. Note that unlike using @|seed|, we can exactly control where
//the random number generator starts by saving the state.
//@<
//rand('state',0) % restores us to startup conditions
//a = rand(1,3) % random sequence 1
//b = rand('state'); % capture the state vector
//c = rand(1,3) % random sequence 2
//rand('state',b); % restart the random generator so...
//c = rand(1,3) % we get random sequence 2 again
//@>
//!
ArrayVector RandFunction(int nargout, const ArrayVector& arg) {
if (!initialized)
InitializeRandGen();
int i;
Array t;
Dimensions dims;
int32 *dp;
if (arg.size() == 0)
dims.makeScalar();
else {
// Check for state assignment
if (arg[0].isString())
return RandStateControl(arg);
// Case 1 - all of the entries are scalar
bool allScalars;
allScalars = true;
for (i=0;i<arg.size();i++)
allScalars &= arg[i].isScalar();
if (allScalars) {
t = arg[0];
if (arg.size() == 1) {
// If all scalars and only one argument - we want a square zero matrix
dims.set(0,t.getContentsAsIntegerScalar());
dims.set(1,dims.get(0));
} else {
// If all scalars and and multiple arguments, we count dimensions
for (i=0;i<arg.size();i++) {
t = arg[i];
dims.set(i,t.getContentsAsIntegerScalar());
}
}
} else {
if (arg.size() > 1)
throw Exception("Arguments to rand function must be either all scalars or a single vector");
t = arg[0];
t.promoteType(FM_UINT32);
dp = (int*) t.getDataPointer();
for (i=0;i<t.getLength();i++)
dims.set(i,dp[i]);
}
bool allPositive;
allPositive = true;
for (i=0;i<dims.getLength();i++)
allPositive &= (dims.get(i) >= 0);
if (!allPositive)
throw Exception("Rand function requires posItive arguments");
}
double *qp;
qp = (double*) Malloc(sizeof(double)*dims.getElementCount());
int len;
int j;
len = dims.getElementCount();
for (j=0;j<len;j++)
qp[j] = genrand_res53();
return singleArrayVector(Array(FM_DOUBLE,dims,qp));
}
syntax highlighted by Code2HTML, v. 0.9.1