.
#======================================================================
# L I B R A N D . P L
# doc: Thu Nov 19 14:27:19 2015
# dlm: Wed Sep 6 10:18:42 2017
# (c) 2015 A.M. Thurnherr
# uE-Info: 17 61 NIL 0 0 70 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
# Nov 19, 2015: - created
# Sep 6, 2017: - finally implemented gauss_rand()
#----------------------------------------------------------------------------------------------------
# From info found at [http://www.design.caltech.edu/erik/Misc/Gaussian.html]
#
# verified with:
# plot '<Cat -Lrand -f =1,1,1e5 -F r=gauss_rand() | Hist r
#----------------------------------------------------------------------------------------------------
{ my($cached); # NB: cached values is normalized
sub gauss_rand(@)
{
my($mu,$sigma) = &antsFunUsage(-1,'ff','[mu[,sigma]]',@_)
if (@_);
$sigma = 1 unless defined($sigma);
$mu = 0 unless defined($mu);
if (defined($cached)) {
my($Y) = $cached * $sigma + $mu;
undef($cached);
return $Y;
}
my($X1,$X2,,$w);
do {
$X1 = 2*rand() - 1;
$X2 = 2*rand() - 1;
$w = $X1**2 + $X2**2;
} while ($w >= 1);
$w = sqrt((-2 * log($w)) / $w);
my($cached) = $X2 * $w;
return $X1 * $w * $sigma + $mu;
}
}
#----------------------------------------------------------------------------------------------------
# From info found at [http://www.mathworks.com/matlabcentral/newsreader/view_thread/301276]
#
# verified with:
# plot '<Cat -Lrand -f =1,1,1e5 -F r=pwrlaw_rand(-2) | Hist -s 100 r | Cat -S $2>2' lt 3,x**-2*1e7
# plot '<Cat -Lrand -f =1,1,1e5 -F r=pwrlaw_rand(-3) | Hist r',x**-3*7e3
# plot '<Cat -Lrand -f =1,1,1e5 -F r=pwrlaw_rand(0) | Hist -s 0.01 r'
# plot '<Cat -Lrand -f =1,1,1e5 -F r=pwrlaw_rand(1) | Hist -s 0.01 r'
# plot '<Cat -Lrand -f =1,1,1e5 -F r=pwrlaw_rand(2) | Hist -s 0.01 r'
#----------------------------------------------------------------------------------------------------
sub pwrlaw_rand($)
{
my($p) = &antsFunUsage(1,'f','exponent',@_);
return rand() ** (1/($p+1));
}
1;