Statistics & Probability for the HP-41

Overview
 

 1°)  One Variable ( Grouped Data )

   1-a)  Moments of a Distribution  ( Program#1 )
   1-b)  Moments of a Distribution  ( Program#2 )
   1-c)  Arithmetic, Geometric and Harmonic Means

 2°)  Two Variables ( Grouped Data )

   2-a)  Standard Deviation, Covariance, Coefficient of Correlation
   2-b)  Linear Regression

 3°)  Probability Functions

   3-a)  Chi2-Probability Function
   3-b)  Gaussian Probability Function
   3-c)  F-Distribution Function

    3-c1)  Program#1
    3-c2)  Inverse FDF function

   3-d)  Student's t-Distribution & UTPT Functions

    3-d1)  Program#1
    3-d2)  UTPT function
    3-d3)  Inverse UTPT function

   3-e)  Poisson distribution
   3-f)  Binomial Distribution
   3-g) Negative Binomial Distribution
   3-h) Pascal Distribution
   3-i)  Geometric Distribution
   3-j)  Hypergeometric Distribution

 4°)  Classic Problems

  4-a)  Binomial Coefficients ( Focal Program & M-Code Routine )
  4-b)  Probability of no repetition in a sample
  4-c)  Hypergeometric Distribution ( Urn Problem )
  4-d)  Multivariate Hypergeometric Distribution
  4-c)  Permutations with repetitions
 

-See also "Quantiles for the HP-41" & "Least-Squares Approximation for the HP-41"

-The latest addition is paragraph 3-c2)
 

1°)  One Variable ( Grouped Data )
 

     a) Moments of a Distribution ( Program#1 )
 

-This program calculates the arithmetic mean, the standard deviation, the skewness and the kurtosis
  of a given set of p points  x₁ , x₂ , ............ ,  x_p
  with respective frequencies  n₁ , n₂ , ............ ,  n_p
 

Data Registers:    R00 = Sum n_i                                         when the routine stops,  R05 = m = arithmetic mean

                               R01 = Sum n_i x_i    R03 = Sum n_i x_i³
                               R02 = Sum n_i x_i²   R04 = Sum n_i x_i⁴

Flag:   F01                  CF 01  if you want to use the sample standard deviation
Subroutines: /            SF 01  if you want to use the population standard deviation
 
 

01  LBL A 02  ST+ 00 03  RCL Y 04  STO T 05  * 06  ST+ 01 07  * 08  ST+ 02 09  * 10  ST+ 03 11  * 12  ST+ 04 13  RDN 14  RTN 15  LBL "STAT"

16  RCL 01 17  RCL 00 18  / 19  ENTER^ 20  STO 05 21  6 22  * 23  RCL 02  24  * 25  RCL 03 26  4 27  * 28  - 29  * 30  RCL 04

31  + 32  RCL 00 33  / 34  R^ 35  X^2 36  X^2 37  3 38  * 39  - 40  RCL 03  41  R^ 42  RCL 02        43  * 44  3 45  *

46  - 47  RCL 00 48  / 49  RCL 05 50  3 51  Y^X 52  ST+ X 53  + 54  RCL 05 55  X^2 56  RCL 00  57  * 58  RCL 02        59  X<>Y 60  -

61  RCL 00 62  FC? 01 63  DSE X 64  / 65  ST/ Z 66  ST/ Z 67  ST/ Y 68  SQRT 69  ST/ Y 70  RCL 05        71  END

 
   ( 95 bytes / SIZE 006 )
 
 

STACK	INPUTS#i	OUTPUTS#i	final OUTPUTS
T	/	/	µ4
Z	/	/	µ3
Y	xi	/	s
X	ni	xi	m
L	/	/	V = s2

      m = arithmetic mean                                               V = s² = variance
       s = sample standard deviation if CF 01                        µ₃ = skewness
       s = population standard deviation if SF 01                   µ₄ = kurtosis

Example:    Given the following data
 
 

xi	2	7	9	12	15	19
ni	3	5	8	14	7	2

 
-Clear registers R00 thru R04  ( for instance execute CLS  if your statistics registers are R00 thru R05 )
-Then,  x_i  ENTER^   n_i S+  ( [A]  in user mode )  for all the points

   2  ENTER^      7  ENTER^      9  ENTER^     12  ENTER^    15  ENTER^     19  ENTER^
   3  S+               5  S+                8  S+              14  S+               7   S+               2   S+                in user mode

-Finally, simply press  R/S  ( or XEQ "STAT" )  it yields:

                                   CF 01                                                               SF 01

                                  m = 10.8718                                                       m = 10.8718
                       RDN    s  =  4.0012                                           RDN    s  =  3.9496
                       RDN   µ₃ = -0.3406 = skewness                        RDN   µ₃ = -0.3542 = skewness
                       RDN   µ₄ =  3.0605 =  kurtosis                          RDN   µ₄ =  3.2238 =  kurtosis
 

-To correct wrong inputs, key in  x_i  ENTER^   n_i  CHS  S+  ( [A] in user mode )  for the wrong ( x_i , n_i )
  or add   LBL a   CHS   just before  LBL A  and key in  x_i  ENTER^   n_iS-  ( [a]  in user mode )

-For a normal distribution,  µ₃ = 0  &  µ₄ =  3
 

     b) Moments of a Distribution ( Program#2 )
 

-This routine computes the centered moment  µ_k =  [ (1/n) Sum n_i.( x_i - m )^k ]/s^k   for a given k-value
  where  m is the arithmetic mean and s is the standard deviation.
  µ_k is a dimensionless quantity.
 

Data Registers:       •  R00 = p = number of points ,                  ( Registers R00 thru R2p  are to be initialized before executing "RCMK" )

                                  •  R01 = x₁     •  R03 = x₂      .......    •  R2p-1 = x_p
                                  •  R02 = n₁     •  R04 = n₂      .......    •  R2p   = n_p

Flag:   F01                  CF 01  if you want to use the sample standard deviation
Subroutines: /            SF 01  if you want to use the population standard deviation
 
 

01  LBL "RCMK" 02  CLA 03  STO O 04  RCL 00 05  ST+ X 06  STO M 07  ENTER^ 08  CLX 09  LBL 01 10  RCL IND Y 11  ST+ N 12  DSE Z

13  RCL IND Z 14  * 15  + 16  DSE Y 17  GTO 01 18  RCL N 19  / 20  X<>Y 21  LBL 02 22  RCL IND M  23  DSE M 24  RCL IND M

25  R^ 26  ST- Y 27  RDN 28  X^2 29  X<>Y 30  ST* Y 31  X<>Y 32  ST+ P 33  X<> L 34  X<>Y 35  X<> O           36  Y^X

37  LASTX 38  X<> O           39  * 40  + 41  DSE M 42  GTO 02 43  X<>Y 44  RCL P 45  RCL N  46  ST/ T 47  FC? 01 48  DSE X

49  / 50  SQRT 51  ENTER^ 52  X<> O           53  Y^X 54  ST/ Z 55  X<> O  56  X<> Z 57  CLA 58  END

 
   ( 97 bytes / SIZE 2p+1 )
 
 

STACK	INPUTS	OUTPUTS
T	/	mk = Sum ni ( xi - m )k/n
Z	/	s
Y	/	m
X	k	µk = mk/sk
L	/	k

         m = arithmetic mean
          s = sample standard deviation if CF 01
          s = population standard deviation if SF 01

Example:
 
 

xi	2	7	9	12	15	19
ni	3	5	8	14	7	2

 
-Store these 12 numbers into

           R01         R03         ..............................            R11
           R02         R04         ..............................            R12         respectively

-There are  p = 6  data points whence   6  STO 00

   CF 01    4   XEQ "RCMK"  >>>>  µ₄ =  3.0605   RDN   m = 10.8718   RDN   s =  4.0012   RDN   m₄ = 784.4261       (  µ₄ = kurtosis )
   SF 01    4           R/S            >>>>  µ₄ =  3.2238   RDN   m = 10.8718   RDN   s =  3.9496   RDN   m₄ = 784.4261

   CF 01    5           R/S           >>>>  µ₅ =  -2.2513   RDN   m = 10.8718   RDN   s =  4.0012   RDN   m₅ = -2308.7606    ( in 10 seconds )
   SF 01    5           R/S           >>>>  µ₅ =  -2.4024   RDN   m = 10.8718   RDN   s =  3.9496   RDN   m₅ = -2308.7606
 

     c)  Arithmetic, Geometric and Harmonic Means
 

