| Subject: | Function random_unit_vector Not Uniform; Weighted Towards First Quadrant. [with patch] |
Regarding Random-PoissonDisc-0.01:
The function Random::PoissonDisc::random_unit_vector does not return uniformly distributed unit vectors. The returned vectors are biased towards the first quadrant (for 2 dimensions).
A test script and a patch is supplied.
I marked this bug as critical since the correctness of the main function of the module depend on the distribution of the generated vectors.
==========(output of attached test script)
Perl version is v5.18.1
Random::PoissonDisc version: v0.01
Histogram (nominal 10000 per bin) [range: count; error]:
[0.0000 .. 0.3142): 17196; +71.96%
[0.3142 .. 0.6283): 19530; +95.30%
[0.6283 .. 0.9425): 21383; +113.83%
[0.9425 .. 1.2566): 20068; +100.68%
[1.2566 .. 1.5708): 17232; +72.32%
[1.5708 .. 1.8850): 13535; +35.35%
[1.8850 .. 2.1991): 10449; +4.49%
[2.1991 .. 2.5133): 7750; -22.50%
[2.5133 .. 2.8274): 6102; -38.98%
[2.8274 .. 3.1416): 4745; -52.55%
[3.1416 .. 3.4558): 4238; -57.62%
[3.4558 .. 3.7699): 3702; -62.98%
[3.7699 .. 4.0841): 3502; -64.98%
[4.0841 .. 4.3982): 3520; -64.80%
[4.3982 .. 4.7124): 4164; -58.36%
[4.7124 .. 5.0265): 4822; -51.78%
[5.0265 .. 5.3407): 6048; -39.52%
[5.3407 .. 5.6549): 7909; -20.91%
[5.6549 .. 5.9690): 10290; +2.90%
[5.9690 .. 6.2832): 13815; +38.15%
==========(end script output)
The output above shows the distribution of the angle of returned 2d vectors (wrapped to the interval "[0 .. 2*pi)" ). Notice that it is biased towards pi/4 (=0.79, middle of first quadrant) and away from 5*pi/4 (=3.93, middle of third quadrant).
This behaviour is due to incorrect scaling/offset of the gaussian distributed vector numbers (before normalizing). The vectors are generated by:
@vec = map { 1-2*gaussian() } 1..$dimensions;
However, gaussian() returns a number from a normal distribution with mean 0 and standard deviation 1. Thus, the scaling/offset will in effect give a distribution with mean 1 and standard deviation 2.
For correct generation of uniformly distributed (w.r.t. angle) unit vectors the coordinates (before normalizing) should be sampled from a normal distribution with mean 0. (Any value can be used for the standard deviation, since we are normalizing anyway.)
The correct code row should thus be:
@vec = map { gaussian() } 1..$dimensions;
Making this change and running the test script again gives errors of max +/- a few percent for all bins. This is consistent with uniformly distributed angles.
Other information:
$ uname -a
Linux johan-901 3.5.0-41-generic #64-Ubuntu SMP Wed Sep 11 15:40:48 UTC 2013 i686 i686 i686 GNU/Linux
| Subject: | poissondisc-randvect-histogram.pl |
#!/usr/bin/env perl
use common::sense;
use POSIX qw( floor ceil log10 );
use Math::Trig qw( pi2 );
use Random::PoissonDisc;
printf "Perl version is v%vd\n", $^V;
say "Random::PoissonDisc version: v$Random::PoissonDisc::VERSION";
use constant HISTSIZE => 20;
use constant HISTW => pi2;
use constant BIN_NOM_COUNT => 10000;
use constant SAMPLES => HISTSIZE * BIN_NOM_COUNT;
sub rand_angle {
# use this instead for testing the histogram
#return rand pi2;
my ($x,$y) = @{Random::PoissonDisc::random_unit_vector(2)};
return atan2($y,$x);
}
my @hist = ( (0) x HISTSIZE );
for my $i (1..SAMPLES) {
my $a = rand_angle;
$a /= pi2;
$a -= floor $a;
$a *= HISTSIZE;
$hist[int $a]++;
}
sub bin_range ($) {map {pi2 * $_ / HISTSIZE} ($_[0], $_[0]+1)}
sub bin_error_pct ($) {scalar 100 * (($_[0]/BIN_NOM_COUNT) - 1 )}
my $w = ceil(log10(BIN_NOM_COUNT))+1;
say "Histogram (nominal ", BIN_NOM_COUNT, " per bin) [range: count; error]:";
while (my ($i, $n) = each @hist) {
printf "[%6.4f .. %6.4f): %${w}d; %+7.2f%%\n", bin_range($i), $n, bin_error_pct $n;
}
| Subject: | poissondisc-randvect.patch |
diff --git i/lib/Random/PoissonDisc.pm w/lib/Random/PoissonDisc.pm
index e3b5902..7577422 100644
--- i/lib/Random/PoissonDisc.pm
+++ w/lib/Random/PoissonDisc.pm
@@ -276,9 +276,8 @@ sub random_unit_vector {
my (@vec,$len);
# Create normal distributed coordinates
- # around 0 in (-1,1)
RETRY: {
- @vec = map { 1-2*gaussian() } 1..$dimensions;
+ @vec = map { gaussian() } 1..$dimensions;
$len = norm(@vec);
redo RETRY unless $len;
};