Re: Discrete random with given distribution ?

Re: Discrete random with given distribution ?

[Gautier direct_replies_not_read]

tmoran:

>Why not just use your second method (summing probabilities till
>you hit a certain level) to calculate, at run time, the population
>of a lookup table, then use the first method to speedily index into
>it to pick a random element.

And why not replace the loop by the traversal of a binary
tree ? (I just don't like that lookup table...)
________________________________________________________
Gautier  --  http://www.mysunrise.ch/users/gdm/gsoft.htm

NB: For a direct answer, address on the Web site!

_________________________________________________________________
Chat with friends online, try MSN Messenger: http://messenger.msn.com

Re: Discrete random with given distribution ?

[tmoran]

>And why not replace the loop by the traversal of a binary
>tree ?
  The initialization loop?  It's presumably executed relatively
rarely so replace it with whatever you like, though I can't imagine
anything much simpler than a loop.
> (I just don't like that lookup table...)
  I love lookup tables.  With the price of RAM these days the dramatic
speed advantages of a lookup table are usually worth the memory cost.

Discrete random with given distribution ?

[Reinert Korsnes]

Hi,

Is it under Ada any "natural" way to generate random
elements of enumeration type (for example (a,b,c))
and with a given non-uniform probability distribution ?

I.e. so for example "a" is produced by probability 1/2
and "b" and "c" both with probability 1/4 ?

(Yes, I can do it via Float_Random, but the code looks ugly). 

reinert

Re: Discrete random with given distribution ?

[Larry Kilgallen]

In article <aea288$cen$1@dolly.uninett.no>, Reinert Korsnes <reinert.korsnes@chello.no> writes:
> Hi,
> 
> Is it under Ada any "natural" way to generate random
> elements of enumeration type (for example (a,b,c))
> and with a given non-uniform probability distribution ?

If you want truly random data, you need a connection to some physical
source of randomness, such as thermistor noise.

If you really meant "pseudo-random", continue the discussion.

> I.e. so for example "a" is produced by probability 1/2
> and "b" and "c" both with probability 1/4 ?
> 
> (Yes, I can do it via Float_Random, but the code looks ugly). 

Re: Discrete random with given distribution ?

[Reinert Korsnes]

Larry Kilgallen wrote:

> In article <aea288$cen$1@dolly.uninett.no>, Reinert Korsnes
> <reinert.korsnes@chello.no> writes:
>> Hi,
>> 
>> Is it under Ada any "natural" way to generate random
>> elements of enumeration type (for example (a,b,c))
>> and with a given non-uniform probability distribution ?
> 
> If you want truly random data, you need a connection to some physical
> source of randomness, such as thermistor noise.
> 
> If you really meant "pseudo-random", continue the discussion.

Pseudo.  The point this time was to make pretty code - not to be
subtle about random numbers.

reinert

> 
>> I.e. so for example "a" is produced by probability 1/2
>> and "b" and "c" both with probability 1/4 ?
>> 
>> (Yes, I can do it via Float_Random, but the code looks ugly).

Re: Discrete random with given distribution ?

[Wes Groleau]

> Is it under Ada any "natural" way to generate random
> elements of enumeration type (for example (a,b,c))
> and with a given non-uniform probability distribution ?
> 
> I.e. so for example "a" is produced by probability 1/2
> and "b" and "c" both with probability 1/4 ?

I don't know what's "natural" but decide for yourself
whether either of these are ugly:

  type Rand_Val is (A1, A2, B, C);

  Lookup : array (Rand_Val) of Out_Type :=
                 (A1 => A,
                  A2 => A,
                  B  => B,
                  C  => C);

  package Rand is
     new Ada.Numerics.Discrete_Random (Rand_Val);

   ....

   return Lookup (Rand.Random (Gen));

-------------

   -- FR has been created as random float)

   Threshold : array (Out_Type) of Float :=
     (A => 0.50,
      B => 0.75,   -- previous plus Probability (B)
      C => 1.00);  -- previous plus Probability (C)

   ....

   if Threshold (Threshold'Last) /= 1.0 then
      raise Programmer_Made_A_Mistake;
   end if;

   F := FR.Random;
   Loop_With_Jump_Out:
     for I in Out_Type loop
        if F > Threshold (I) then
           return I;
        end if;
     end loop Loop_With_Jump_Out;

-----

The first is awkward if the probabilities
are complicated.  The second could be a pain
doing all the addition, but one could easily
write a quarter-page program to prompt for
the probabilities and output the aggregate text.

If the list of items is long, and you want good
average speed, use an intermediate type arranged
biggest probabilities first, and

   return Lookup (I);

-- 
Wes Groleau
http://freepages.rootsweb.com/~wgroleau

Re: Discrete random with given distribution ?

[Kevin Rigotti]

Wes Groleau <wesgroleau@despammed.com> wrote in message
news:3D08AB63.6AF60F95@despammed.com...
>
>    F := FR.Random;
>    Loop_With_Jump_Out:
>      for I in Out_Type loop
>         if F > Threshold (I) then
>            return I;
>         end if;
>      end loop Loop_With_Jump_Out;

If any of the enumeration values has probability zero then you want to start
looking at the other end of the table and work backwards, until you find a
cumulative probability value that is greater than its predecessor, ie.

    for I in reverse Out_Type range 2 .. Out_Type'Last ....

Kevin
--
ATC Systems Group, QinetiQ, St Andrews Road, Malvern, Worcestershire WR14
3PS, UK
Phone +44 (0)1684 89 69 11, fax +44(0)1684 89 41 09
QinetiQ disclaimers and restrictions apply, details on request

Re: Discrete random with given distribution ?

[Jeffrey Carter]

Wes Groleau wrote:
> 
>    Threshold : array (Out_Type) of Float :=
>      (A => 0.50,
>       B => 0.75,   -- previous plus Probability (B)
>       C => 1.00);  -- previous plus Probability (C)
> 
>    F := FR.Random;
>    Loop_With_Jump_Out:
>      for I in Out_Type loop
>         if F > Threshold (I) then
>            return I;
>         end if;
>      end loop Loop_With_Jump_Out;

This can never return C. F <= 1.0, and for many generators, F < 1.0.
Also, if F = 0.9, this returns A.

I suspect you need

F < Threshold (I)

-- 
Jeff Carter
"Oh Lord, bless this thy hand grenade, that with it thou
mayst blow thine enemies to tiny bits, in thy mercy."
Monty Python and the Holy Grail

Re: Discrete random with given distribution ?

[Wes Groleau]

> This can never return C. F <= 1.0, and for many generators, F < 1.0.
> Also, if F = 0.9, this returns A.
> 
> I suspect you need
> 
> F < Threshold (I)

I suspect you're right.  Coding error.
The (intended) method is valid.

-- 
Wes Groleau
http://freepages.rootsweb.com/~wgroleau

Re: Discrete random with given distribution ?

[Gautier]

> >    F := FR.Random;

Jeffrey Carter:

> This can never return C. F <= 1.0, and for many generators, F < 1.0.
> Also, if F = 0.9, this returns A.
> 
> I suspect you need
> 
> F < Threshold (I)

I'm afraid this is not the main problem, since the probability of
F obtaining any value (-12.0, 0.6 or 1.0) is 0...

G.

Re: Discrete random with given distribution ?

[tmoran]

>  type Rand_Val is (A1, A2, B, C);
>
>  Lookup : array (Rand_Val) of Out_Type :=
>                 (A1 => A,
>                  A2 => A,
>                  B  => B,
>                  C  => C);
>
>  package Rand is
>     new Ada.Numerics.Discrete_Random (Rand_Val);
>
>   ....
>
>   return Lookup (Rand.Random (Gen));
>The first is awkward if the probabilities are complicated.
  If your probabilites are expressed with three digit accuracy, for
instance, just make a 1000 element Lookup table with the right number of
A's, B's, and C's at initialization time.  Few probabilities are
more precise than that.  (If they appear to be, .3333... for instance,
create Lookup with an appropriate LCM size other than a power of ten.)

Re: Discrete random with given distribution ?

[Reinert Korsnes]

tmoran@acm.org wrote:

>>  type Rand_Val is (A1, A2, B, C);
>>
>>  Lookup : array (Rand_Val) of Out_Type :=
>>                 (A1 => A,
>>                  A2 => A,
>>                  B  => B,
>>                  C  => C);
>>
>>  package Rand is
>>     new Ada.Numerics.Discrete_Random (Rand_Val);
>>
>>   ....
>>
>>   return Lookup (Rand.Random (Gen));
>>The first is awkward if the probabilities are complicated.
>   If your probabilites are expressed with three digit accuracy, for
> instance, just make a 1000 element Lookup table with the right number of
> A's, B's, and C's at initialization time.  Few probabilities are
> more precise than that.  (If they appear to be, .3333... for instance,
> create Lookup with an appropriate LCM size other than a power of ten.)

Yep, the quickest method may be to "discretize" the probabilities.
And the elements and their probabilities must in my case be given 
during runtime so the approach above must be somehow generalised.

reinert

Re: Discrete random with given distribution ?

[Wes Groleau]

> Yep, the quickest method may be to "discretize" the probabilities.
> And the elements and their probabilities must in my case be given
> during runtime so the approach above must be somehow generalised.

My first method must be determined before compile time
(Actually, there is a way to approximate it at run-time,
but it is a bit clumsy)

If the probabilities are to be determined at runtime,
my second method will work, if you

1. Fix the typo two other folks detected,

2. Write code to assign the the probabilities:

   Sum := 0.0

   for I in Threshold'Range loop

      Sum := Sum + Probability (I);
      Threshold (I) := Sum;

   end loop;

   if abs (Sum - 1.0) > Tolerance then
      raise Something_Is_Screwed_Up_Here;
   end if;

-- 
Wes Groleau
http://freepages.rootsweb.com/~wgroleau

Re: Discrete random with given distribution ?

[tmoran]

>My first method must be determined before compile time
>(Actually, there is a way to approximate it at run-time,
>but it is a bit clumsy)

>If the probabilities are to be determined at runtime,
>my second method will work, if you
  Why not just use your second method (summing probabilities till you
hit a certain level) to calculate, at run time, the population of
a lookup table, then use the first method to speedily index into
it to pick a random element.

Re: Discrete random with given distribution ?

[Wes Groleau]

> >My first method must be determined before compile time
> >(Actually, there is a way to approximate it at run-time,
> >but it is a bit clumsy)
> 
> >If the probabilities are to be determined at runtime,
> >my second method will work, if you
>   Why not just use your second method (summing probabilities till you
> hit a certain level) to calculate, at run time, the population of
> a lookup table, then use the first method to speedily index into
> it to pick a random element.

That's what I meant by a bit clumsy.  And it could be a fairly
large lookup table, depending on how accurate you try to be.

Size := 10 ** Digits_Of_Accuracy;

Lookup : array ( 1 .. Size ) of Discrete_Things;

Index := Lookup'First;
for I in Discrete_Things loop
   Span := Integer (Probability (I));
   for J in Index .. Index + Span - 1 loop
      Lookup (J) := I;
   end loop;
   Index := Index + Span;
end loop;

Hmmm, that's not as clumsy as I thought it would be.  :-)

-- 
Wes Groleau
http://freepages.rootsweb.com/~wgroleau

Re: Discrete random with given distribution ?

[tmoran]

> Hmmm, that's not as clumsy as I thought it would be.  :-)
  Unfortunately, roundoff complicates it a bit.   Consider:
Size = 10
type Discrete_Things is (Apples, Bananas, Cantaloupes, Dates, Eggs);
Probability := (Apples .. Dates => => 0.04, Eggs => 0.84);
Integer( Probability(I) * Size) = 0 for I in Apples .. Dates
  which is fine for a one digit accuracy, but
Integer( Probability(I) * Size) = 8 for I = Eggs
  leaving Lookup(9 .. 10) undefined.  And if you just select
from Lookup(1 .. 8), then Probability of random selection = Eggs is 100%

Re: Discrete random with given distribution ?

[Wes Groleau]

> > Hmmm, that's not as clumsy as I thought it would be.  :-)
>   Unfortunately, roundoff complicates it a bit.   Consider:
> Size = 10

If size is ten, you're only asking for one digit of accuracy.
You could compute three and round back down to two instead of
one.  Or you could go with method two.

> type Discrete_Things is (Apples, Bananas, Cantaloupes, Dates, Eggs);
> Probability := (Apples .. Dates => => 0.04, Eggs => 0.84);
> Integer( Probability(I) * Size) = 0 for I in Apples .. Dates
>   which is fine for a one digit accuracy, but
> Integer( Probability(I) * Size) = 8 for I = Eggs
>   leaving Lookup(9 .. 10) undefined.  And if you just select
> from Lookup(1 .. 8), then Probability of random selection = Eggs is 100%

-- 
Wes Groleau
http://freepages.rootsweb.com/~wgroleau

Re: Discrete random with given distribution ?

[tmoran]

> >   Unfortunately, roundoff complicates it a bit.   Consider:
> > Size = 10
>
> If size is ten, you're only asking for one digit of accuracy.
  Size = 10 was just a small example.  Consider:
Size : constant := 1000;
type Discrete_Things is range 1 .. 50;
Probability : constant array(Discrete_Things) of Float
  := (1 .. 49 => 1.0/1024.0, 50 => 1.0-49.0/1024.0);
the first 49, 1000*(1.0/1024.0), will each round to 1
the 50th, 1000*(1.0-49.0/1024.0), rounds to 952.
So, ignoring the roundoff problem, the simple algorithm will try to put
49+952= 1001 entries in Lookup, which only has room for 1000 entries.
While a programmer would, of course, recognize the problem when he wrote
"1.0/1024.0" and knew he was dealing with 3 digit accuracy, the problem
could easily creep in if the Probability array was calculated at run time.

Re: Discrete random with given distribution ?

[Wes Groleau]

> While a programmer would, of course, recognize the problem when he wrote
> "1.0/1024.0" and knew he was dealing with 3 digit accuracy, the problem
> could easily creep in if the Probability array was calculated at run time.

OK, well, perhaps a little extra code would solve it.
I did originally say it was clumsy.  I was surprised
at how simple it _seemed_ to be when I actually
wrote it out.

-- 
Wes Groleau
http://freepages.rootsweb.com/~wgroleau

Re: Discrete random with given distribution ?

[tmoran]

type ps is delta 0.001 range 0.0 .. 1000.0;
for ps'small use 0.001;
> > While a programmer would, of course, recognize the problem when he wrote
> > "1.0/1024.0" and knew he was dealing with 3 digit accuracy, the problem
> > could easily creep in if the Probability array was calculated at run time.
> ... at how simple it _seemed_ to be when I actually wrote it out.
  I think all you need is
  type Probability_Type is delta 0.001 range 0.0 .. 1000.0;
  for Probability_Type'small use 0.001;
  ...
  -- add all the Probability(I) and then divide them by their sum
  -- to guarantee they all add to 1.0 even after rounding to 3 digits.
In this example, of course, 1.0/1024.0 => 0.0 and the last Discrete_Type
is chosen 100% of the time, which may or may not be desirable behavior.

Re: Discrete random with given distribution ?

[tmoran]

Oops, ignore the two extraneous lines:
>type ps is delta 0.001 range 0.0 .. 1000.0;
>for ps'small use 0.001;

Re: Discrete random with given distribution ?

[Robert I. Eachus]

tmoran@acm.org wrote:

> type ps is delta 0.001 range 0.0 .. 1000.0;
> for ps'small use 0.001;


Actually better is to use Ada 95 decimal types:

type ps is delta 0.001 digits 3;

It is one of the things they are provided for...

Re: Discrete random with given distribution ?

[Gautier]

Reinert Korsnes:

> Is it under Ada any "natural" way to generate random
> elements of enumeration type (for example (a,b,c))
> and with a given non-uniform probability distribution ?

Maybe you'll be interested by gnatchopping the following...
____________________________________________________________
Gautier  --  http://www.mysunrise.ch/users/gdm/index.htm#Ada

NB: For a direct answer, address on the Web site!

---8<--gnatchop---8<--gnatchop---8<--gnatchop---8<--gnatchop---8<--

--From: Reinert Korsnes (reinert.korsnes@chello.no)
--Subject: Discrete random with given distribution ? 
--Newsgroups: comp.lang.ada
--View this article only 
--Date: 2002-06-13 05:25:17 PST 

--Hi,

--Is it under Ada any "natural" way to generate random
--elements of enumeration type (for example (a,b,c))
--and with a given non-uniform probability distribution ?

--I.e. so for example "a" is produced by probability 1/2
--and "b" and "c" both with probability 1/4 ?

generic
  type A_float is digits <>;
  type Thing is (<>);
  type Proba_array is array(Thing) of A_float;
  with function Uniform_random return A_float;
function Finite_distributed_random(proba: Proba_array) return Thing;

function Finite_distributed_random(proba: Proba_array) return Thing is
  subtype AF01 is A_float range 0.0 .. 1.0;
  U: constant AF01:= Uniform_random;
  p: AF01:= 0.0;
  Probabilities_dont_sum_up_to_1: exception;
begin
  for th in Thing loop
    p:= p + AF01(proba(th));
    if U < p then
      return th;
    end if;
  end loop;
  if (p-1.0) > A_float'epsilon then
    raise Probabilities_dont_sum_up_to_1;
  end if;
  return Thing'last;
end Finite_distributed_random;

with Finite_distributed_random;

with Ada.Text_IO;       use Ada.Text_IO;
with Ada.Float_Text_IO; use Ada.Float_Text_IO;
with Ada.Numerics.Float_Random;

procedure Test_Finite_distributed_random is

  G: Ada.Numerics.Float_Random.Generator;

  function URand return Float is
  begin
    return Ada.Numerics.Float_Random.Random(G);
  end;

  type ABC is (a,b,c);
  
  type Proba_array is array(ABC) of Float;

  function ABC_Rand is new
    Finite_distributed_random(
      A_float        => Float,
      Thing          => ABC,
      Proba_array    => Proba_array,
      Uniform_random => URand);
  
  procedure Test( tirages: Positive; proba: Proba_array ) is
    sample: array(ABC) of Natural;
    res: ABC;
  begin
    Ada.Numerics.Float_Random.Reset(G);
    sample:= (others => 0);
    for i in 1..tirages loop
      res:= ABC_Rand( proba );
      sample(res):= sample(res) + 1;
    end loop;
    for x in ABC loop
      Put( ABC'image(x) & " :    prob = ");
      Put( proba(x), 1,4,0);
      Put(                " ;    stat = ");
      Put( Float( sample(x) ) / Float( tirages ), 1,4,0);
      New_Line;
    end loop;
    New_Line;
  end Test;
  
begin
  Test( 100_000, (A=> 0.50, B=> 0.25, C=> 0.25) );
  Test( 100_000_000, (A=> 0.123, B=> 0.456, C=> 0.421) );
end;