-The arithmetic mean m is defined by   m = ( n₁ x₁ + n₂ x₂ + ...... + n_p x_p ) / ( n₁ + n₂ + .....  + n_p )
-The geometric mean m' is defined by  m' = [ ( x₁^n₁ ).( x₂^n₂ ) ...... ( x_p^n_p ) ] ^1/n                                        with   n =  n₁ + n₂ + .....  + n_p
-The harmonic mean m" is defined by   m" = ( n₁ + n₂ + .....  + n_p ) / ( n₁/x₁ + n₂/x₂ + ...... + n_p/x_p )
 

Data Registers:       •  R00 = p = number of points                   ( Registers R00 thru R2p  are to be initialized before executing "MEANS" )

                                  •  R01 = x₁     •  R03 = x₂      .......    •  R2p-1 = x_p
                                  •  R02 = n₁     •  R04 = n₂      .......    •  R2p   = n_p
Flags: /
Subroutines: /
 
 

01  LBL "MEANS" 02  RCL 00 03  ST+ X 04  1 05  CLA 06  LBL 01 07  RCL IND Y

08  STO O 09  ST+ P 10  DSE Z 11  RCL IND Z     12  / 13  ST+ N 14  CLX

15  RCL O  16  LASTX 17  * 18  ST+ M  19  X<> L 20  RCL O             21  Y^X

22  * 23  DSE Y 24  GTO 01           25  RCL P  26  ST/ M 27  ST/ N 28  1/X

29  Y^X 30  RCL N             31  1/X 32  X<>Y 33  RCL M  34  CLA 35  END

 
   ( 63 bytes / SIZE 2p+1 )
 
 

STACK	INPUTS	OUTPUTS
Z	/	m"
Y	/	m'
X	/	m

        m = arithmetic mean
        m' = geometric mean
        m" = harmonic mean

Example:
 
 

xi	2	7	9	12	15	19
ni	3	5	8	14	7	2

 
-Store these 12 numbers into

          R01         R03         ..............................            R11
          R02         R04         ..............................            R12         respectively

-There are 6 points whence   6  STO 00

   XEQ "MEANS"  >>>>   m  = 10.8718        ( execution time = 7 seconds )
                              RDN    m' =  9.8020
                              RDN    m" =  8.0549
 

2°)  Two Variables ( Grouped Data )
 

     a) Standard Deviation, Covariance, Coefficient of Correlation
 
 

-This program calculates the arithmetic mean, the standard deviations and the covariance
  of a given set of p points  ( x₁ , y₁ )  , ( x₂ , y₂ )  ............ ,  ( x_p , y_p )
  with respective frequencies  n₁ , n₂ , ............ ,  n_p
 

Data Registers:    R00 = Sum n_i x_i     R02 = Sum n_i y_i      R04 = Sum n_i x_i y_i    and when the routine stops,  R06 = m_x = arithmetic mean
                               R01 = Sum n_i x_i²    R03 = Sum n_i y_i²    R05 = n = Sum n_i                                                  R07 = m_y = arithmetic mean

Flag:   F01                  CF 01  if you want to use the sample standard deviation & covariance
Subroutines: /            SF 01  if you want to use the population standard deviation & covariance

-If your statistics registers are R00 thru R05, lines 01 to 22 may be replaced by

     LBL A
     X<> Z
     LBL 01
     S+                      ( Sigma+ )
     LASTX
     DSE Z
     GTO 01
     RTN
 
 

01  LBL A 02  ST+ 05 03  STO 06 04  X<>Y 05  * 06  STO 07 07  ST+ 02 08  LASTX          09  *  10  ST+ 03  11  X<>Y 12  ENTER^

13  ST* 07  14  RCL 06  15  *  16  ST+ 00 17  *  18  ST+ 01 19  RCL 07 20  ST+ 04 21  RCL 06 22  RTN 23  LBL "STAT2" 24  RCL 04

25  RCL 02 26  STO 07 27  RCL 00 28  RCL 05 29  ST/ 07 30  / 31  STO 06          32  * 33  - 34  RCL 01 35  RCL 06  36  X^2

37  RCL 05 38  * 39  - 40  RCL 07 41  X^2 42  RCL 05 43  * 44  RCL 03          45  X<>Y 46  - 47  RCL 05  48  FC? 01

49  DSE X 50  ST/ T 51  ST/ Z 52  / 53  SQRT  54  ST/ T 55  X<>Y 56  SQRT             57  ST/ T 58  END

 
    ( 83 bytes / SIZE 008 )
 
 

STACK	INPUTS#i	OUTPUTS#i	final OUTPUTS
T	/	/	r
Z	xi	/	cov(x,y)
Y	yi	/	sy
X	ni	ni	sx
L	/	/	Vx = sx2

       s = sample standard deviation if CF 01                             cov(x,y) = sample covariance if CF 01                   r = coefficient of correlation
       s = population standard deviation if SF 01                        cov(x,y) = population covariance if SF 01

      V = s² = variance       ,     R06 = m_x = arithmetic mean   ,    R07 = m_y = arithmetic mean

Example:    With the following data
 
 

xi	2	5	10	16	21
yi	4	7	12	19	24
ni	3	4	7	5	2

 
-Clear registers R00 thru R05  ( for instance execute CLS  if your statistics registers are R00 thru R05 )
-Then,  x_i  ENTER^   y_i  ENTER^    n_i   S+  ( [A]  in user mode )  for all the points

   2  ENTER^      5  ENTER^      10  ENTER^     16  ENTER^    21  ENTER^
   4  ENTER^      7  ENTER^      12  ENTER^     19  ENTER^    24  ENTER^
   3  S+                4  S+                7   S+                5   S+              2   S+            in user mode

