File : randdist.sp


$ spar randdist
aleph  2.00260000000000E-01 2.00000000000000E-01
beth  1.66376000000000E-01 1.66666666666667E-01
gimel  1.42698000000000E-01 1.42857142857143E-01
daleth  1.25408000000000E-01 1.25000000000000E-01
he  1.11311000000000E-01 1.11111111111111E-01
waw  9.97980000000000E-02 1.00000000000000E-01
zayin  9.09570000000000E-02 9.09090909090909E-02
heth  6.31920000000000E-02 rest


#!/usr/local/bin/spar

pragma annotate( summary, "randdist" )
              @( description, "Given a mapping between items and their required" )
              @( description, "probability of occurrence, generate a million items" )
              @( description, "randomly subject to the given probabilities and compare" )
              @( description, "the target probability of occurrence versus the" )
              @( description, "generated values." )
              @( description, "" )
              @( description, "The total of all the probabilities should equal one." )
              @( description, "(Because floating point arithmetic is involved this is" )
              @( description, "subject to rounding errors).  Use the following mapping" )
              @( description, "to test your programs: aleph 1/5.0, beth 1/6.0," )
              @( description, "gimel 1/7.0, daleth 1/8.0, he 1/9.0,  waw 1/10.0" )
              @( description, "zayin 1/11.0, heth 1759/27720 adjusted so that" )
              @( description, "probabilities add to 1" )
              @( see_also, "http://rosettacode.org/wiki/Probabilistic_choice" )
              @( author, "Ken O. Burtch" );
pragma license( unrestricted );

pragma restriction( no_external_commands );

procedure randdist is
  trials : constant positive := 1_000_000;
  type outcome is (aleph, beth, gimel, daleth, he, waw, zayin, heth);
  pr : constant array(aleph..heth) of float :=
     (1/5, 1/6, 1/7, 1/8, 1/9, 1/10, 1/11, 1 );
  samples : array(aleph..heth) of natural := (0, 0, 0, 0, 0, 0, 0, 0);
  random_value : float;
begin
  for try in 1..trials loop
    random_value := numerics.random;
    for i in arrays.first( pr )..arrays.last( pr ) loop
       if random_value <= pr(i) then
         samples(i) := samples(i) + 1;
         exit;
       else
         random_value := @ - pr(i);
       end if;
    end loop;
  end loop;
  -- Show results
  for i in arrays.first( pr )..arrays.last( pr ) loop
    put( i ) @ ( " " ) @ ( float( samples( i ) ) / float( trials ) );
    if i = heth then
       put_line( " rest" );
    else
       put_line( pr(i) );
    end if;
  end loop;
end randdist;

-- VIM editor formatting instructions
-- vim: ft=spar