{{task|Probability and statistics}} Given a mapping between items and their required probability of occurrence, generate a million items ''randomly'' subject to the given probabilities and compare the target probability of occurrence versus the generated values.

The total of all the probabilities should equal one. (Because floating point arithmetic is involved, this is subject to rounding errors).

Use the following mapping to test your programs:

aleph   1/5.0
beth    1/6.0
gimel   1/7.0
daleth  1/8.0
he      1/9.0
waw     1/10.0
zayin   1/11.0
heth    1759/27720 # adjusted so that probabilities add to 1

;Related task:

• [[Random number generator (device)]]

Ada

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

procedure Random_Distribution is
   Trials : constant := 1_000_000;
   type Outcome is (Aleph, Beth, Gimel, Daleth, He, Waw, Zayin, Heth);
   Pr : constant array (Outcome) of Uniformly_Distributed :=
        (1.0/5.0, 1.0/6.0, 1.0/7.0, 1.0/8.0, 1.0/9.0, 1.0/10.0, 1.0/11.0, 1.0);
   Samples : array (Outcome) of Natural := (others => 0);
   Value   : Uniformly_Distributed;
   Dice    : Generator;
begin
   for Try in 1..Trials loop
      Value := Random (Dice);
      for I in Pr'Range loop
         if Value <= Pr (I) then
            Samples (I) := Samples (I) + 1;
            exit;
         else
            Value := Value - Pr (I);
         end if;
      end loop;
   end loop;
      -- Printing the results
   for I in Pr'Range loop
      Put (Outcome'Image (I) & Character'Val (9));
      Put (Float'Image (Float (Samples (I)) / Float (Trials)) & Character'Val (9));
      if I = Heth then
         Put_Line (" rest");
      else
         Put_Line (Uniformly_Distributed'Image (Pr (I)));
      end if;
   end loop;
end Random_Distribution;

Sample output:

ALEPH    2.00167E-01     2.00000E-01
BETH     1.67212E-01     1.66667E-01
GIMEL    1.42290E-01     1.42857E-01
DALETH   1.24186E-01     1.25000E-01
HE       1.11455E-01     1.11111E-01
WAW      1.00325E-01     1.00000E-01
ZAYIN    9.10220E-02     9.09091E-02
HETH     6.33430E-02     rest

ALGOL 68

{{trans|C}}

{{works with|ALGOL 68|Revision 1 - no extensions to language used}}

{{works with|ALGOL 68G|Any - tested with release [http://sourceforge.net/projects/algol68/files/algol68g/algol68g-1.18.0/algol68g-1.18.0-9h.tiny.el5.centos.fc11.i386.rpm/download 1.18.0-9h.tiny]}} {{wont work with|ELLA ALGOL 68|Any (with appropriate job cards) - tested with release [http://sourceforge.net/projects/algol68/files/algol68toc/algol68toc-1.8.8d/algol68toc-1.8-8d.fc9.i386.rpm/download 1.8-8d] - due to extensive use of FORMATted transput}}

INT trials = 1 000 000;

MODE LREAL = LONG REAL;

MODE ITEM = STRUCT(
  STRING name,
  INT prob count,
  LREAL expect,
        mapping
);
INT col width = 9;
FORMAT real repr = $g(-col width+1, 6)$,
       item repr = $"Name: "g", Prob count: "g(0)", Expect: "f(real repr)", Mapping: ", f(real repr)l$;

[8]ITEM items := (
  ( "aleph",  0, ~, ~ ),
  ( "beth",   0, ~, ~ ),
  ( "gimel",  0, ~, ~ ),
  ( "daleth", 0, ~, ~ ),
  ( "he",     0, ~, ~ ),
  ( "waw",    0, ~, ~ ),
  ( "zayin",  0, ~, ~ ),
  ( "heth",   0, ~, ~ )
);

main:
(
  LREAL offset = 5; # const #

# initialise items #
  LREAL total sum := 0;
  FOR i FROM LWB items TO UPB items - 1 DO
    expect OF items[i] := 1/(i-1+offset);
    total sum +:= expect OF items[i]
  OD;
  expect OF items[UPB items] := 1 - total sum;

  mapping OF items[LWB items] := expect OF items[LWB items];
  FOR i FROM LWB items + 1 TO UPB items DO
    mapping OF items[i] := mapping OF items[i-1] + expect OF items[i]
  OD;

  # printf((item repr, items)) #

# perform the sampling #
  PROC sample = (REF[]LREAL mapping)INT:(
    INT out;
    LREAL rand real = random;
    FOR j FROM LWB items TO UPB items DO
      IF rand real < mapping[j] THEN
        out := j;
	done
      FI
    OD;
    done: out
  );

  FOR i TO trials DO
      prob count OF items[sample(mapping OF items)] +:= 1
  OD;

  FORMAT indent = $17k$;

# print the results #
  printf(($"Trials: "g(0)l$, trials));
  printf(($"Items:"$,indent));
  FOR i FROM LWB items TO UPB items DO printf(($gn(col width)k" "$, name OF items[i])) OD;
  printf(($l"Target prob.:"$, indent, $f(real repr)" "$, expect OF items));
  printf(($l"Attained prob.:"$, indent));
  FOR i FROM LWB items TO UPB items DO printf(($f(real repr)" "$, prob count OF items[i]/trials)) OD;
  printf($l$)
)

Sample output:

Trials: 1000000
Items:          aleph    beth     gimel    daleth   he       waw      zayin    heth
Target prob.:   0.200000 0.166667 0.142857 0.125000 0.111111 0.100000 0.090909 0.063456
Attained prob.: 0.199987 0.166917 0.142531 0.124203 0.111338 0.099702 0.091660 0.063662

AutoHotkey

contributed by Laszlo on the ahk [http://www.autohotkey.com/forum/post-276279.html#276279 forum]

s1 := "aleph",   p1 := 1/5.0                       ; Input
s2 := "beth",    p2 := 1/6.0
s3 := "gimel",   p3 := 1/7.0
s4 := "daleth",  p4 := 1/8.0
s5 := "he",      p5 := 1/9.0
s6 := "waw",     p6 := 1/10.0
s7 := "zayin",   p7 := 1/11.0
s8 := "heth",    p8 := 1-p1-p2-p3-p4-p5-p6-p7
n := 8, r0 := 0, r%n% := 1                         ; auxiliary data

Loop % n-1
   i := A_Index-1, r%A_Index% := r%i% + p%A_Index% ; cummulative distribution

Loop 1000000 {
   Random R, 0, 1.0
   Loop %n%                                        ; linear search
      If (R < r%A_Index%) {
          c%A_Index%++
          Break
      }
}
                                                   ; Output
Loop %n%
   t .= s%A_Index% "`t" p%A_Index% "`t" c%A_Index%*1.0e-6 "`n"
Msgbox %t%

/*
output:
---------------------------
aleph  0.200000   0.199960
beth   0.166667   0.166146
gimel  0.142857   0.142624
daleth 0.125000   0.124924
he     0.111111   0.111226
waw    0.100000   0.100434
zayin  0.090909   0.091344
heth   0.063456   0.063342
---------------------------
*/

AWK

#!/usr/bin/awk -f

BEGIN {
    ITERATIONS = 1000000
    delete symbMap
    delete probMap
    delete counts
    initData();

    for (i = 0; i < ITERATIONS; i++) {
        distribute(rand())
    }
    showDistributions()

    exit
}

function distribute(rnd,    cnt, symNum, sym, symPrb) {
    cnt = length(symbMap)
    for (symNum = 1; symNum <= cnt; symNum++) {
        sym = symbMap[symNum];
        symPrb = probMap[sym];
        rnd -= symPrb;
        if (rnd <= 0) {
            counts[sym]++
            return;
        }
    }
}

function showDistributions(    s, sym, prb, actSum, expSum, totItr) {
    actSum = 0.0
    expSum = 0.0
    totItr = 0
    printf "%-7s  %-7s  %-5s  %-5s\n", "symb", "num.", "act.", "expt."
    print  "-------  -------  -----  -----"
    for (s = 1; s <= length(symbMap); s++) {
        sym = symbMap[s]
        prb = counts[sym]/ITERATIONS
        actSum += prb
        expSum += probMap[sym]
        totItr += counts[sym]
        printf "%-7s  %7d  %1.3f  %1.3f\n", sym, counts[sym], prb, probMap[sym]
    }
    print  "-------  -------  -----  -----"
    printf "Totals:  %7d  %1.3f  %1.3f\n", totItr, actSum, expSum
}

function initData(    sym) {
    srand()

    probMap["aleph"]  = 1.0 / 5.0
    probMap["beth"]   = 1.0 / 6.0
    probMap["gimel"]  = 1.0 / 7.0
    probMap["daleth"] = 1.0 / 8.0
    probMap["he"]     = 1.0 / 9.0
    probMap["waw"]    = 1.0 / 10.0
    probMap["zyin"]   = 1.0 / 11.0
    probMap["heth"]   = 1759.0 / 27720.0

    symbMap[1] = "aleph"
    symbMap[2] = "beth"
    symbMap[3] = "gimel"
    symbMap[4] = "daleth"
    symbMap[5] = "he"
    symbMap[6] = "waw"
    symbMap[7] = "zyin"
    symbMap[8] = "heth"

    for (sym in probMap)
        counts[sym] = 0;
}

Example output:

symb num. act. expt.

aleph 200598 0.201 0.200 beth 166317 0.166 0.167 gimel 142391 0.142 0.143 daleth 125051 0.125 0.125 he 110658 0.111 0.111 waw 100464 0.100 0.100 zyin 90649 0.091 0.091 heth 63872 0.064 0.063

Totals: 1000000 1.000 1.000

Rounding off makes the results look perfect.

BBC BASIC

      DIM item$(7), prob(7), cnt%(7)
      item$() = "aleph","beth","gimel","daleth","he","waw","zayin","heth"
      prob()  = 1/5.0, 1/6.0, 1/7.0, 1/8.0, 1/9.0, 1/10.0, 1/11.0, 1759/27720
      IF ABS(SUM(prob())-1) > 1E-6 ERROR 100, "Probabilities don't sum to 1"

      FOR trial% = 1 TO 1E6
        r = RND(1)
        p = 0
        FOR i% = 0 TO DIM(prob(),1)
          p += prob(i%)
          IF r < p cnt%(i%) += 1 : EXIT FOR
        NEXT
      NEXT

      @% = &2060A
      PRINT "Item        actual    theoretical"
      FOR i% = 0 TO DIM(item$(),1)
        PRINT item$(i%), cnt%(i%)/1E6, prob(i%)
      NEXT

'''Output:'''

Item        actual    theoretical
aleph       0.200306  0.200000
beth        0.165963  0.166667
gimel       0.143089  0.142857
daleth      0.125387  0.125000
he          0.111057  0.111111
waw         0.100098  0.100000
zayin       0.091031  0.090909
heth        0.063069  0.063456

C

#include <stdio.h>
#include <stdlib.h>

/* pick a random index from 0 to n-1, according to probablities listed
   in p[] which is assumed to have a sum of 1. The values in the probablity
   list matters up to the point where the sum goes over 1 */
int rand_idx(double *p, int n)
{
	double s = rand() / (RAND_MAX + 1.0);
	int i;
	for (i = 0; i < n - 1 && (s -= p[i]) >= 0; i++);
	return i;
}

#define LEN 8
#define N 1000000
int main()
{
	const char *names[LEN] = { "aleph", "beth", "gimel", "daleth",
			  "he", "waw", "zayin", "heth" };
	double s, p[LEN] = { 1./5, 1./6, 1./7, 1./8, 1./9, 1./10, 1./11, 1e300 };
	int i, count[LEN] = {0};

	for (i = 0; i < N; i++) count[rand_idx(p, LEN)] ++;

	printf("  Name  Count    Ratio Expected\n");
	for (i = 0, s = 1; i < LEN; s -= p[i++])
		printf("%6s%7d %7.4f%% %7.4f%%\n",
			names[i], count[i],
			(double)count[i] / N * 100,
			((i < LEN - 1) ? p[i] : s) * 100);

	return 0;
}

output Name Count Ratio Expected aleph 199928 19.9928% 20.0000% beth 166489 16.6489% 16.6667% gimel 143211 14.3211% 14.2857% daleth 125257 12.5257% 12.5000% he 110849 11.0849% 11.1111% waw 99935 9.9935% 10.0000% zayin 91001 9.1001% 9.0909% heth 63330 6.3330% 6.3456%

## C++


```cpp
#include <cstdlib>
#include <iostream>
#include <vector>
#include <utility>
#include <algorithm>
#include <ctime>
#include <iomanip>

int main( ) {
   typedef std::vector<std::pair<std::string, double> >::const_iterator SPI ;
   typedef std::vector<std::pair<std::string , double> > ProbType ;
   ProbType probabilities ;
   probabilities.push_back( std::make_pair( "aleph" , 1/5.0 ) ) ;
   probabilities.push_back( std::make_pair( "beth" , 1/6.0 ) ) ;
   probabilities.push_back( std::make_pair( "gimel" , 1/7.0 ) ) ;
   probabilities.push_back( std::make_pair( "daleth" , 1/8.0 ) ) ;
   probabilities.push_back( std::make_pair( "he" , 1/9.0 ) ) ;
   probabilities.push_back( std::make_pair( "waw" , 1/10.0 ) ) ;
   probabilities.push_back( std::make_pair( "zayin" , 1/11.0 ) ) ;
   probabilities.push_back( std::make_pair( "heth" , 1759/27720.0 ) ) ;
   std::vector<std::string> generated ; //for the strings that are generatod
   std::vector<int> decider ; //holds the numbers that determine the choice of letters
   for ( int i = 0 ; i < probabilities.size( ) ; i++ ) {
      if ( i == 0 ) {
	 decider.push_back( 27720 * (probabilities[ i ].second) ) ;
      }
      else {
	 int number = 0 ;
	 for ( int j = 0 ; j < i ; j++ ) {
	    number +=  27720 * ( probabilities[ j ].second ) ;
	 }
	 number += 27720 * probabilities[ i ].second ;
	 decider.push_back( number ) ;
      }
   }
   srand( time( 0 ) ) ;
   for ( int i = 0 ; i < 1000000 ; i++ ) {
      int randnumber = rand( ) % 27721 ;
      int j = 0 ;
      while ( randnumber > decider[ j ] )
	 j++ ;
      generated.push_back( ( probabilities[ j ]).first ) ;
   }
   std::cout << "letter  frequency attained   frequency expected\n" ;
   for ( SPI i = probabilities.begin( ) ; i != probabilities.end( ) ; i++ ) {
      std::cout << std::left << std::setw( 8 ) << i->first ;
      int found = std::count ( generated.begin( ) , generated.end( ) , i->first ) ;
      std::cout << std::left << std::setw( 21 ) << found / 1000000.0 ;
      std::cout << std::left << std::setw( 17 ) << i->second << '\n' ;
   }
   return 0 ;
}

Output:

letter  frequency attained   frequency expected
aleph   0.200089             0.2
beth    0.16695              0.166667
gimel   0.142693             0.142857
daleth  0.124859             0.125
he      0.111258             0.111111
waw     0.099665             0.1
zayin   0.090654             0.0909091
heth    0.063832             0.063456

USE: rosettacode.proba 1000000 data example ``` outputs Key Value expected heth: 0.063469 0.063456 waw: 0.100226 0.100000 daleth: 0.125844 0.125000 beth: 0.166264 0.166667 zayin: 0.090806 0.090909 he: 0.110562 0.111111 aleph: 0.199868 0.200000 gimel: 0.142961 0.142857 ``` ## Forth ```forth include random.fs \ common factors of desired probabilities (1/5 .. 1/11) 2 2 * 2 * 3 * 3 * 5 * 7 * 11 * constant denom \ 27720 \ represent each probability as the numerator with 27720 as the denominator : ,numerators ( max min -- ) do denom i / , loop ; \ final item is 27720 - sum(probs) : ,remainder ( denom addr len -- ) cells bounds do i @ - 1 cells +loop , ; create probs 12 5 ,numerators denom probs 7 ,remainder create bins 8 cells allot : choose ( -- 0..7 ) denom random 8 0 do probs i cells + @ - dup 0< if drop i unloop exit then loop abort" can't get here" ; : trials ( n -- ) 0 do 1 bins choose cells + +! loop ; : str-table create ( c-str ... n -- ) 0 do , loop does> ( n -- str len ) swap cells + @ count ; here ," heth" here ," zayin" here ," waw" here ," he" here ," daleth" here ," gimel" here ," beth" here ," aleph" 8 str-table names : .header cr ." Name" #tab emit ." Prob" #tab emit ." Actual" #tab emit ." Error" ; : .result ( n -- ) cr dup names type #tab emit dup cells probs + @ s>f denom s>f f/ fdup f. #tab emit dup cells bins + @ s>f 1e6 f/ fdup f. #tab emit f- fabs fs. ; : .results .header 8 0 do i .result loop ; ``` ```txt bins 8 cells erase 3 set-precision 1000000 trials .results Name Prob Actual Error aleph 0.2 0.2 9.90E-5 beth 0.167 0.167 4.51E-4 gimel 0.143 0.142 4.99E-4 daleth 0.125 0.125 1.82E-4 he 0.111 0.111 2.10E-4 waw 0.1 0.1 3.30E-5 zayin 0.0909 0.0912 2.77E-4 heth 0.0635 0.0636 9.70E-5 ok ``` ## Fortran {{works with|Fortran|90 and later}} ```fortran PROGRAM PROBS IMPLICIT NONE INTEGER, PARAMETER :: trials = 1000000 INTEGER :: i, j, probcount(8) = 0 REAL :: expected(8), mapping(8), rnum CHARACTER(6) :: items(8) = (/ "aleph ", "beth ", "gimel ", "daleth", "he ", "waw ", "zayin ", "heth " /) expected(1:7) = (/ (1.0/i, i=5,11) /) expected(8) = 1.0 - SUM(expected(1:7)) mapping(1) = 1.0 / 5.0 DO i = 2, 7 mapping(i) = mapping(i-1) + 1.0/(i+4.0) END DO mapping(8) = 1.0 DO i = 1, trials CALL RANDOM_NUMBER(rnum) DO j = 1, 8 IF (rnum < mapping(j)) THEN probcount(j) = probcount(j) + 1 EXIT END IF END DO END DO WRITE(*, "(A,I10)") "Trials: ", trials WRITE(*, "(A,8A10)") "Items: ", items WRITE(*, "(A,8F10.6)") "Target Probability: ", expected WRITE(*, "(A,8F10.6)") "Attained Probability:", REAL(probcount) / REAL(trials) ENDPROGRAM PROBS ``` Sample Output: ```txt Trials: 1000000 Items: aleph beth gimel daleth he waw zayin heth Target Probability: 0.200000 0.166667 0.142857 0.125000 0.111111 0.100000 0.090909 0.063456 Attained Probability: 0.199631 0.166907 0.142488 0.124920 0.110906 0.099943 0.091775 0.063430 ``` ## FreeBASIC ```freebasic ' FB 1.05.0 Win64 Dim letters (0 To 7) As String = {"aleph", "beth", "gimel", "daleth", "he", "waw", "zayin", "heth"} Dim actual (0 To 7) As Integer '' all zero by default Dim probs (0 To 7) As Double = {1/5.0, 1/6.0, 1/7.0, 1/8.0, 1/9.0, 1/10.0, 1/11.0} Dim cumProbs (0 To 7) As Double cumProbs(0) = probs(0) For i As Integer = 1 To 6 cumProbs(i) = cumProbs(i - 1) + probs(i) Next cumProbs(7) = 1.0 probs(7) = 1.0 - cumProbs(6) Randomize Dim rand As Double Dim n As Double = 1000000 Dim sum As Double = 0.0 For i As Integer = 1 To n rand = Rnd '' random number where 0 <= rand < 1 Select case rand Case Is <= cumProbs(0) actual(0) += 1 Case Is <= cumProbs(1) actual(1) += 1 Case Is <= cumProbs(2) actual(2) += 1 Case Is <= cumProbs(3) actual(3) += 1 Case Is <= cumProbs(4) actual(4) += 1 Case Is <= cumProbs(5) actual(5) += 1 Case Is <= cumProbs(6) actual(6) += 1 Case Else actual(7) += 1 End Select Next Dim sumActual As Double = 0 Print "Letter", " Actual", "Expected" Print "------", "--------", "--------" For i As Integer = 0 To 7 Print letters(i), Print Using "#.######"; actual(i)/n; sumActual += actual(i)/n Print , Using "#.######"; probs(i) Next Print , "--------", "--------" Print , Using "#.######"; sumActual; Print , Using "#.######"; 1.000000 Print Print "Press any key to quit" Sleep ``` {{out}} ```txt Letter Actual Expected ------ -------- -------- aleph 0.199987 0.200000 beth 0.166663 0.166667 gimel 0.143134 0.142857 daleth 0.125132 0.125000 he 0.110772 0.111111 waw 0.100236 0.100000 zayin 0.090664 0.090909 heth 0.063412 0.063456 -------- -------- 1.000000 1.000000 ``` ## Go ```go package main import ( "fmt" "math/rand" "time" ) type mapping struct { item string pr float64 } func main() { // input mapping m := []mapping{ {"aleph", 1 / 5.}, {"beth", 1 / 6.}, {"gimel", 1 / 7.}, {"daleth", 1 / 8.}, {"he", 1 / 9.}, {"waw", 1 / 10.}, {"zayin", 1 / 11.}, {"heth", 1759 / 27720.}} // adjusted so that probabilities add to 1 // cumulative probability cpr := make([]float64, len(m)-1) var c float64 for i := 0; i < len(m)-1; i++ { c += m[i].pr cpr[i] = c } // generate const samples = 1e6 occ := make([]int, len(m)) rand.Seed(time.Now().UnixNano()) for i := 0; i < samples; i++ { r := rand.Float64() for j := 0; ; j++ { if r < cpr[j] { occ[j]++ break } if j == len(cpr)-1 { occ[len(cpr)]++ break } } } // report fmt.Println(" Item Target Generated") var totalTarget, totalGenerated float64 for i := 0; i < len(m); i++ { target := m[i].pr generated := float64(occ[i]) / samples fmt.Printf("%6s %8.6f %8.6f\n", m[i].item, target, generated) totalTarget += target totalGenerated += generated } fmt.Printf("Totals %8.6f %8.6f\n", totalTarget, totalGenerated) } ``` Output: ```txt Item Target Generated aleph 0.200000 0.199509 beth 0.166667 0.167194 gimel 0.142857 0.143293 daleth 0.125000 0.124869 he 0.111111 0.110896 waw 0.100000 0.099849 zayin 0.090909 0.090789 heth 0.063456 0.063601 Totals 1.000000 1.000000 ``` ## Haskell ```haskell import System.Random (newStdGen, randomRs) dataBinCounts :: [Float] -> [Float] -> [Int] dataBinCounts thresholds range = let sampleSize = length range xs = ((-) sampleSize . length . flip filter range . (<)) <$> thresholds in zipWith (-) (xs ++ [sampleSize]) (0 : xs) main :: IO () main = do g <- newStdGen let fractions = recip <$> [5 .. 11] :: [Float] expected = fractions ++ [1 - sum fractions] actual = ((/ 1000000.0) . fromIntegral) <$> dataBinCounts (scanl1 (+) expected) (take 1000000 (randomRs (0, 1) g)) piv n = take n . (++ repeat ' ') putStrLn " expected actual" mapM_ putStrLn $ zipWith3 (\l s c -> piv 7 l ++ piv 13 (show s) ++ piv 12 (show c)) ["aleph", "beth", "gimel", "daleth", "he", "waw", "zayin", "heth"] expected actual ``` {{Out}} Sample ```txt expected actual aleph 0.2 0.200597 beth 0.16666667 0.167192 gimel 0.14285715 0.142781 daleth 0.125 0.124556 he 0.11111111 0.111128 waw 0.1 9.9671e-2 zayin 9.090909e-2 9.0294e-2 heth 6.345594e-2 6.3781e-2 ``` ## HicEst ```HicEst REAL :: trials=1E6, n=8, map(n), limit(n), expected(n), outcome(n) expected = 1 / ($ + 4) expected(n) = 1 - SUM(expected) + expected(n) map = expected map = map($) + map($-1) DO i = 1, trials random = RAN(1) limit = random > map item = INDEX(limit, 0) outcome(item) = outcome(item) + 1 ENDDO outcome = outcome / trials DLG(Text=expected, Text=outcome, Y=0) ``` Exported from the spreadsheet-like DLG function: ```txt 0.2 0.199908 0.1666667 0.166169 0.1428571 0.142722 0.125 0.124929 0.1111111 0.111706 0.1 0.099863 0.0909091 0.090965 0.063456 0.063738 ``` =={{header|Icon}} and {{header|Unicon}}== ```Icon record Item(value, probability) procedure find_item (items, v) sum := 0.0 every item := !items do { if v < sum+item.probability then return item.value else sum +:= item.probability } fail # v exceeded 1.0 end # -- helper procedures # count the number of occurrences of i in list l, # assuming the items are strings procedure count (l, i) result := 0.0 every x := !l do if x == i then result +:= 1 return result end procedure rand_float () return ?1000/1000.0 end # -- test the procedure procedure main () items := [ Item("aleph", 1/5.0), Item("beth", 1/6.0), Item("gimel", 1/7.0), Item("daleth", 1/8.0), Item("he", 1/9.0), Item("waw", 1/10.0), Item("zayin", 1/11.0), Item("heth", 1759/27720.0) ] # collect a sample of results sample := [] every (1 to 1000000) do push (sample, find_item(items, rand_float ())) # return comparison of expected vs actual probability every item := !items do write (right(item.value, 7) || " " || left(item.probability, 15) || " " || left(count(sample, item.value)/*sample, 6)) end ``` Output: ```txt aleph 0.2 0.1988 beth 0.1666666667 0.1676 gimel 0.1428571429 0.1431 daleth 0.125 0.1249 he 0.1111111111 0.1112 waw 0.1 0.0996 zayin 0.09090909091 0.0908 heth 0.06345598846 0.0636 ``` ## J ```J main=: verb define hdr=. ' target actual ' lbls=. ; ,:&.> ;:'aleph beth gimel daleth he waw zayin heth' prtn=. +/\ pt=. (, 1-+/)1r1%5+i.7 da=. prtn I. ?y # 0 pa=. y%~ +/ da =/ i.8 hdr, lbls,. 9j6 ": |: pt,:pa ) Note 'named abbreviations' hdr (header) lbls (labels) pt (target proportions) prtn (partitions corresponding to target proportions) da (distribution of actual values among partitions) pa (actual proportions) ) ``` Example use: ```j main 1e6 target actual aleph 0.200000 0.200344 beth 0.166667 0.166733 gimel 0.142857 0.142611 daleth 0.125000 0.124458 he 0.111111 0.111455 waw 0.100000 0.099751 zayin 0.090909 0.091121 heth 0.063456 0.063527 ``` Note that there is no rounding error in summing the proportions, as they are represented as rational numbers, not floating-point approximations. ```J pt=. (, 1-+/)1r1%5+i.7 pt 1r5 1r6 1r7 1r8 1r9 1r10 1r11 1759r27720 +/pt 1 ``` ## Java {{trans|C}} ```java public class Prob{ static long TRIALS= 1000000; private static class Expv{ public String name; public int probcount; public double expect; public double mapping; public Expv(String name, int probcount, double expect, double mapping){ this.name= name; this.probcount= probcount; this.expect= expect; this.mapping= mapping; } } static Expv[] items= {new Expv("aleph", 0, 0.0, 0.0), new Expv("beth", 0, 0.0, 0.0), new Expv("gimel", 0, 0.0, 0.0), new Expv("daleth", 0, 0.0, 0.0), new Expv("he", 0, 0.0, 0.0), new Expv("waw", 0, 0.0, 0.0), new Expv("zayin", 0, 0.0, 0.0), new Expv("heth", 0, 0.0, 0.0)}; public static void main(String[] args){ int i, j; double rnum, tsum= 0.0; for(i= 0, rnum= 5.0;i < 7;i++, rnum+= 1.0){ items[i].expect= 1.0 / rnum; tsum+= items[i].expect; } items[7].expect= 1.0 - tsum; items[0].mapping= 1.0 / 5.0; for(i= 1;i < 7;i++){ items[i].mapping= items[i - 1].mapping + 1.0 / ((double)i + 5.0); } items[7].mapping= 1.0; for(i= 0;i < TRIALS;i++){ rnum= Math.random(); for(j= 0;j < 8;j++){ if(rnum < items[j].mapping){ items[j].probcount++; break; } } } System.out.printf("Trials: %d\n", TRIALS); System.out.printf("Items: "); for(i= 0;i < 8;i++) System.out.printf("%-8s ", items[i].name); System.out.printf("\nTarget prob.: "); for(i= 0;i < 8;i++) System.out.printf("%8.6f ", items[i].expect); System.out.printf("\nAttained prob.: "); for(i= 0;i < 8;i++) System.out.printf("%8.6f ", (double)(items[i].probcount) / (double)TRIALS); System.out.printf("\n"); } } ``` Output: ```txt Trials: 1000000 Items: aleph beth gimel daleth he waw zayin heth Target prob.: 0.200000 0.166667 0.142857 0.125000 0.111111 0.100000 0.090909 0.063456 Attained prob.: 0.199615 0.167517 0.142612 0.125211 0.110970 0.099614 0.091002 0.063459 ``` {{works with|Java|1.5+}} ```java5 import java.util.EnumMap; public class Prob { public static long TRIALS= 1000000; public enum Glyph{ ALEPH, BETH, GIMEL, DALETH, HE, WAW, ZAYIN, HETH; } public static EnumMap probs = new EnumMap(Glyph.class){{ put(Glyph.ALEPH, 1/5.0); put(Glyph.BETH, 1/6.0); put(Glyph.GIMEL, 1/7.0); put(Glyph.DALETH, 1/8.0); put(Glyph.HE, 1/9.0); put(Glyph.WAW, 1/10.0); put(Glyph.ZAYIN, 1/11.0); put(Glyph.HETH, 1759./27720); }}; public static EnumMap counts = new EnumMap(Glyph.class){{ put(Glyph.ALEPH, 0.);put(Glyph.BETH, 0.); put(Glyph.GIMEL, 0.);put(Glyph.DALETH, 0.); put(Glyph.HE, 0.);put(Glyph.WAW, 0.); put(Glyph.ZAYIN, 0.);put(Glyph.HETH, 0.); }}; public static void main(String[] args){ System.out.println("Target probabliities:\t" + probs); for(long i = 0; i < TRIALS; i++){ Glyph choice = getChoice(); counts.put(choice, counts.get(choice) + 1); } //correct the counts to probablities in (0..1] for(Glyph glyph:counts.keySet()){ counts.put(glyph, counts.get(glyph) / TRIALS); } System.out.println("Actual probabliities:\t" + counts); } private static Glyph getChoice() { double rand = Math.random(); for(Glyph item:Glyph.values()){ if(rand < probs.get(item)){ return item; } rand -= probs.get(item); } return null; } } ``` Output: ```txt Target probabliities: {ALEPH=0.2, BETH=0.16666666666666666, GIMEL=0.14285714285714285, DALETH=0.125, HE=0.1111111111111111, WAW=0.1, ZAYIN=0.09090909090909091, HETH=0.06345598845598846} Actual probabliities: {ALEPH=0.200794, BETH=0.165916, GIMEL=0.143286, DALETH=0.124727, HE=0.110818, WAW=0.100168, ZAYIN=0.090878, HETH=0.063413} ``` ## JavaScript ### ES5 Fortunately, iterating over properties added to an object maintains the insertion order. ```javascript var probabilities = { aleph: 1/5.0, beth: 1/6.0, gimel: 1/7.0, daleth: 1/8.0, he: 1/9.0, waw: 1/10.0, zayin: 1/11.0, heth: 1759/27720 }; var sum = 0; var iterations = 1000000; var cumulative = {}; var randomly = {}; for (var name in probabilities) { sum += probabilities[name]; cumulative[name] = sum; randomly[name] = 0; } for (var i = 0; i < iterations; i++) { var r = Math.random(); for (var name in cumulative) { if (r <= cumulative[name]) { randomly[name]++; break; } } } for (var name in probabilities) // using WSH WScript.Echo(name + "\t" + probabilities[name] + "\t" + randomly[name]/iterations); ``` output: ```txt aleph 0.2 0.200597 beth 0.16666666666666666 0.166527 gimel 0.14285714285714285 0.142646 daleth 0.125 0.124613 he 0.1111111111111111 0.111342 waw 0.1 0.099888 zayin 0.09090909090909091 0.091141 heth 0.06345598845598846 0.063246 ``` ### ES6 By functional composition: {{Trans|Haskell}} ```JavaScript (() => { 'use strict'; // GENERIC FUNCTIONS ----------------------------------------------------- // transpose :: [[a]] -> [[a]] const transpose = xs => xs[0].map((_, iCol) => xs.map(row => row[iCol])); // justifyLeft :: Int -> Char -> Text -> Text const justifyLeft = (n, cFiller, strText) => n > strText.length ? ( (strText + cFiller.repeat(n)) .substr(0, n) ) : strText; // 2 or more arguments // curry :: Function -> Function const curry = (f, ...args) => { const go = xs => xs.length >= f.length ? (f.apply(null, xs)) : function () { return go(xs.concat([].slice.apply(arguments))); }; return go([].slice.call(args, 1)); }; // zipWith :: (a -> b -> c) -> [a] -> [b] -> [c] const zipWith = (f, xs, ys) => { const ny = ys.length; return (xs.length <= ny ? xs : xs.slice(0, ny)) .map((x, i) => f(x, ys[i])); }; // subtract :: (Num a) => a -> a -> a const subtract = (x, y) => y - x; // scanl1 :: (a -> a -> a) -> [a] -> [a] const scanl1 = (f, xs) => xs.length > 0 ? scanl(f, xs[0], xs.slice(1)) : []; // scanl :: (b -> a -> b) -> b -> [a] -> [b] const scanl = (f, startValue, xs) => xs.reduce((a, x) => { const v = f(a.acc, x); return { acc: v, scan: a.scan.concat(v) }; }, { acc: startValue, scan: [startValue] }) .scan; // unwords :: [String] -> String const unwords = xs => xs.join(' '); // PROBABILISTIC CHOICE -------------------------------------------------- // samples :: Int -> Int -> [Float] const samples = n => Array.from({ length: n }, Math.random); // thresholds :: Float const thresholds = scanl1( (a, b) => a + b, [5, 6, 7, 8, 9, 10, 11].map(x => 1 / x) ) .concat(1); // expected :: Float -> Float const expected = limits => limits.map((x, i, xs) => i > 0 ? (x - xs[i - 1]) : x); // dataBinCounts :: [Float] -> [Float] -> [Int] const dataBinCounts = (thresholds, samples) => { const lng = samples.length, xs = thresholds .map(x => lng - samples.filter(v => v > x) .length); return zipWith(subtract, [0].concat(xs), xs.concat(lng)); }; // intSamples :: Integer const intSamples = 1000000; // aligned :: a -> String const aligned = x => justifyLeft(12, ' ', isNaN(x) ? x : x.toFixed(7)); return transpose([ ['', 'Aleph', 'Beit', 'Gimel', 'Dalet', 'He', 'Vav', 'Zayin', 'Chet'] .map(curry(justifyLeft)(7, ' ')), ['Expected'].concat(expected(thresholds)) .map(aligned), ['Observed'].concat(dataBinCounts(thresholds, samples(intSamples)) .map(x => x / intSamples)) .map(aligned) ]) .map(unwords) .join('\n'); })(); ``` {{Out}} Sample: ```txt Expected Observed Aleph 0.2000000 0.2002440 Beit 0.1666667 0.1665330 Gimel 0.1428571 0.1433880 Dalet 0.1250000 0.1244630 He 0.1111111 0.1112830 Vav 0.1000000 0.0998390 Zayin 0.0909091 0.0909630 Chet 0.0634560 0.0632870 ``` ## Julia I made the solution to this task more difficult than I had anticipated by using the Hebrew characters (rather than their anglicised names) as labels for the sampled collection of objects. In doing so, I encountered an interesting subtlety of bidirectional text in Unicode. Namely, that strong right-to-left characters, such as those of Hebrew, override the directionality of European digits, which have weak directionality. Because of this property of Unicode, my table of items and yields had its lines of data interpreted as if it were entirely Hebrew and output in reverse order (from my English speaking perspective). I was able to get the table to display as I liked on my terminal by preceding the the Hebrew characters by the Unicode RLI (right-to-left isolate) control character (U+2067). However, when I pasted this output into this Rosetta Code entry, the display reverted to the "backwards" version. Rather than continue the struggle, trying to force this entry to display as it does on my terminal, I created an alternative version of the table. This "Displayable Here" table adds "yields" to to each line, and this strong left-to-right text makes the whole line display as left-to-right (without the need for a RLI characer). ```Julia p = [1/i for i in 5:11] plen = length(p) q = [0.0, [sum(p[1:i]) for i = 1:plen]] plab = [char(i) for i in 0x05d0:(0x05d0+plen)] hi = 10^6 push!(p, 1.0 - sum(p)) plen += 1 accum = zeros(Int, plen) for i in 1:hi accum[sum(rand() .>= q)] += 1 end r = accum/hi println("Rates at which items are selected (", hi, " trials).") println(" Item Expected Actual") for i in 1:plen println(@sprintf(" \u2067%s %8.6f %8.6f", plab[i], p[i], r[i])) end println() println("Rates at which items are selected (", hi, " trials).") println(" Item Count Expected Actual") for i in 1:plen println(@sprintf(" %s yields %6d %8.6f %8.6f", plab[i], accum[i], p[i], r[i])) end ``` {{out}} '''Original''' This table displays properly on my terminal, but the lines of data are reversed in this display. ```txt Rates at which items are selected (1000000 trials). Item Expected Actual ⁧א 0.200000 0.199872 ⁧ב 0.166667 0.166618 ⁧ג 0.142857 0.143302 ⁧ד 0.125000 0.125040 ⁧ה 0.111111 0.110602 ⁧ו 0.100000 0.099833 ⁧ז 0.090909 0.091313 ⁧ח 0.063456 0.063420 ``` '''Displayable Here''' This is the same data, less elegantly presented but accurately displayed on both my terminal and here at Rosetta Code. ```txt Rates at which items are selected (1000000 trials). Item Count Expected Actual א yields 199872 0.200000 0.199872 ב yields 166618 0.166667 0.166618 ג yields 143302 0.142857 0.143302 ד yields 125040 0.125000 0.125040 ה yields 110602 0.111111 0.110602 ו yields 99833 0.100000 0.099833 ז yields 91313 0.090909 0.091313 ח yields 63420 0.063456 0.063420 ``` ## Kotlin {{trans|FreeBASIC}} ```scala // version 1.0.6 fun main(args: Array) { val letters = arrayOf("aleph", "beth", "gimel", "daleth", "he", "waw", "zayin", "heth") val actual = IntArray(8) val probs = doubleArrayOf(1/5.0, 1/6.0, 1/7.0, 1/8.0, 1/9.0, 1/10.0, 1/11.0, 0.0) val cumProbs = DoubleArray(8) cumProbs[0] = probs[0] for (i in 1..6) cumProbs[i] = cumProbs[i - 1] + probs[i] cumProbs[7] = 1.0 probs[7] = 1.0 - cumProbs[6] val n = 1000000 (1..n).forEach { val rand = Math.random() when { rand <= cumProbs[0] -> actual[0]++ rand <= cumProbs[1] -> actual[1]++ rand <= cumProbs[2] -> actual[2]++ rand <= cumProbs[3] -> actual[3]++ rand <= cumProbs[4] -> actual[4]++ rand <= cumProbs[5] -> actual[5]++ rand <= cumProbs[6] -> actual[6]++ else -> actual[7]++ } } var sumActual = 0.0 println("Letter\t Actual Expected") println("------\t-------- --------") for (i in 0..7) { val generated = actual[i].toDouble() / n println("${letters[i]}\t${String.format("%8.6f %8.6f", generated, probs[i])}") sumActual += generated } println("\t-------- --------") println("\t${"%8.6f".format(sumActual)} 1.000000") } ``` {{out}} ```txt Letter Actual Expected ------ -------- -------- aleph 0.199427 0.200000 beth 0.166862 0.166667 gimel 0.142756 0.142857 daleth 0.125442 0.125000 he 0.110868 0.111111 waw 0.100405 0.100000 zayin 0.090799 0.090909 heth 0.063441 0.063456 -------- -------- 1.000000 1.000000 ``` ## Liberty BASIC ```lb names$="aleph beth gimel daleth he waw zayin heth" dim sum(8) dim counter(8) s = 0 for i = 1 to 7 s = s+1/(i+4) sum(i)=s next N =1000000 ' number of throws for i =1 to N rand =rnd( 1) for j = 1 to 7 if sum(j)> rand then exit for next counter(j)=counter(j)+1 next print "Observed", "Intended" for i = 1 to 8 print word$(names$, i), using( "#.#####", counter(i) /N), using( "#.#####", 1/(i+4)) next ``` ## Lua ```lua items = {} items["aleph"] = 1/5.0 items["beth"] = 1/6.0 items["gimel"] = 1/7.0 items["daleth"] = 1/8.0 items["he"] = 1/9.0 items["waw"] = 1/10.0 items["zayin"] = 1/11.0 items["heth"] = 1759/27720 num_trials = 1000000 samples = {} for item, _ in pairs( items ) do samples[item] = 0 end math.randomseed( os.time() ) for i = 1, num_trials do z = math.random() for item, _ in pairs( items ) do if z < items[item] then samples[item] = samples[item] + 1 break; else z = z - items[item] end end end for item, _ in pairs( items ) do print( item, samples[item]/num_trials, items[item] ) end ``` Output ```txt gimel 0.142606 0.14285714285714 heth 0.063434 0.063455988455988 beth 0.166788 0.16666666666667 zayin 0.091097 0.090909090909091 daleth 0.124772 0.125 aleph 0.200541 0.2 he 0.1107 0.11111111111111 waw 0.100062 0.1 ``` ## Mathematica Built-in function can already do a weighted random choosing. Example for making a million random choices would be: ```Mathematica choices={{"aleph", 1/5},{"beth", 1/6},{"gimel", 1/7},{"daleth", 1/8},{"he", 1/9},{"waw", 1/10},{"zayin", 1/11},{"heth", 1759/27720}}; data=RandomChoice[choices[[All,2]]->choices[[All,1]],10^6]; ``` To compare the data we use the following code to make a table: ```Mathematica Grid[{#[[1]],N[Count[data,#[[1]]]/10^6],N[#[[2]]]}&/@choices] ``` gives back (item, attained prob., target prob.): ```txt aleph 0.200036 0.2 beth 0.166591 0.166667 gimel 0.142699 0.142857 daleth 0.125018 0.125 he 0.111306 0.111111 waw 0.100433 0.1 zayin 0.090671 0.0909091 heth 0.063246 0.063456 ``` ## MATLAB {{works with|MATLAB|with Statistics Toolbox}} ```MATLAB function probChoice choices = {'aleph' 'beth' 'gimel' 'daleth' 'he' 'waw' 'zayin' 'heth'}; w = [1/5 1/6 1/7 1/8 1/9 1/10 1/11 1759/27720]; R = randsample(length(w), 1e6, true, w); T = tabulate(R); fprintf('Value\tCount\tPercent\tGoal\n') for k = 1:size(T, 1) fprintf('%6s\t%.f\t%.2f%%\t%.2f%%\n', ... choices{k}, T(k, 2), T(k, 3), 100*w(k)) end end ``` {{out}} ```txt Value Count Percent Goal aleph 199635 19.96% 20.00% beth 166427 16.64% 16.67% gimel 143342 14.33% 14.29% daleth 125014 12.50% 12.50% he 111031 11.10% 11.11% waw 99920 9.99% 10.00% zayin 91460 9.15% 9.09% heth 63171 6.32% 6.35% ``` {{works with|MATLAB|without toolboxes}} ```MATLAB function probChoice choices = {'aleph' 'beth' 'gimel' 'daleth' 'he' 'waw' 'zayin' 'heth'}; w = [1/5 1/6 1/7 1/8 1/9 1/10 1/11 1759/27720]; nSamp = 1e6; nChoice = length(w); R = rand(nSamp, 1); wCS = cumsum(w); results = zeros(1, nChoice); fprintf('Value\tCount\tPercent\tGoal\n') for k = 1:nChoice choiceKIdxs = R < wCS(k); R(choiceKIdxs) = k; results(k) = sum(choiceKIdxs); fprintf('%6s\t%.f\t%.2f%%\t%.2f%%\n', ... choices{k}, results(k), 100*results(k)/nSamp, 100*w(k)) end end ``` {{out}} ```txt Value Count Percent Goal aleph 200327 20.03% 20.00% beth 166318 16.63% 16.67% gimel 143040 14.30% 14.29% daleth 125136 12.51% 12.50% he 111251 11.13% 11.11% waw 99946 9.99% 10.00% zayin 90974 9.10% 9.09% heth 63008 6.30% 6.35% ``` ## Nim ```nim import tables, math, strutils, times const num_trials = 1000000 precsn = 6 var start = cpuTime() var probs = initTable[string,float](16) probs.add("aleph", 1/5.0) probs.add("beth", 1/6.0) probs.add("gimel", 1/7.0) probs.add("daleth", 1/8.0) probs.add("he", 1/9.0) probs.add("waw", 1/10.0) probs.add("zayin", 1/11.0) probs.add("heth", 1759/27720) var samples = initTable[string,int](16) for i, j in pairs(probs): samples.add(i,0) randomize() for i in 1 .. num_trials: var z = random(1.0) for j,k in pairs(probs): if z < probs[j]: samples[j] = samples[j] + 1 break else: z = z - probs[j] var s1, s2: float echo("Item ","\t","Target ","\t","Results ","\t","Difference") echo("==== ","\t","====== ","\t"," ### = ","\t"," ### ==== ") for i, j in pairs(probs): s1 += samples[i]/num_trials*100.0 s2 += probs[i]*100.0 echo( i, "\t", formatFloat(probs[i],ffDecimal,precsn), "\t", formatFloat(samples[i]/num_trials,ffDecimal,precsn), "\t", formatFloat(100.0*(1.0-(samples[i]/num_trials)/probs[i]),ffDecimal,precsn),"%") echo("======","\t"," ### = ","\t"," ### == ") echo("Total:","\t",formatFloat(s2,ffDecimal,2)," \t",formatFloat(s1,ffDecimal,2)) echo("\n",formatFloat(cpuTime()-start,ffDecimal,2)," secs") ``` {{out}} ```txt Item Target Results Difference ==== ====== ### = ### ==== he 0.111111 0.110760 0.316000% heth 0.063456 0.063777 -0.505881% beth 0.166667 0.166386 0.168400% aleph 0.200000 0.200039 -0.019500% zayin 0.090909 0.090923 -0.015300% waw 0.100000 0.100513 -0.513000% gimel 0.142857 0.142691 0.116300% daleth 0.125000 0.124911 0.071200% ====== ### = ### == Total: 100.00 100.00 7.06 secs ``` ## OCaml ```ocaml let p = [ "Aleph", 1.0 /. 5.0; "Beth", 1.0 /. 6.0; "Gimel", 1.0 /. 7.0; "Daleth", 1.0 /. 8.0; "He", 1.0 /. 9.0; "Waw", 1.0 /. 10.0; "Zayin", 1.0 /. 11.0; "Heth", 1759.0 /. 27720.0; ] let rec take k = function | (v, p)::tl -> if k < p then v else take (k -. p) tl | _ -> invalid_arg "take" let () = let n = 1_000_000 in Random.self_init(); let h = Hashtbl.create 3 in List.iter (fun (v, _) -> Hashtbl.add h v 0) p; let tot = List.fold_left (fun acc (_, p) -> acc +. p) 0.0 p in for i = 1 to n do let sel = take (Random.float tot) p in let n = Hashtbl.find h sel in Hashtbl.replace h sel (succ n) (* count the number of each item *) done; List.iter (fun (v, p) -> let d = Hashtbl.find h v in Printf.printf "%s \t %f %f\n" v p (float d /. float n) ) p ``` Output: ```txt Aleph 0.200000 0.200272 Beth 0.166667 0.166381 Gimel 0.142857 0.142497 Daleth 0.125000 0.125005 He 0.111111 0.111272 Waw 0.100000 0.100069 Zayin 0.090909 0.091136 Heth 0.063456 0.063368 ``` ## PARI/GP ```parigp pc()={ my(v=[5544,10164,14124,17589,20669,23441,25961,27720],u=vector(8),e); for(i=1,1e6, my(r=random(27720)); for(j=1,8, if(r 1e6; sub prob_choice_picker { my %options = @_; my ($n, @a) = 0; while (my ($k,$v) = each %options) { $n += $v; push @a, [$n, $k]; } return sub { my $r = rand; ( first {$r <= $_->[0]} @a )->[1]; }; } my %ps = (aleph => 1/5, beth => 1/6, gimel => 1/7, daleth => 1/8, he => 1/9, waw => 1/10, zayin => 1/11); $ps{heth} = 1 - sum values %ps; my $picker = prob_choice_picker %ps; my %results; for (my $n = 0 ; $n < TRIALS ; ++$n) { ++$results{$picker->()}; } print "Event Occurred Expected Difference\n"; foreach (sort {$results{$b} <=> $results{$a}} keys %results) { printf "%-6s %f %f %f\n", $_, $results{$_}/TRIALS, $ps{$_}, abs($results{$_}/TRIALS - $ps{$_}); } ``` Sample output: ```txt Event Occurred Expected Difference aleph 0.198915 0.200000 0.001085 beth 0.166804 0.166667 0.000137 gimel 0.142992 0.142857 0.000135 daleth 0.125155 0.125000 0.000155 he 0.111160 0.111111 0.000049 waw 0.100229 0.100000 0.000229 zayin 0.091014 0.090909 0.000105 heth 0.063731 0.063456 0.000275 ``` ## Perl 6 {{works with|Rakudo|2018.10}} ```perl6 constant TRIALS = 1e6; constant @event = ; constant @P = flat (1 X/ 5 .. 11), 1759/27720; constant @cP = [\+] @P; my atomicint @results[+@event]; (^TRIALS).race.map: { @results[ @cP.first: { $_ > once rand }, :k ]⚛++; } say 'Event Occurred Expected Difference'; for ^@results { my ($occurred, $expected) = @results[$_], @P[$_] * TRIALS; printf "%-9s%8.0f%9.1f%12.1f\n", @event[$_], $occurred, $expected, abs $occurred - $expected; } ``` {{out}} ```txt Event Occurred Expected Difference aleph 200369 200000.0 369.0 beth 167005 166666.7 338.3 gimel 142690 142857.1 167.1 daleth 125061 125000.0 61.0 he 110563 111111.1 548.1 waw 100214 100000.0 214.0 zayin 90617 90909.1 292.1 heth 63481 63456.0 25.0 ``` ## Phix ```Phix constant {names, probs} = columnize({{"aleph", 1/5}, {"beth", 1/6}, {"gimel", 1/7}, {"daleth", 1/8}, {"he", 1/9}, {"waw", 1/10}, {"zayin", 1/11}, {"heth", 1759/27720}}) sequence results = repeat(0,length(names)) atom r constant lim = 1000000 for j=1 to lim do r = rnd() for i=1 to length(probs) do r -= probs[i] if r<=0 then results[i]+=1 exit end if end for end for printf(1," Name Actual Expected\n") for i=1 to length(probs) do printf(1,"%6s %8.6f %8.6f\n",{names[i],results[i]/lim,probs[i]}) end for ``` {{out}} ```txt Name Actual Expected aleph 0.201010 0.200000 beth 0.166311 0.166667 gimel 0.143354 0.142857 daleth 0.124841 0.125000 he 0.110544 0.111111 waw 0.100228 0.100000 zayin 0.090270 0.090909 heth 0.063442 0.063456 ``` ## PicoLisp ```PicoLisp (let (Count 1000000 Denom 27720 N Denom) (let Probs (mapcar '((I S) (prog1 (cons N (*/ Count I) 0 S) (dec 'N (/ Denom I)) ) ) (range 5 12) '(aleph beth gimel daleth he waw zayin heth) ) (do Count (inc (cddr (rank (rand 1 Denom) Probs T))) ) (let Fmt (-6 12 12) (tab Fmt NIL "Probability" "Result") (for X Probs (tab Fmt (cdddr X) (format (cadr X) 6) (format (caddr X) 6) ) ) ) ) ) ``` Output: ```txt Probability Result aleph 0.200000 0.199760 beth 0.166667 0.166878 gimel 0.142857 0.142977 daleth 0.125000 0.124983 he 0.111111 0.111200 waw 0.100000 0.100173 zayin 0.090909 0.090591 heth 0.083333 0.063438 ``` ## PL/I ```pli probch: Proc Options(main); Dcl prob(8) Dec Float(15) Init((1/5.0), /* aleph */ (1/6.0), /* beth */ (1/7.0), /* gimel */ (1/8.0), /* daleth */ (1/9.0), /* he */ (1/10.0), /* waw */ (1/11.0), /* zayin */ (1759/27720));/* heth */ Dcl what(8) Char(6) Init('aleph ','beth ','gimel ','daleth', 'he ','waw ','zayin ','heth '); Dcl ulim(0:8) Dec Float(15) Init((9)0); Dcl i Bin Fixed(31); Dcl ifloat Dec Float(15); Dcl one Dec Float(15) Init(1); Dcl num Dec Float(15) Init(1759); Dcl denom Dec Float(15) Init(27720); Dcl x Dec Float(15) Init(0); Dcl pr Dec Float(15) Init(0); Dcl (n,nn) Bin Fixed(31); Dcl cnt(8) Bin Fixed(31) Init((8)0); nn=1000000; Do i=1 To 8; ifloat=i+4; If i<8 Then prob(i)=one/ifloat; Else prob(i)=num/denom; Ulim(i)=ulim(i-1)+prob(i); /* Put Skip list(i,prob(i),ulim(i));*/ End; Do n=1 To nn; x=random(); Do i=1 To 8; If x 1 then some choices may be impossible (define (probabalistic-choice choices) (let-values (((_ choice) ;; the fold provides two values, we only need the second ;; the first will always be a negative number showing that ;; I've run out of random steam (for/fold ((rnd (random)) (choice #f)) (((v p) choices) #:break (<= rnd 0)) (values (- rnd p) v)))) choice)) ;;; ditto, but all probabilities must be exact rationals ;;; the optional lcd ;;; ;;; not the most efficient, since it provides a wrapper (and demo) ;;; for p-c/i-w below (define (probabalistic-choice/exact choices #:gcd (GCD (/ (apply gcd (hash-values choices))))) (probabalistic-choice/integer-weights (for/hash (((k v) choices)) (values k (* v GCD))) #:sum-of-weights GCD)) ;;; this proves useful in Rock-Paper-Scissors (define (probabalistic-choice/integer-weights choices #:sum-of-weights (sum-of-weights (apply + (hash-values choices)))) (let-values (((_ choice) (for/fold ((rnd (random sum-of-weights)) (choice #f)) (((v p) choices) #:break (< rnd 0)) (values (- rnd p) v)))) choice)) (module+ test (define test-samples (make-parameter 1000000)) (define (test-p-c-function f w) (define test-selection (make-hash)) (for* ((i (in-range 0 (test-samples))) (c (in-value (f w)))) (when (zero? (modulo i 100000)) (eprintf "~a," (quotient i 100000))) (hash-update! test-selection c add1 0)) (printf "~a~%choice\tcount\texpected\tratio\terror~%" f) (for* (((k v) (in-hash test-selection)) (e (in-value (* (test-samples) (hash-ref w k))))) (printf "~a\t~a\t~a\t~a\t~a%~%" k v e (/ v (test-samples)) (real->decimal-string (exact->inexact (* 100 (/ (- v e) e))))))) (define test-weightings/rosetta (hash 'aleph 1/5 'beth 1/6 'gimel 1/7 'daleth 1/8 'he 1/9 'waw 1/10 'zayin 1/11 'heth 1759/27720; adjusted so that probabilities add to 1 )) (define test-weightings/50:50 (hash 'woo 1/2 'yay 1/2)) (define test-weightings/1:2:3 (hash 'woo 1 'yay 2 'foo 3)) (test-p-c-function probabalistic-choice test-weightings/50:50) (test-p-c-function probabalistic-choice/exact test-weightings/50:50) (test-p-c-function probabalistic-choice test-weightings/rosetta) (test-p-c-function probabalistic-choice/exact test-weightings/rosetta)) ``` Output (note that the progress counts, which go to standard error, are interleaved with the output on standard out) ```txt 0,1,2,3,4,5,6,7,8,9,# choice count expected ratio error yay 499744 500000 15617/31250 -0.05% woo 500256 500000 15633/31250 0.05% 0,1,2,3,4,5,6,7,8,9,# choice count expected ratio error yay 499852 500000 124963/250000 -0.03% woo 500148 500000 125037/250000 0.03% 0,1,2,3,4,5,6,7,8,9,# choice count expected ratio error daleth 124964 125000 31241/250000 -0.03% zayin 90233 1000000/11 90233/1000000 -0.74% gimel 142494 1000000/7 71247/500000 -0.25% heth 64045 43975000/693 12809/200000 0.93% aleph 199690 200000 19969/100000 -0.15% beth 166861 500000/3 166861/1000000 0.12% waw 100075 100000 4003/40000 0.07% he 111638 1000000/9 55819/500000 0.47% 0,1,2,3,4,5,6,7,8,9,# choice count expected ratio error beth 166423 500000/3 166423/1000000 -0.15% heth 63462 43975000/693 31731/500000 0.01% daleth 125091 125000 125091/1000000 0.07% waw 99820 100000 4991/50000 -0.18% aleph 200669 200000 200669/1000000 0.33% gimel 142782 1000000/7 71391/500000 -0.05% zayin 90478 1000000/11 45239/500000 -0.47% he 111275 1000000/9 4451/40000 0.15% ``` ## REXX Note: REXX can generate random numbers up to 100,000. ```rexx /*REXX program displays results of probabilistic choices, gen random #s per probability.*/ parse arg trials digs seed . /*obtain the optional arguments from CL*/ if trials=='' | trials=="," then trials=1000000 /*Not specified? Then use the default.*/ if digs=='' | digs=="," then digs=15 /* " " " " " " */ if datatype(seed, 'W') then call random ,,seed /*allows repeatability for RANDOM nums.*/ numeric digits digs /*use a specific number of decimal digs*/ names= 'aleph beth gimel daleth he waw zayin heth ───totals───►' /*names of the cells.*/ HI=100000 /*max REXX RANDOM num*/ z=words(names); #=z - 1 /*#≡the number of actual/useable names.*/ $=0 /*initialize sum of the probabilities. */ do n=1 for #; prob.n=1 / (n+4); if n==# then prob.n= 1759 / 27720 $=$ + prob.n; Hprob.n=prob.n * HI end /*n*/ prob.z=$ /*define the value of the ───totals───.*/ @.=0 /*initialize all counters in the range.*/ @.z=trials /*define the last counter of " " */ do j=1 for trials; r=random(HI) /*gen TRIAL number of random numbers.*/ do k=1 for # /*for each cell, compute percentages. */ if r<=Hprob.k then @.k=@.k + 1 /* " " " range, bump the counter*/ end /*k*/ end /*j*/ _= '═' /*_: a literal used for CENTER BIF pad*/ w=digs + 6 /*W: display width for the percentages*/ d=4 + max( length(trials), length('count') ) /* [↓] display a formatted top header.*/ say center('name',15,_) center('count',d,_) center('target %',w,_) center('actual %',w,_) do cell=1 for z /*display each of the cells and totals.*/ say ' ' left( word(names, cell), 13) right(@.cell, d-2) " ", left( format( prob.cell * 100, d), w-2), left( format( @.cell/trials * 100, d), w-2) if cell==# then say center(_,15,_) center(_,d,_) center(_,w,_) center(_,w,_) end /*c*/ /* [↑] display a formatted foot header*/ /*stick a fork in it, we are all done.*/ ``` {{out|output|text= when using the default input:}} ```txt ═════name══════ ═══count═══ ══════target %═══════ ══════actual %═══════ aleph 200135 20 20.0135 beth 166912 16.6666666 16.6912 gimel 143222 14.2857142 14.3222 daleth 124991 12.5 12.4991 he 111259 11.1111111 11.1259 waw 100049 10 10.0049 zayin 90978 9.0909090 9.0978 heth 63278 6.3455988 6.3278 ═══════════════ ═══════════ ═════════════════════ ═════════════════════ ───totals───► 1000000 100 100 ``` ## Ring ```ring # Project : Probabilistic choice cnt = list(8) item = ["aleph","beth","gimel","daleth","he","waw","zayin","heth"] prob = [1/5.0, 1/6.0, 1/7.0, 1/8.0, 1/9.0, 1/10.0, 1/11.0, 1759/27720] for trial = 1 to 1000000 r = random(10)/10 p = 0 for i = 1 to len(prob) p = p + prob[i] if r < p cnt[i] = cnt[i] + 1 loop ok next next see "item actual theoretical" + nl for i = 1 to len(item) see "" + item[i] + " " + cnt[i]/1000000 + " " + prob[i] + nl next ``` Output: ```txt item actual theoretical aleph 0.091307 0.200000 beth 0.181073 0.166667 gimel 0.181884 0.142857 daleth 0.090985 0.125000 he 0.090958 0.111111 waw 0.091064 0.100000 zayin 0.091061 0.090909 heth 0 0.063456 ``` ## Ruby ```ruby probabilities = { "aleph" => 1/5.0, "beth" => 1/6.0, "gimel" => 1/7.0, "daleth" => 1/8.0, "he" => 1/9.0, "waw" => 1/10.0, "zayin" => 1/11.0, } probabilities["heth"] = 1.0 - probabilities.each_value.inject(:+) ordered_keys = probabilities.keys sum, sums = 0.0, {} ordered_keys.each do |key| sum += probabilities[key] sums[key] = sum end actual = Hash.new(0) samples = 1_000_000 samples.times do r = rand for k in ordered_keys if r < sums[k] actual[k] += 1 break end end end puts "key expected actual diff" for k in ordered_keys act = Float(actual[k]) / samples val = probabilities[k] printf "%-8s%.8f %.8f %6.3f %%\n", k, val, act, 100*(act-val)/val end ``` {{out}} ```txt key expected actual diff aleph 0.20000000 0.19949200 -0.254 % beth 0.16666667 0.16689900 0.139 % gimel 0.14285714 0.14309300 0.165 % daleth 0.12500000 0.12494200 -0.046 % he 0.11111111 0.11037800 -0.660 % waw 0.10000000 0.10030100 0.301 % zayin 0.09090909 0.09162700 0.790 % heth 0.06345599 0.06326800 -0.296 % ``` ## Rust ```Rust extern crate rand; use rand::distributions::{IndependentSample, Sample, Weighted, WeightedChoice}; use rand::{weak_rng, Rng}; const DATA: [(&str, f64); 8] = [ ("aleph", 1.0 / 5.0), ("beth", 1.0 / 6.0), ("gimel", 1.0 / 7.0), ("daleth", 1.0 / 8.0), ("he", 1.0 / 9.0), ("waw", 1.0 / 10.0), ("zayin", 1.0 / 11.0), ("heth", 1759.0 / 27720.0), ]; const SAMPLES: usize = 1_000_000; /// Generate a mapping to be used by `WeightedChoice` fn gen_mapping() -> Vec> { DATA.iter() .enumerate() .map(|(i, &(_, p))| Weighted { // `WeightedChoice` requires `u32` weights rather than raw probabilities. For each // probability, we convert it to a `u32` weight, and associate it with an index. We // multiply by a constant because small numbers such as 0.2 when casted to `u32` // become `0`. This conversion decreases the accuracy of the mapping, which is why we // provide an implementation which uses `f64`s for the best accuracy. weight: (p * 1_000_000_000.0) as u32, item: i, }) .collect() } /// Generate a mapping of the raw probabilities fn gen_mapping_float() -> Vec { // This does the work of `WeightedChoice::new`, splitting a number into various ranges. The // `item` of `Weighted` is represented here merely by the probability's position in the `Vec`. let mut running_total = 0.0; DATA.iter() .map(|&(_, p)| { running_total += p; running_total }) .collect() } /// An implementation of `WeightedChoice` which uses probabilities rather than weights. Refer to /// the `WeightedChoice` source for serious usage. struct WcFloat { mapping: Vec, } impl WcFloat { fn new(mapping: &[f64]) -> Self { Self { mapping: mapping.to_vec(), } } // This is roughly the same logic as `WeightedChoice::ind_sample` (though is likely slower) fn search(&self, sample_prob: f64) -> usize { let idx = self.mapping .binary_search_by(|p| p.partial_cmp(&sample_prob).unwrap()); match idx { Ok(i) | Err(i) => i, } } } impl IndependentSample for WcFloat { fn ind_sample(&self, rng: &mut R) -> usize { // Because we know the total is exactly 1.0, we can merely use a raw float value. // Otherwise caching `Range::new(0.0, running_total)` and sampling with // `range.ind_sample(&mut rng)` is recommended. let sample_prob = rng.next_f64(); self.search(sample_prob) } } impl Sample for WcFloat { fn sample(&mut self, rng: &mut R) -> usize { self.ind_sample(rng) } } fn take_samples(rng: &mut R, wc: &T) -> [usize; 8] where T: IndependentSample, { let mut counts = [0; 8]; for _ in 0..SAMPLES { let sample = wc.ind_sample(rng); counts[sample] += 1; } counts } fn print_mapping(counts: &[usize]) { println!("Item | Expected | Actual "); println!("-------+----------+----------"); for (&(name, expected), &count) in DATA.iter().zip(counts.iter()) { let real = count as f64 / SAMPLES as f64; println!("{:6} | {:.6} | {:.6}", name, expected, real); } } fn main() { let mut rng = weak_rng(); println!(" ~~~ U32 METHOD ~~~"); let mut mapping = gen_mapping(); let wc = WeightedChoice::new(&mut mapping); let counts = take_samples(&mut rng, &wc); print_mapping(&counts); println!(); println!(" ~~~ FLOAT METHOD ~~~"); // initialize the float version of `WeightedChoice` let mapping = gen_mapping_float(); let wc = WcFloat::new(&mapping); let counts = take_samples(&mut rng, &wc); print_mapping(&counts); } ``` {{out}} ```txt ~~~ U32 METHOD ~~~ Item | Expected | Actual -------+----------+---------- aleph | 0.200000 | 0.200195 beth | 0.166667 | 0.166182 gimel | 0.142857 | 0.142502 daleth | 0.125000 | 0.125503 he | 0.111111 | 0.110820 waw | 0.100000 | 0.100166 zayin | 0.090909 | 0.090927 heth | 0.063456 | 0.063705 ~~~ FLOAT METHOD ~~~ Item | Expected | Actual -------+----------+---------- aleph | 0.200000 | 0.199984 beth | 0.166667 | 0.166634 gimel | 0.142857 | 0.143218 daleth | 0.125000 | 0.124956 he | 0.111111 | 0.111047 waw | 0.100000 | 0.099805 zayin | 0.090909 | 0.090513 heth | 0.063456 | 0.063843 ``` ## Scala This algorithm consists of a concise two-line tail-recursive loop (def weighted). The rest of the code is for API robustness, testing and display. weightedProb is for the task as stated (0 < p < 1), and weightedFreq is the equivalent based on integer frequencies (f >= 0). ```Scala object ProbabilisticChoice extends App { import scala.collection.mutable.LinkedHashMap def weightedProb[A](prob: LinkedHashMap[A,Double]): A = { require(prob.forall{case (_, p) => p > 0 && p < 1}) assume(prob.values.sum == 1) def weighted(todo: Iterator[(A,Double)], rand: Double, accum: Double = 0): A = todo.next match { case (s, i) if rand < (accum + i) => s case (_, i) => weighted(todo, rand, accum + i) } weighted(prob.toIterator, scala.util.Random.nextDouble) } def weightedFreq[A](freq: LinkedHashMap[A,Int]): A = { require(freq.forall{case (_, f) => f >= 0}) require(freq.values.sum > 0) def weighted(todo: Iterator[(A,Int)], rand: Int, accum: Int = 0): A = todo.next match { case (s, i) if rand < (accum + i) => s case (_, i) => weighted(todo, rand, accum + i) } weighted(freq.toIterator, scala.util.Random.nextInt(freq.values.sum)) } // Tests: val probabilities = LinkedHashMap( 'aleph -> 1.0/5, 'beth -> 1.0/6, 'gimel -> 1.0/7, 'daleth -> 1.0/8, 'he -> 1.0/9, 'waw -> 1.0/10, 'zayin -> 1.0/11, 'heth -> 1759.0/27720 ) val frequencies = LinkedHashMap( 'aleph -> 200, 'beth -> 167, 'gimel -> 143, 'daleth -> 125, 'he -> 111, 'waw -> 100, 'zayin -> 91, 'heth -> 63 ) def check[A](original: LinkedHashMap[A,Double], results: Seq[A]) { val freq = results.groupBy(x => x).mapValues(_.size.toDouble/results.size) original.foreach{case (k, v) => val a = v/original.values.sum val b = freq(k) val c = if (Math.abs(a - b) < 0.001) "ok" else "**" println(f"$k%10s $a%.4f $b%.4f $c") } println(" "*10 + f" ${1}%.4f ${freq.values.sum}%.4f") } println("Checking weighted probabilities:") check(probabilities, for (i <- 1 to 1000000) yield weightedProb(probabilities)) println println("Checking weighted frequencies:") check(frequencies.map{case (a, b) => a -> b.toDouble}, for (i <- 1 to 1000000) yield weightedFreq(frequencies)) } ``` {{out}} ```txt Checking weighted probabilities: 'aleph 0.2000 0.2001 ok 'beth 0.1667 0.1665 ok 'gimel 0.1429 0.1430 ok 'daleth 0.1250 0.1248 ok 'he 0.1111 0.1112 ok 'waw 0.1000 0.1000 ok 'zayin 0.0909 0.0911 ok 'heth 0.0635 0.0632 ok 1.0000 1.0000 Checking weighted frequencies: 'aleph 0.2000 0.2000 ok 'beth 0.1670 0.1672 ok 'gimel 0.1430 0.1432 ok 'daleth 0.1250 0.1243 ok 'he 0.1110 0.1105 ok 'waw 0.1000 0.1002 ok 'zayin 0.0910 0.0913 ok 'heth 0.0630 0.0632 ok 1.0000 1.0000 ``` ## Seed7 To reduce the runtime this program should be compiled. ```seed7 $ include "seed7_05.s7i"; include "float.s7i"; const type: letter is new enum aleph, beth, gimel, daleth, he, waw, zayin, heth end enum; const func string: str (in letter: aLetter) is return [] ("aleph", "beth", "gimel", "daleth", "he", "waw", "zayin", "heth") [succ(ord(aLetter))]; enable_output(letter); const array [letter] integer: table is [letter] ( 5544, 4620, 3960, 3465, 3080, 2772, 2520, 1759); const func letter: randomLetter is func result var letter: resultLetter is aleph; local var integer: number is 0; begin number := rand(1, 27720); while number > table[resultLetter] do number -:= table[resultLetter]; incr(resultLetter); end while; end func; const proc: main is func local var integer: count is 0; var letter: aLetter is aleph; var array [letter] integer: occurrence is letter times 0; begin for count range 1 to 1000000 do aLetter := randomLetter; incr(occurrence[aLetter]); end for; writeln("Name Count Ratio Expected"); for aLetter range letter.first to letter.last do writeln(aLetter rpad 7 <& occurrence[aLetter] lpad 6 <& flt(occurrence[aLetter]) / 10000.9 digits 4 lpad 8 <& "%" <& 100.0 * flt(table[aLetter]) / 27720.0 digits 4 lpad 8 <& "%"); end for; end func; ``` Outout: ```txt Name Count Ratio Expected aleph 199788 19.9770% 20.0000% beth 166897 16.6882% 16.6667% gimel 143103 14.3090% 14.2857% daleth 125060 12.5049% 12.5000% he 110848 11.0838% 11.1111% waw 99550 9.9541% 10.0000% zayin 90918 9.0910% 9.0909% heth 63836 6.3830% 6.3456% ``` ## Scheme Using guile scheme 2.0.11. ```scheme (use-modules (ice-9 format)) (define (random-choice probs) (define choice (random 1.0)) (define (helper val prob-lis) (let ((nval (- val (cadar prob-lis)))) (if (< nval 0) (caar prob-lis) (helper nval (cdr prob-lis))))) (helper choice probs)) (define (add-result result delta table) (cond ((null? table) (list (list result delta))) ((eq? (caar table) result) (cons (list result (+ (cadar table) delta)) (cdr table))) (#t (cons (car table) (add-result result delta (cdr table)))))) (define (choices trials probs) (define (helper trial-num freq-table) (if (= trial-num trials) freq-table (helper (+ trial-num 1) (add-result (random-choice probs) (/ 1 trials) freq-table)))) (helper 0 '())) (define (format-results probs results) (for-each (lambda (x) (format #t "~10a~10,5f~10,5f~%" (car x) (cadr x) (cadr (assoc (car x) results)))) probs)) (define probs '((aleph 1/5) (beth 1/6) (gimel 1/7) (daleth 1/8) (he 1/9) (waw 1/10) (zayin 1/11) (heth 1759/27720))) (format-results probs (choices 1000000 probs)) ``` Example output: aleph 0.20000 0.20051 beth 0.16667 0.16680 gimel 0.14286 0.14231 daleth 0.12500 0.12538 he 0.11111 0.11136 waw 0.10000 0.09955 zayin 0.09091 0.09096 heth 0.06346 0.06313 ## Sidef {{trans|Perl}} ```ruby define TRIALS = 1e4;   func prob_choice_picker(options) { var n = 0; var a = []; options.each { |k,v| n += v; a << [n, k]; } func { var r = 1.rand; a.first{|e| r <= e[0] }[1]; } }   var ps = Hash( aleph => 1/5, beth => 1/6, gimel => 1/7, daleth => 1/8, he => 1/9, waw => 1/10, zayin => 1/11 )   ps{:heth} = (1 - ps.values.sum)   var picker = prob_choice_picker(ps) var results = Hash()   TRIALS.times { results{picker()} := 0 ++; }   say "Event Occurred Expected Difference"; for k,v in (results.sort_by {|k| results{k} }.reverse) { printf("%-6s  %f  %f  %f\n", k, v/TRIALS, ps{k}, abs(v/TRIALS - ps{k}) ); } ``` {{out}} ```txt Event Occurred Expected Difference aleph 0.196300 0.200000 0.003700 beth 0.165600 0.166667 0.001067 gimel 0.143700 0.142857 0.000843 daleth 0.123900 0.125000 0.001100 he 0.111800 0.111111 0.000689 waw 0.101900 0.100000 0.001900 zayin 0.088100 0.090909 0.002809 heth 0.068800 0.063456 0.005344 ``` ## Stata ```stata clear mata letters="aleph","beth","gimel","daleth","he","waw","zayin","heth" a=letters[rdiscrete(10000,1,(1/5,1/6,1/7,1/8,1/9,1/10,1/11,1759/27720))]' st_addobs(10000) st_addvar("str10","a") st_sstore(.,.,a) end ``` ## Tcl ```tcl package require Tcl 8.5 set map [dict create] set sum 0.0 foreach name {aleph beth gimel daleth he waw zayin} \ prob {1/5.0 1/6.0 1/7.0 1/8.0 1/9.0 1/10.0 1/11.0} \ { set prob [expr $prob] set sum [expr {$sum + $prob}] dict set map $name [dict create probability $prob limit $sum count 0] } dict set map heth [dict create probability [expr {1.0 - $sum}] limit 1.0 count 0] set samples 1000000 for {set i 0} {$i < $samples} {incr i} { set n [expr {rand()}] foreach name [dict keys $map] { if {$n <= [dict get $map $name limit]} { set count [dict get $map $name count] dict set map $name count [incr count] break } } } puts "using $samples samples:" puts [format "%-10s %-21s %-9s %s" "" expected actual difference] dict for {name submap} $map { dict with submap { set actual [expr {$count * 1.0 / $samples}] puts [format "%-10s %-21s %-9s %4.2f%%" $name $probability $actual \ [expr {abs($actual - $probability)/$probability*100.0}] ] } } ``` ```txt using 1000000 samples: expected actual difference aleph 0.2 0.199641 0.18% beth 0.16666666666666666 0.1674 0.44% gimel 0.14285714285714285 0.143121 0.18% daleth 0.125 0.124864 0.11% he 0.1111111111111111 0.111036 0.07% waw 0.1 0.100021 0.02% zayin 0.09090909090909091 0.09018 0.80% heth 0.06345598845598843 0.063737 0.44% ``` ## Ursala The stochasm library function used here constructs a weighted non-deterministic choice of a set of functions. The pseudo-random number generator is a 64 bit Mersenne twistor implemented by the run time system. ```Ursala #import std #import nat #import flo outcomes = <'aleph ','beth ','gimel ','daleth','he ','waw ','zayin ','heth '> probabilities = ^lrNCT(~&,minus/1.+ plus:-0) div/*1. float* skip/5 iota12 simulation = ^(~&rn,div+ float~~rmPlX)^*D/~& iota; ^A(~&h,length)*K2+ * stochasm@p/probabilities !* outcomes format = :/' frequency probability'+ * ^lrlrTPT/~&n (printf/'%12.8f')^~/~&m outcomes-$probabilities@n #show+ results = format simulation 1000000 ``` output: ```txt frequency probability daleth 0.12484500 0.12500000 beth 0.16680600 0.16666667 aleph 0.19973700 0.20000000 waw 0.10016900 0.10000000 gimel 0.14293100 0.14285714 he 0.11131100 0.11111111 zayin 0.09104700 0.09090909 heth 0.06315400 0.06345599 ``` ## VBScript Derived from the BBC BASIC version ```vb item = Array("aleph","beth","gimel","daleth","he","waw","zayin","heth") prob = Array(1/5.0, 1/6.0, 1/7.0, 1/8.0, 1/9.0, 1/10.0, 1/11.0, 1759/27720) Dim cnt(7) 'Terminate script if sum of probabilities <> 1. sum = 0 For i = 0 To UBound(prob) sum = sum + prob(i) Next If sum <> 1 Then WScript.Quit End If For trial = 1 To 1000000 r = Rnd(1) p = 0 For i = 0 To UBound(prob) p = p + prob(i) If r < p Then cnt(i) = cnt(i) + 1 Exit For End If Next Next WScript.StdOut.Write "item" & vbTab & "actual" & vbTab & vbTab & "theoretical" WScript.StdOut.WriteLine For i = 0 To UBound(item) WScript.StdOut.Write item(i) & vbTab & FormatNumber(cnt(i)/1000000,6) & vbTab & FormatNumber(prob(i),6) WScript.StdOut.WriteLine Next ``` {{out}} ```txt item actual theoretical aleph 0.199755 0.200000 beth 0.166861 0.166667 gimel 0.143240 0.142857 daleth 0.124474 0.125000 he 0.110879 0.111111 waw 0.100341 0.100000 zayin 0.090745 0.090909 heth 0.063705 0.063456 ``` ## XPL0 ```XPL0 include c:\cxpl\codes; def Size = 10_000_000; int Tbl(12+1); int I, J, N; real X, S0, S1; [for J:= 5 to 12 do Tbl(J):= 0; for I:= 0 to 1_000_000-1 do \generate one million items [N:= Ran(Size); for J:= 5 to 11 do [N:= N - Size/J; if N < 0 then [Tbl(J):= Tbl(J)+1; J:= 100]; ]; if J=12 then Tbl(12):= Tbl(12)+1; ]; S0:= 0.0; S1:= 0.0; for J:= 5 to 11 do [X:= 1.0/float(J); RlOut(0, X); S0:= S0+X; X:= float(Tbl(J)) / 1_000_000.0; RlOut(0, X); S1:= S1+X; CrLf(0); ]; X:= 1759.0 / 27720.0; RlOut(0, X); S0:= S0+X; X:= float(Tbl(12)) / 1_000_000.0; RlOut(0, X); S1:= S1+X; CrLf(0); Text(0, " ------- ------- "); RlOut(0, S0); RlOut(0, S1); ] ``` Output: ```txt 0.20000 0.20012 0.16667 0.16679 0.14286 0.14305 0.12500 0.12510 0.11111 0.11113 0.10000 0.09990 0.09091 0.09077 0.06346 0.06313 ------- ------- 1.00000 1.00000 ``` ## zkl {{trans|C}} ```zkl var names=T("aleph", "beth", "gimel", "daleth", "he", "waw", "zayin", "heth"); var ptable=T(5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0).apply('/.fp(1.0)); ptable=ptable.append(1.0-ptable.sum(0.0)); // add last weight to sum to 1.0 var [const] N=ptable.len(); fcn ridx{ i:=0; s:=(0.0).random(1); while((s-=ptable[i]) > 0) { i+=1 } i } const M=0d1_000_000; var r=(0).pump(N,List,T(Ref,0)); // list of references to int 0 (0).pump(M,Void,fcn{r[ridx()].inc()}); // 1,000,000 weighted random #s r=r.apply("value").apply("toFloat"); // (reference to int)-->int-->float println(" Name Count Ratio Expected"); foreach i in (N){ "%6s%7d %7.4f%% %7.4f%%".fmt(names[i], r[i], r[i]/M*100, ptable[i]*100).println(); } ``` {{out}} ```txt Name Count Ratio Expected aleph 200214 20.0214% 20.0000% beth 166399 16.6399% 16.6667% gimel 143100 14.3100% 14.2857% daleth 125197 12.5197% 12.5000% he 111167 11.1167% 11.1111% waw 100097 10.0097% 10.0000% zayin 90692 9.0692% 9.0909% heth 63162 6.3162% 6.3456% ```