-Finally, simply press  R/S  ( or XEQ "STAT2" )  it yields:

     if CF 01                         s_x  =  5.9622                                       if SF 01                               s_x = 5.8185
                       RDN            s_y =  6.3808                                                              RDN            s_y  =  6.2270
                       RDN   cov(x,y) = 38.0143                                                              RDN   cov(x,y) = 36.2041
                       RDN            r   =  0.9992                                                               RDN            r   =  0.9992

-And we have also

       R06 =  m_x = 10.3810
       R07 =  m_y = 12.7143
 

-To correct wrong inputs, key in  x_i  ENTER^  y_i  ENTER^   n_i  CHS  S+  ( [A] in user mode )  for the wrong ( x_i , y_i , n_i )
  or add   LBL a   CHS   just before  LBL A  and key in  x_i  ENTER^  y_i  ENTER^  n_i  S-  ( [a]  in user mode )
 

     b) Linear Regression
 
 

Data Registers:    R00 = Sum n_i x_i     R02 = Sum n_i y_i      R04 = Sum n_i x_i y_i        R06 = m_x = arithmetic mean
                               R01 = Sum n_i x_i²    R03 = Sum n_i y_i²    R05 = n = Sum n_i         R07 = m_y = arithmetic mean

Flag:   F01 is used by "STAT2"
Subroutine:  "STAT2"
 
 

01  LBL "LR+"  02  XEQ "STAT2"  03  R^  04  R^  05  LASTX  06  /  07  RCL 06  08  RCL Y  09  *  10  RCL 07  11  X<>Y  12  -  13  X<>Y  14  END

 
( 29 bytes / SIZE 008 )
 
 

STACK	INPUTS	OUTPUTS
Z	/	r
Y	/	p
X	/	m

  where   y = m x + p  is the least-squares line and  r = linear correlation coefficient

Example:    With the same data as in paragraph 2-a)

   XEQ "LR+"  gives   m = 1.0694
                       RDN   p =  1.6130
                       RDN   r  =  0.9992

Notes:

-Though "STAT2" produces different outputs if flag F01 is set or clear, "LR+" yields the same results ( with perhaps a small difference in the last decimal )
-See also "LR" & "LSL" listed in "Least-Squares Approximation for the HP-41"
 

3°)  Probability Functions
 

     a)  Chi2-Probability Function
 

-This program computes P( X² | m ) = 2 ^-m/2 [ 1/Gamma(m/2) ] §₀^{X^2}   t ^-1+m/2 e ^-t/2 dt

-It uses the series exansion:    P( X² | m ) = (X²/2)^m/2 e ^{-(X^2)/2} [ 1/Gamma(m/2) ]  [ 1 + Sum_k=1,2,...  X^2k/[ (m+2)(m+4) ..... (m+2k) ]
 

Data Registers: /
Flags: /
Subroutine:  "GAM+" or "GAM" or "1/G"   ( cf "Gamma , Digamma , Polygamma & Beta Functions for the HP-41" )
 
 

01  LBL "PXN" 02  STO N 03  X<>Y 04  STO O 05  2 06  ST+ Y 07  / 08  XEQ "GAM+" 09  RCL N 10  2

11  / 12  E^X 13  ST* Y 14  X<> L 15  RCL O            16  2 17  / 18  Y^X 19  X<>Y 20  /

21  STO M 22  1 23  RCL O            24  1 25  ENTER^  26  LBL 01 27  X<> T 28  RCL N 29  * 30  R^

31  2 32  + 33  STO T 34  / 35  STO T  36  X<>Y 37  ST+ Y 38  X#Y? 39  GTO 01          40  RCL M

41  * 42  RCL N 43  SIGN 44  X<> O            45  X<>Y 46  CLA 47  END

 
   ( 78 bytes / SIZE 000 )
 
 

STACK	INPUTS	OUTPUTS
Y	m	m
X	X2	P( X2 | m )
L	/	X2

  where m is a positive integer ( m > 0 )

Example:

    7   ENTER^
   12  XEQ "PXN"  >>>>   P (12 | 7 ) = 0.899441131  ( in 14 seconds )
 

     b)  Gaussian Probability Function
 

 "GPF"  calculates  P(x) = (1/sqrt(2.pi)) § _-infinity^x  exp (-t²/2).dt    via the error function if x >= 0 and via the complementary error function if x < 0

Formulae:

          P(x) = [ 1 + erf ( x/sqrt(2) ) ] / 2      if x >= 0
          P(x) = (1/2) erfc ( x/sqrt(2) )            if x < 0

          Q(x) = 1 - P(x)  verifies  Q(x) = P(-x)
 

Data Registers: /
Flag:  F10
Subroutine:  "ERF"  ( cf "Error Function for the HP-41" )
 
 

01  LBL "GPF"  02  CF 10  03  X<0?  04  SF 10  05  2  06  SQRT  07  /  08  XEQ "ERF"  09  FS? 10  10  X<>Y  11  .5  12  ST* Y  13  FS?C 10  14  CLX  15  +  16  END

 
( 34 bytes / SIZE 000 )
 
 

STACK	INPUT	OUTPUT
X	x	P(x)

 
Examples:

  3.14  XEQ "GPF"  >>>>  P(3.14) = 0.999155261 = Q(-3.14)                  ( in 8 seconds )
   -1         R/S          >>>>    P(-1)   = 0.158655254 =   Q(1)                       ( in 7 seconds )
 

   3-c)  F-Distribution Function
 

       3-c1)  Program#1
 

-This program calculates    Q( x | n₁ , n₂ ) =  (n₁/n₂)^(n1)/2  Gamma[(n₁+n₂)/2] / [ Gamma(n₁/2) Gamma(n₂/2) ] §_x^infinity  t ^-1+(n1)/2 ( 1 + t.n₁/n₂ ) ^-(n1+n2)/2 dt

  via the incomplete beta function:   Q( x | n₁ , n₂ ) = I_y( n₂/2 , n₁/2 )   with  y = n₂/(n₂+n₁.x)
 

Data Registers:   R00 thru  R05
Flag:   F10
Subroutines:   "BETA"  "IBETA"  ( cf  "Gamma , Digamma , Polygamma & Beta Functions for the HP-41" )
 
 

01  LBL "FDF"      02  RCL Z 03  STO 01 04  * 05  STO Z 06  RCL Y 07  STO 05

08  + 09  ST/ Z 10  / 11  CF 10             12  X>Y? 13  SF 10 14  X>Y?

15  X<>Y 16  STO 00  17  2 18  ST/ 01 19  ST/ 05 20  RCL 01          21  RCL 05

22  XEQ "BETA" 23  X<> 05 24  RCL 01 25  FS? 10 26  X<>Y 27  RCL 00 28  XEQ "IBETA"

29  FS? 10 30  CHS 31  RCL 05 32  FS?C 10          33  ST+ Y 34  / 35  END

 
   ( 67 bytes / SIZE 006 )
 
 

