various weighted random number algorithms in perl

random_weights.pl

#!/usr/bin/perl

=comment
  various attempts at getting weighted random number distributions
  
  first,second are  dice-roll convergences 
  (if you roll multiple die multiple times the results will converge on a normal curve)
  http://en.wikipedia.org/wiki/Central_limit_theorem
  boxmuller gives a normal/gaussian distribution (classic bell curve)
  rayleigh and weibull have more variables to play with
  they both tend to weight the peak to the lower part of the range, which is what we are after
  (also try using the log result of boxmuller)

=cut

use Math::Trig;
$samples = 10000;
$max=100;
$min=0;
srand();

boxmuller();
output();

sub weibull {
       print "\n weibull \n ";
 		$alpha = 2.5;
	   $beta = 100.0;
	   $scale = 100.0;
	   print "alpha $alpha beta $beta scale $scale\n";
	for ($n=1;$n<$samples;$n++){
	   $U = rand(1);
	   #$X = -log(1-$U);
	   $X=$scale*((-1/$alpha)*log(1-$U))^(1/$beta);
	   print " $X ";
	  #unless ($X > $max) {
			push (@results, int($X));
		#}
	}
}

sub rayleigh {
		$sigma = 0.5;

       print "\n rayleigh sigma:$sigma \n ";
		for ($n=1;$n<$samples;$n++){
			$U1 = rand(1);
			$f1 = $sigma * (sqrt(-2 * log($U1)) );
			print " $f1 ";
	#		unless ($f1 > $max) {
				push (@results, int($f1*100));
	#			}
		}
}
sub boxmuller {
       print "\n boxmuller \n ";

$scale = 20;
	for ($n=1;$n<$samples;$n++){
	$W1 = 2;
  	while ( $W1>1 or $W1 eq 0) {
		$U1 = 2*rand(1)-1;
		$U2 = 2*rand(1)-1;

		$W1 = $U1*$U1+$U2*$U2;
		}
		
	$W2 = sqrt((-2*log($W1))/$W1);

	$Y1 = int($U1*$W2*$scale);
	$Y2 = int($U2*$W2*$scale);
	#$Y1 = ($Y1 eq 0)?0:log ($Y1*$Y1);
	#$Y2 = ($Y2 eq 0)?0:log ($Y2*$Y2);
	 print "$Y1 ";
	push (@results,$Y1); 
	print "$Y2 ";
	push (@results,$Y2); 
}
}

sub second {
       print "\n second \n ";

	for ($n=1;$n<$samples;$n++){
		 $result = 0;
		 $divisor = 0.0;
		 $randomInt = 0;
		 $sumOfRandoms = 0;
		 $diff = ($max-$min);
		for ($idx=$diff; $idx>0 ; $idx--) {
			$randomInt = rand($diff);
			$sumOfRandoms += $randomInt;
		}
		if ($diff>0) {
		$result = int($sumOfRandoms/$diff);
		}
		 push (@results, $result);
		printf (" %d ",$result);
	}
}

sub first {
       print "\n first \n ";

	for ($n=1;$n<$samples;$n++){
		 $result = 0;
		 $divisor = 0.0;
		 $randomInt = 0;
		 $sum = 0;
		 $diff = ($max-$min);
		for ($idx=$diff; $idx>0 ; $idx--) {
			$randomInt += rand($diff);
			$sum += $idx;
		}
		if ($sum>0 && $diff>0) {
			$divisor = $sum/$diff;
		   
		$result = $min+int($randomInt/$divisor);
		 }
		 push (@results, $result);
		printf (" %d ",$result);
	}
}

sub output {
   # @results = @$_[0];
	@sorted = reverse sort {$a <=> $b} @results;
	$oldval=0;
	while (@sorted) {
		$newval = pop(@sorted);
		if ($newval ne $oldval) {
			printf ("\n %2.d ",$newval);
		}
		print "*";
		$oldval = $newval;
		$newval = 0;
	}
}