STACK	INPUTS	OUTPUTS
Z	n1	/
Y	n2	/
X	x >= 0	Q( x | n1 , n2 )

 
Examples:

     7    ENTER^
     6    ENTER^
   4.3   XEQ "FDF"  >>>>  Q ( 4.3 |  7 , 6 ) = 0.04764079963  ( in 24 seconds )

     7    ENTER^
     9    ENTER^
   0.4      R/S       >>>>       Q ( 0.4 | 7 , 9 ) =  0.879873501  ( in 27 seconds )

Notes:

-We have used the relation   I_y( a , b ) = 1 - I_1-y( b , a )  to reduce execution time and roundoff errors when y > 0.5
-  Q( 0 | n₁ , n₂ ) =  1
 

       3-c2)  Inverse FDF function
 

-"IFDF" employs Halley's iterative method to solve  Q( y | n₁ , n₂ ) = x   i-e   to find  y = IFDF ( x | n₁ , n₂ )

       t_k+1 = t_k - 2 f(t_k) f '(t_k) / [ 2 ( f '(t_k) )² - f(t_k) f ''(t_k) ]      where   f = FDF
 

Data Registers:   R00 thru R11: temp
Flags: /
Subroutines:   "BETA"         ( cf  "Gamma , Digamma , Polygamma & Beta Functions for the HP-41" )
                          "FDF"           ( cf paragraph above )

-Line 21,  1  is the 1st guess whch could probably be replaced by a better estimation
-Lines 59-60 are probably superfluous
 
 

01  LBL "IFDF"  02  STO 06  03  RDN  04  STO 08  05  X<>Y  06  STO 07  07  2  08  ST/ Z  09  /  10  XEQ "BETA"  11  RCL 07  12  RCL 08  13  /  14  RCL 07  15  2  16  /  17  Y^X

18  X<>Y  19  /  20  STO 09  21  1  22  STO 10  23  LBL 01  24  VIEW 10  25  RCL 07  26  RCL 08  27  RCL 10  28  XEQ "FDF"    29  RCL 06  30  -  31  STO 11  32  RCL 10  33  RCL 07  34  2

35  /  36  1  37  -  38  STO 03  39  Y^X  40  RCL 07  41  RCL 10  42  *  43  RCL 08  44  ST+ Y  45  /  46  STO 04          47  RCL 07  48  RCL 08  49  +  50  2  51  /

52  STO 05  53  Y^X  54  /  55  RCL 09  56  *  57  RCL 03  58  RCL 10  59  X=0?  60  SIGN  61  /  62  RCL 07  63  RCL 08          64  /  65  RCL 05  66  *  67  RCL 04  68  /

69  -  70  RCL 11  71  STO T  72  *  73  2  74  /  75  +  76  /  77  ST+ 10  78  ABS  79   E-4  80  X<Y?  81  GTO 01  82  RCL 10          83  CLD  84  END

 
    ( 110 bytes / SIZE 012 )
 
 

STACK	INPUTS	OUTPUTS
Z	n1	/
Y	n2	/
X	x	ifdf ( x | n1 , n2 )

  Where  y = ifdf ( x | n₁ , n₂ )  is the solution of   FDF ( y | n₁ , n₂ ) = x

Examples:

     7     ENTER^
     6     ENTER^
   0.05  XEQ "IFDF"  >>>>   IFDF ( 0.05 | 7 , 6 ) = 4.206658486                     ---Execution time = 72s---

     3     ENTER^
     5     ENTER^
  0.001 XEQ "IFDF"  >>>>  IFDF ( 0.001 | 3 , 5 ) = 33.20246321                    ---Execution time = 93s---

Notes:

-The execution times suppose that "IBETA" calls the M-Code routine HGF+
-Line 79  E-4  is a small number to control the precision of the result.
-Since Halley's method is of order 3, almost all the decimals are already correct !
 

   3-d)  Student's t-Distribution & UTPT Functions
 

       3-d1)  Program#1
 

- "STD"  computes  A( t | n ) = n ^-1/2 / Beta(0.5 ; n/2 )  §_-t^+t  ( 1 + x²/n ) ^-(n+1)/2  dx

-This routine also computes   n ^-1/2 / Beta(0.5 ; n/2 )  §_t^+infinity  ( 1 + x²/n ) ^-(n+1)/2  dx = [ 1 - A ( t | n ) ] / 2 = UTPT(t,n)  on the HP-48

Formula:   A( t | n ) = 1 - I_y( n/2 , 1/2 )   with   y = n/(n+t²)
 

Data Registers:   R00 thru  R05
Flag:  F10
Subroutines:   "GAM+" ( or "GAM" )  "IBETA"  ( cf  "Gamma , Digamma , Polygamma & Beta Functions for the HP-41" )
 
 

01  LBL "STD"     02  X^2 03  STO Z 04  RCL Y 05  STO 00 06  + 07  ST/ Z 08  / 09  CF 10 10  X>Y?

11  SF 10 12  X>Y? 13  X<>Y 14  X<> 00 15  .5 16  * 17  STO 05 18  LASTX 19  + 20  XEQ "GAM+"

21  STO 04 22  RCL 05 23  XEQ "GAM+" 24  PI 25  SQRT 26  * 27  RCL 04 28  / 29  X<> 05 30  .5

31  FS? 10 32  X<>Y 33  RCL 00 34  XEQ "IBETA" 35  ENTER^ 36  CHS 37  RCL 05 38  + 39  FC?C 10 40  X<>Y

41  LASTX           42  ST/ Z 43  ST+ X 44  / 45  X<>Y 46  END

 
   ( 83 bytes / SIZE 006 )
 
 

STACK	INPUTS	OUTPUTS
Y	n	UTPT(t,n)
X	t >= 0	A( t | n )

 
Examples:

   4   ENTER^
   1   XEQ "STD"   >>>>            A( 1 | 4 ) = 0.626099033        ( in 21 seconds )
                              RDN     UTPT( 1 , 4 ) = 0.186950483

  4   ENTER^
 10    R/S         >>>>            A( 10 | 4 ) = 0.999437997          ( in 16 seconds )
                        RDN     UTPT( 10 , 4 ) = 0.0002810018113

-We have used the symmetry relation   I_y( a , b ) = 1 - I_1-y( b , a )  to reduce execution time and roundoff errors when y > 0.5
- A( 0 | n ) = 0
 

       3-d2)  UTPT function
 

-There is however a possible loss of significant digits with "STD" if y > 0.5
-To overcome that difficulty, the following program employs a quadratic transformation to calculate the hypergeometric function F. It yields:

   UTPT (t | n)  = { 1/[ n.B(n/2,1/2) ] }  y^n/2  F ( 1 , n ; 1+n/2 ; y / [ 2 + 2 (1-y)^1/2 ] )     where  B = beta function  and  y = n/(n+t²)
 

Data Registers:   R00 thru  R06
Flags:  /
Subroutines:   "1/G+" ( or "GAM+" )    ( cf  "Gamma , Digamma , Polygamma & Beta Functions for the HP-41" )
                          "HGF"                                           ( cf  "Hypergeometric functions for the HP-41" )
 
 

01  LBL "UTPT"  02  STO 04  03  X<>Y  04  STO 02  05  2  06  /  07  STO 03  08  XEQ "1/G+"  09  STO 06  10  RCL 03

11  .5  12  +  13  XEQ "1/G+"  14  PI  15  SQRT  16  *  17  RCL 06  18  /  19  STO 06  20  1

21  STO 01  22  ST+ 03  23  RCL 02   24  RCL 02  25  RCL 04         26  X^2  27  STO T  28  +  29  ST/ Z  30  /

31  STO 05  32  X<>Y  33  SQRT  34  1  35  +  36  ST+ X  37  /  38  XEQ "HGF"  39  RCL 05  40  RCL 02

41  2  42  /  43  Y^X  44  *  45  RCL 02   46  RCL 06         47  *  48  /  49  END

 
    ( 77 bytes / SIZE 007 )
 
 

STACK	INPUTS	OUTPUTS
Y	n	/
X	t >= 0	UTPT(t,n)

 
Examples:

   •  4    ENTER^
     2.8  XEQ "UTPT"  >>>>   UTPT ( 2.8 | 4 ) = 0.02440577530           ---Execution time = 10s---

   •  7    ENTER^
     10     R/S        >>>>   UTPT ( 10 | 4 ) = 1.069710145 E-5            ---Execution time = 8s---

Note:

-If you to compute  A(t | n) , simply add   ENTER^   ST+ X  1  -  CHS  X<>Y   after line 48.
 

       3-d3)  Inverse UTPT function
 

-IUTPT uses the same expression to calculte UTPT(t,n)  and find a solution of the equation  UTPT(t,n) - A = 0  where  A and n are given
-It employs Halley's method which is even better than Newton's method:

       t_k+1 = t_k - 2 f(t_k) f '(t_k) / [ 2 ( f '(t_k) )² - f(t_k) f ''(t_k) ]     with   f =  UTPT
 

Data Registers:   R00 thru  R09
Flags:  /
Subroutines:   "1/G+" ( or "GAM+" )    ( cf  "Gamma , Digamma , Polygamma & Beta Functions for the HP-41" )
                          "HGF"                                           ( cf  "Hypergeometric functions for the HP-41" )
 
 
 

01  LBL "IUTPT"  02  STO 06  03  RDN  04  STO 02  05  1  06  STO 01  07  STO 03  08  +  09  2  10  /  11  STO 08  12  X<>Y  13  STO 04  14  RCL 02  15  2  16  /  17  ST+ 03

18  STO 05  19  XEQ "1/G+"  20  STO 07  21  RCL 08  22  XEQ "1/G+"   23  PI  24  SQRT  25  *  26   RCL 07  27  /  28  ST* 04  29  LBL 01  30  VIEW 06  31  RCL 02  32  RCL 02  33  RCL 06  34  X^2

35  STO T  36  +  37  STO 09  38  ST/ Z  39  /  40  STO 07  41  X<>Y  42  SQRT  43  1  44  +  45  ST+ X  46  /  47  XEQ "HGF"   48  RCL 07  49  RCL 05  50  Y^X  51  *

52  RCL 02  53  /  54  RCL 04  55  -  56  RCL 07  57  RCL 08  58  Y^X  59  RCL 02  60  SQRT  61  /  62  RCL 02   63  RCL 06  64  *  65  RCL 08          66  *  67  RCL 09  68  X^2

69  /  70  R^  71  *  72  RCL 07   73  /  74  -  75  /  76  ST+ 06  77  ABS  78   E-8  79  X<Y?  80  GTO 01  81  RCL 06          82  CLD  83  END

 
     ( 117 bytes / SIZE 010 )
 
 

STACK	INPUTS	OUTPUTS
Z	A
Y	n	/
X	t0	IUTPT(A,n)

  where  t₀  is an initial guess  and  t = IUTPT(A,n) is  the solution of the equation  UTPT(t,n) = A

Example1:          A = 0.0001 ,  n = 4  ,  with the initial approximation    t₀ = 10

     0.0001   ENTER^
         4        ENTER^
        10       XEQ "IUTPT"  >>>>   t = 13.03367172                   ---Execution time = 22s---
 

Example2:          A = 0.025 ,  n = 21  ,  with the initial approximation    t₀ = 1

       0.025  ENTER^
         21     ENTER^
          1      XEQ "IUTPT"   >>>>   t = 2.079613847                   ---Execution time = 63s---
 

Example3:          A = 0.001 ,  n = 21  ,  with the initial approximation    t₀ = 1

       0.001  ENTER^
         21     ENTER^
          1      XEQ "IUTPT"   >>>>   t = 3.527153671                   ---Execution time = 66s---
 

   3-e)  Poisson Distribution
 

-The function is defined by   P(X=k) = exp(-L) L^k / k!        where  L > 0  and  k is a non-negative integer
 

Data Registers: /
Flags: /
Subroutines: /
 
 

01  LBL "POI"  02  RCL Y  03  X<>Y  04  Y^X  05  LASTX  06  FACT  07  /  08  X<>Y  09  E^X  10  /  11  END

 
( 20 bytes / SIZE 000 )
 
 

STACK	INPUTS	OUTPUTS
Y	L	/
X	k	P( X=k )

 
Example:           L = 5  ,  k = 6

   5   ENTER^
   6   XEQ "POI"   >>>>   P(X=6) = 0.146222808

-The cumulative distribution is given by   P( X <= k ) = exp(-L) SUM_j=0,1,...,k  L^j / j!
 

Data Register:  R00 = k , k-1 , ..... , 0
Flags: /
Subroutines: /
 
 

01  LBL "POI+"  02  X=0?  03  GTO 02  04  STO 00  05  CLX  06  LBL 01  07  RCL Y  08  RCL 00  09  Y^X  10  LASTX  11  FACT  12  /  13  +  14  DSE 00  15  GTO 01  16  LBL 02  17  1  18  +  19  RCL Y  20  E^X  21  /  22  END

 
( 36 bytes / SIZE 001 )
 
 

STACK	INPUTS	OUTPUTS
Y	L	/
X	k	P( X<=k )

 
Example:           L = 5  ,  k = 6

   5   ENTER^
   6   XEQ "POI+"   >>>>   P(X<=6) = 0.762183463
 

   3-f)  Binomial Distribution
 

-This routine computes  P(X=k) = C_n^k p^k (1-p)^n-k    where   0 <= k <= n  are integers and  0 < p < 1
 

Data Registers:  R00 thru R02: temp
Flags: /
Subroutine:  "CNK"  or an M-Code CNK  cf §4-a) below
 
 

01  LBL "BNP"  02  STO 00  03  X<>Y  04  STO 01  05  X<>Y  06  Y^X  07  X<>Y  08  STO 02  09  RCL 00  10  XEQ "CNK"  11  *  12  1  13  RCL 01  14  -  15  RCL 02  16  RCL 00  17  -  18  Y^X  19  *  20  END

 
( 32 bytes / SIZE 003 )
 
 

STACK	INPUTS	OUTPUTS
Z	n	/
Y	p	/
X	k	P( X=k )

Example:        n = 100 ,  p = 1/6

   100   ENTER^
     6     1/X
    20    XEQ "BNP"   >>>>   P(X=20) = 0.067861995
 

-The cumulative distribution could be calculated by calling "BNP"  (k+1) times.
-The following routine is faster, at least for large k-values:
 

Data Registers:  R00 thru R03: temp
Flags: /
Subroutines: /
 
 

01  LBL "BNP+" 02  STO 00 03  SIGN 04  X<>Y 05  STO 01 06  - 07  ST/ 01 08  X<>Y

09  STO 02 10  Y^X 11  SIGN 12  RCL 00         13  X=0? 14  GTO 02  15  CLX 16  STO 03

17  LASTX 18  LBL 01 19  LASTX  20  RCL 01         21  * 22  RCL 02 23  RCL 03  24  -

25  * 26  RCL 03  27  1 28  + 29  STO 03         30  / 31  + 32  DSE 00

33  GTO 01         34  SIGN 35  LBL 02  36  LASTX 37  END

 
   ( 50 bytes / SIZE 004 )
 
 

STACK	INPUTS	OUTPUTS
Z	n	/
Y	p	/
X	k	P( X<=k )

 
Example:        n = 100 ,  p = 1/6

   100   ENTER^
     6     1/X
    15    XEQ "BNP+"   >>>>   P(X<=15) = 0.387657550
 

   3-g)  Negative Binomial Distribution
 

-Let  r a positive integer, k a non-negative integer and a real p  0 < p < 1

     P(X=k) = C_k+r-1^r-1 p^r (1-p)^k
 

Data Registers:  R00-R01: temp
Flags: /
Subroutine:  "CNK"  or an M-Code CNK  cf §4-a) below
 
 

01  LBL "JRP"  02  STO 00  03  SIGN  04  X<>Y  05  STO 01  06  -  07  RCL 00  08  Y^X  09  RCL 01  10  RCL Z  11  Y^X  12  *  13  RCL 00  14  R^  15  1  16  -  17  ST+ Y  18  XEQ "CNK"  19  *  20  END

 
( 34 bytes / SIZE 002 )
 
 

STACK	INPUTS	OUTPUTS
Z	r	/
Y	p	/
X	k	P( X=k )

 
Example:       r = 5  ,  p = 1/7

    5   ENTER^
    7   1/X
   10  XEQ "JRP"  >>>>   P(X=10) = 0.012748996
 

 "JRP+"  computes the cumulative negative binomial distribution
 

Data Registers:  R00 thru R03 temp
Flags: /
Subroutines: /
 
 

01  LBL "JRP+" 02  STO 00       03  RDN 04  STO 01 05  X<>Y 06  STO 02 07  Y^X

08  ABS 09  RCL 00       10  X=0? 11  GTO 02 12  CLX 13  STO 03 14  X<>Y

15  LBL 01 16  LASTX 17  RCL 02  18  RCL 03       19  + 20  * 21  1

22  ST+ 03 23  RCL 01  24  - 25  * 26  RCL 03       27  / 28  +

29  DSE 00 30  GTO 01       31  SIGN 32  LBL 02  33  LASTX 34  END

 
   ( 47 bytes / SIZE 004 )
 
 

STACK	INPUTS	OUTPUTS
Z	r	/
Y	p	/
X	k	P( X<=k )

 
Example:       r = 5  ,  p = 1/7

    5   ENTER^
    7   1/X
   10  XEQ "JRP+"  >>>>   P(X<=10) = 0.051581690
 

   3-h)  Pascal Distribution
 

-The function is defined as   P(X=k) = C_k-1^r-1 p^r (1-p)^k-1    where   r and k are positive integers ( r <= k )  and  0 < p < 1
 

Data Registers:  R00-R01-R02: temp
Flags: /
Subroutine:  "CNK"  or an M-Code CNK  cf §4-a) below
 
 

01  LBL "PRP"   02  STO 00 03  1 04  - 05  RDN

06  STO 01 07  RDN 08  STO 02       09  SIGN 10  -

11  XEQ "CNK" 12  RCL 01 13  RCL 02 14  Y^X 15  *

16  1 17  RCL 01        18  - 19  RCL 00 20  RCL 02

21  - 22  Y^X              23  * 24  END

 
   ( 36 bytes / SIZE 003 )
 
 

STACK	INPUTS	OUTPUTS
Z	r	/
Y	p	/
X	k	P( X=k )

 
Example:       r = 10  ,  p = 0.6

  10   ENTER^
  0.6  ENTER^
  15   XEQ "PRP"  >>>>  P(X=15) = 0.123958563
 

-Now,  "PRP+"  calculates  P(X<=k)
 

Data Registers:  R00 thru R03: temp
Flags: /
Subroutines: /
 
 

01  LBL "PRP+" 02  RDN 03  STO 01       04  X<>Y 05  STO 02 06  STO 03 07  R^ 08  X<Y?

09  GTO 02  10  X<>Y 11  - 12  STO 00       13  RCL 01 14  RCL 02 15  Y^X 16  ABS

17  ISG 00 18  LBL 01 19  DSE 00  20  X<0? 21  RTN 22  LASTX 23  1 24  RCL 01

25  - 26  * 27  RCL 03  28  ST* Y 29  1 30  + 31  STO 03  32  RCL 02

33  - 34  / 35  + 36  GTO 01       37  LBL 02 38  CLX 39  END

 
   ( 53 bytes / SIZE 004 )
 
 

STACK	INPUTS	OUTPUTS
Z	r	/
Y	p	/
X	k	P( X<=k )

 
Example:       r = 10  ,  p = 0.6

  10   ENTER^
  0.6  ENTER^
  15   XEQ "PRP"  >>>>  P(X<=15) = 0.403215551
 

Note:      P(X=k)  is in L-register
 

   3-i)  Geometric Distribution
 

-This probability function is defined by                                    P(X=k)  =  p (1-p)^k-1                  where   0 < p < 1  and  k is a positive integer
-The cumulative distribution may easily be computed since     P(X<=k) = 1 - (1-p)^k

-So a simple program can return both results.
 

Data Registers: /
Flags: /
Subroutines: /
 
 

01  LBL "GP"  02  1  03  RCL Z  04  -  05  ENTER^  06  X<> Z  07  DSE X  08  INT  09  Y^X  10  ST* Z  11  *  12  1  13  X<>Y  14  -  15  X<>Y  16  END

 
( 27 bytes / SIZE 000 )
 
 

STACK	INPUTS	OUTPUTS
Y	p	P( X<=k )
X	k	P( X=k )

 
Example:    p = 0.6

    0.6    ENTER^
     4      XEQ "GP"   >>>>    P(X=4)  = 0.0384
                                RDN    P(X<=4) = 0.9744
 

   3-j)  Hypergeometric Distribution
 

-Here,   P(X=k)  =  C_m^k C_n-m^a-k / C_n^a    where   n , m , a , k  are integers with  m <= n ,   k <= a <= m  ,  a - k <= n - m    ( cf § 4-c) for a concrete example )
 

Data Register:  R00: temp
Flags: /
Subroutine:    "CNK"  ( cf 4-a) )
 
 

01  LBL "HGD"  02  RDN  03  STO 00  04  X<>Y  05  R^  06  ST- 00  07  XEQ "CNK"  08  RDN  09  XEQ "CNK"  10  ST/ Z  11  X<> T  12  X<>Y  13  -  14  RCL 00  15  XEQ "CNK"  16  *  17  END

 
( 40 bytes / SIZE 001 )
 
 

STACK	INPUTS	OUTPUTS
T	n	/
Z	m	/
Y	a	/
X	k	P(X=k)

 
Example:     n = 30  ,  m = 12  ,  a = 10

  30  ENTER^
  12  ENTER^
  10  ENTER^
   5   XEQ "HGD"  >>>>   P(X=5) = 0.225856303
 

- "HGD+" below computes  P(X<=k)
 

Data Registers:  R00 thru R05: temp
Flags: /
Subroutine:  "CNK"  or an M-Code CNK  cf §4-a)
 
 

01  LBL 'HGD+" 02  STO 00 03  RDN 04  STO 01 05  RDN 06  STO 02 07  X<>Y 08  STO 03 09  - 10  RCL 01 11  + 12  STO 05 13  X<0?

14  CLX 15  STO 04 16  RCL 00 17  - 18  X>0? 19  GTO 02 20  CHS 21  STO 00 22  RCL 02 23  RCL 04 24  XEQ "CNK" 25  RCL 03 26  RCL 02

27  - 28  RCL 01 29  RCL 04 30  - 31  XEQ "CNK" 32  * 33  RCL 03 34  RCL 01 35  XEQ "CNK" 36  / 37  ABS 38  ISG 00 39  LBL 01

40  DSE 00 41  X<0? 42  RTN 43  LASTX 44  RCL 01         45  RCL 04  46  - 47  * 48  RCL 02 49  RCL 04 50  - 51  * 52  RCL 04

53  1 54  + 55  STO 04 56  ST/ Y 57  RCL 05         58  - 59  / 60  + 61  GTO 01 62  LBL 02 63  CLX 64  END

 
   ( 90 bytes / SIZE 006 )
 
 

STACK	INPUTS	OUTPUTS
T	n	/
Z	m	/
Y	a	/
X	k	P(X<=k)

 
Example1:       n = 30  ,  m = 12  ,  a = 10

  30  ENTER^
  12  ENTER^
  10  ENTER^
   5   XEQ "HGD+"  >>>>   P(X<=5) = 0.881728367

Example2:       n = 100  ,  m = 84  ,  a = 41

  100   ENTER^
   84    ENTER^
   41    ENTER^
   37    XEQ "HGD+"   >>>>  P(X<=37) =  0.958543068
 

4°)  Classic Problems
 

     4-a)  Binomial Coefficients
 

-This routine is only a slight modification of a program by Keith Jarret listed in "Synthetic Programming Made Easy"
-I've simply added lines 09-10-20 to avoid a data error if k = 0 or k = n

-"CNK" computes  C_n^k =  n! / [ k! (n-k)! ]
 

Data Registers: /
Flags: /
Subroutines: /
 
 

01  LBL "CNK" 02  ST- Y 03  SIGN 04  STO M 05  CLX

06  LASTX  07  X>Y? 08  X<>Y 09  X=0? 10  GTO 02

11  LBL 01        12  X<>Y 13  ISG X 14  CLX 15  ST* M

16  X<>Y 17  ST/ M 18  DSE X         19  GTO 01  20  LBL 02

21  X<> M        22  X<>Y 23  RDN 24  END

 
   ( 41 bytes / SIZE 000 )
 
 

STACK	INPUTS	OUTPUTS
T	T	n
Z	Z	T
Y	n	Z
X	k	Cnk
L	/	k

 
Example:      n = 100  ,  k = 84

    100  ENTER^
     84   XEQ "CNK"  >>>>   C₁₀₀⁸⁴ = 1.345860628 10¹⁸    ( in 5 seconds )

-One could of course use the buil-in FACT function, but there will be an OUT OF RANGE if n > 69

-The following M-Code routine performs the same calculations.
-Moreover, 0 is returned if n or k or n-k is negative or fractional.
-However, it does not check for alpha data or overflow.
-Synthetic register Q is used.
 

08B   "K"
00E   "N"
003   "C"
2A0   SETDEC
3C8   CLRKEY
0F8   READ 3(X)
128   WRIT 4(L)
0B8   READ 2(Y)
0E8   WRIT 3(X)
04E   C=0 ALL
0A8  WRIT 2(Y)
0F8   READ 3(X)
2FE  ?C#0 MS
0BF  JC+23d=+17h
084   CLRF 5                 C
0ED  ?NCXQ                 =
064   193B                   frc(C)
2EE   ?C#0 ALL
097   JC+18d=+12h
138   READ 4(L)
2FE  ?C#0 MS
07F   JC+15d=+0Fh
084   CLRF 5                 C
0ED  ?NCXQ                 =
064   193B                   frc(C)
2EE   ?C#0 ALL
057   JC+10d=+0Ah
138   READ 4(L)
2BE  C=-C-1 MS         C=-C
070   N=C ALL
10E   A=C ALL
0F8   READ 3(X)
01D  ?NCXQ               C=
060   1807                   A+C
2FE  ?C#0 MS
01B   JNC+03
3A5   ?NCGO          ?NCgoto
052    14E9                  RDN
10E    A=C ALL
0F0    C<>N ALL
01D  ?NCXQ               C=
060   1807                   A+C
2FE  ?C#0 MS
01F   JC+03
138   READ 4(L)
070   N=C ALL
04E   C=0 ALL
268   WRIT 9(Q)
35C  PT=12
050   LD@PT- 1
0A8  WRIT 2(Y)
278   READ 9(Q)
10E   A=C ALL
0B0   C=N ALL
36E   ?A#C ALL
36B   JNC-19d=-13h
1BE   A=A-1 MS
0F8   READ 3(X)
01D  ?NCXQ               C=
060   1807                   A+C
10E   A=C ALL
0B8   READ 2(Y)
135   ?NCXQ               C=
060   184D                  A*C
0A8  WRIT 2(Y)
278   READ 9(Q)
00E   A=0 ALL             A
35C  PT=12                  =
162   A=A+1@PT        1
01D  ?NCXQ               C=
060   1807                   A+C
268   WRIT 9(Q)
10E   A=C ALL
0B8   READ 2(Y)
0AE  A<>C ALL
261   ?NCXQ               C=
060   1898                   A/C
3CC  ?KEY
360    ?C RTN                                        the routine stops here if you press a key ( useful if, for example,  n = 10000000 and k = 1000000 )
31B    JNC-29d=-1Dh
 
 

STACK	INPUTS	OUTPUTS
T	T	n
Z	Z	T
Y	n	Z
X	k	Cnk
L	/	k

 
Example:      n = 400  ,  k = 100

    400  ENTER^
    100   XEQ "CNK"  >>>>   C₄₀₀¹⁰⁰ = 2.241854802 10⁹⁶    ( in 12 seconds )

-The last decimals should be 792 instead of 802

-With, for instance, n = 4 and k = 0.5  "CNK" returns 0.
-This may be useful when "CNK" is used in the cumulative distributions above.
 

     4-b)  Probability of no repetition in a sample
 

-The typical exercise is:  "In a room containing m persons, what is the probability P that nobody has the same birthday?"  ( neglecting the leap years )

-The answer is  P = ( 1 - 1/n ).( 1 - 2/n )........( 1 - (m-1)/n )           with  n = 365
 

Data Registers: /
Flags: /
Subroutines: /
 
 

01  LBL "NRS"  02  X<>Y  03  ENTER^  04  SIGN  05  LBL 01  06  LASTX  07  *  08  X<>Y  09  ST/ Y  10  X<>Y  11  DSE L  12  ""   13  DSE Z  14  GTO 01  15  END

 
( 27 bytes / SIZE 000 )
 
 

STACK	INPUTS	OUTPUTS
Y	n	n
X	m	probability

 
Example:    n = 365 ,  m = 24

   365   ENTER^
    24    XEQ "NRS"  >>>>  probability = 0.4617   ( in 8 seconds )
 

     4-c)  Hypergeometric Distribution ( Urn Problem )
 

 "An urn contains m black marbles and (n-m) white marbles. You take simultaneously ( or successively without replacement )  "a"  marbles.
  What is the probability P that you have drawn exactly k black marbles ( and (a-k) white ones ) ?"

-The answer is:    P = C_m^k C_n-m^a-k / C_n^a

-So we can use the program "HGD" listed in § 3-j above.
 

Example:    The urn contains 60 marbles ( 15 black ones and 45 white ones )  You take 10 marbles. Probability that you have exactly 3 black marbles?

  n = 60 , m = 15 , a = 10 , k = 3

  60  ENTER^
  15  ENTER^
  10  ENTER^
   3   XEQ "HGD"  >>>>   P = 0.273864227
 

     4-d)  Multivariate Hypergeometric Distribution
 

-The previous exercise may be generalized if the urn contains    n₁ marbles of color 1  ,   n₂ marbles of color 2  ,   .............................  ,  n_k marbles of color k
-You take simultaneously - or successively without replacement -  m marbles.
-What is the probability P that you have drawn exactly    p₁ marbles of color 1  ,   p₂ marbles of color 2  ,   .............................  ,  p_k marbles of color k
  ( of course  m =  p₁ +  p₂ + ............  +  p_k )

-The program below computes the answer   P = C_n1^p1 C_n2^p2 ......  C_nk^pk   / C_{n1+n2+....+nk}^{p1+p2+....+pk}
 

Data Registers:       •  R00 = k                                     ( Registers R00 thru R2k  are to be initialized before executing "MHGD" )

                                  •  R01 = n₁     •  R03 = n₂      .......    •  R2k-1 = n_k
                                  •  R02 = p₁     •  R04 = p₂      .......    •  R2k   = p_k
Flags: /
Subroutine:   "CNK"  ( cf 4-a) )
 
 

01  LBL "MHGD" 02  RCL 00 03  ST+ X 04  STO N 05  1 06  LBL 01

07  RCL IND Y 08  DSE Z 09  RCL IND Z 10  X<>Y 11  XEQ "CNK"  12  *

13  DSE Y 14  GTO 01 15  0 16  ENTER^ 17  LBL 02 18  RCL IND N

19  + 20  DSE N 21  X<>Y 22  RCL IND N  23  + 24  X<>Y

25  DSE N 26  GTO 02 27  XEQ "CNK"  28  / 29  END

 
   ( 58 bytes / SIZE 2k+1 )
 
 

STACK	INPUT	OUTPUT
X	/	probability

 
Example:        An urn contains 12 marbles ( 3 blue ones , 5 red ones , 4 green ones )  You take simultaneously 7 marbles
                          Calculate the probability P that you have exactly  2 blue marbles , 3 red marbles and 2 green marbles ?

  k = 3     n₁ = 3     n₂ =  5     n₃ = 4     are to be stored into      R00     R01     R03     R05
               p₁ = 2     p₂ = 3      p₃ = 2     are to be stored into                  R02     R04     R06      respectively

    XEQ "MHGD"  >>>>  probability = 0.227272727 = 5/22    ( in 6 seconds )
 

     4-c)  Permutations with repetitions
 

-This program calculates  P = n! / ( n₁! n₂! ....... n_k! )   where  n = n₁ + n₂ + ...... + n_k
 

Data Registers:     •  R00 = k                                               ( Registers R00 thru Rkk are to be initialized before executing "PERM" )

                                •  R01 = n₁   •  R02 = n₂   .....................   •  Rkk = n_k
Flags: /
Subroutines: /
 
 

01  LBL "PERM" 02  RCL 00 03  1 04  - 05  RCL IND 00

06  LBL 01 07  RCL IND Y  08  X=0? 09  GTO 03 10  LBL 02

11  X<>Y 12  ISG X           13  CLX 14  ST* L 15  X<>Y

16  ST/ L 17  DSE X  18  GTO 02         19  LBL 03  20  RDN

21  DSE Y  22  GTO 01          23  LASTX 24  END

 
   ( 43 bytes / SIZE k+1 )
 
 

STACK	INPUTS	OUTPUTS
Y	/	n
X	/	P

 
Example1:     Calculate  157! / ( 12! 20! 41! 84! )

   ( k = 4 )      n₁ = 12  ,  n₂ = 20  ,  n₃ = 41  ,  n₄ = 84   ( n = 12 + 20 + 41 + 84 = 157 )

   4  STO 00   12  STO 01   20  STO 02   41  STO 03   84  STO 04

   XEQ "PERM"  >>>>  P = 9.078363541 10⁷⁴  ( in 21 seconds )
 

Example2:     How many anagrams can you write with  the 16-letter "word"   AAABBCCCCEEEEELL

  ( k = 5 )     n₁ = 3 A's  ,  n₂ = 2 B's  ,  n₃ = 4 C's  ,  n₄ = 5 E's  ,  n₅ = 2 L's    (  n = 16 )

  5  STO 00      3  STO 01     2  STO 02      4  STO 03     5  STO 04     2  STO 05

  XEQ "PERM"   yields   302,702,400  words ( including AAABBCCCCEEEEELL itself )  So, this "word" has 302,702,399  ( meaningless ) anagrams
 

Hint:   Store the largest n_i into the last register, this will reduce execution time.
 

References:

[1]  Abramowitz and Stegun - "Handbook of Mathematical Functions" -  Dover Publications -  ISBN  0-486-61272-4
[2]  Press , Vetterling , Teukolsky , Flannery - "Numerical Recipes in C++" - Cambridge University Press - ISBN  0-521-75033-4