Search code examples
juliaprobability

How do you do arbitrary precision calculations in Julia?


I come from a Python background and an important point using Python’s Decimal module for arbitrary precision calculations is knowing when you do and don’t have to specify numbers of type Decimal in, for example, probability calculations.

How often do you have to specify arbitrary precision types for numbers in a Julia calculation?

I learn best by examples and the example I want to solve follows.

Calculate the probability of at least one success when randomly sampling without replacement a bucket of 26 ABC to Z blocks 26!, 3(26)!, and 5(26)! times.

Given:

Success on one trial = randomly drawing the blocks out in the correct order: A, B, C, to Z.

The number of random trials is n = 26! or 3n or 5n.

The probability of success on one random trial, p = 1/n

The probability of failure on one random trial is f = 1 – p.

Calculate (using type specification as little as necessary):

The confidence level, CL, of at least one success in n trials using CL = 1 – f^n

The CL for 3n trials using CL = 1 - f^(3n)

The CL for 5n trials using CL = 1 - f^(5n)


Solution

  • Julia works very well at flowing the "big" number type specification through the calculations. For example, in the following code I only specify one number as big"26", everything else is automatic.

    I have just started using Julia, but have been doing arbitrary precision calculations in various ways for years. Julia provides, hands-down, the most pleasant experience with this kind of thing I have ever had.

    # Set precision to 150 bits which should be adequate precision.
    setprecision(150)
    
    # Note that we only have to specify a "big" number here.
    n = factorial(big"26")
    println("n = factorial(big\"26\") = ", n)
    println("Note that we never have to use \"big\" again in the following code.")
    println("typeof(n) = ", typeof(n), "\n")  
    
    # p is the probability of success on 1 trial.
    p = 1/n
    println("p = 1/n = ", p)
    # Note we did not have to specify the type of p.
    println("typeof(p) = ", typeof(p), "\n")
    
    # f is the probability of failure on 1 trial.
    f = 1 - p
    println("f = 1 - p = ", f)
    println("typeof(f) = ", typeof(f), "\n")   
    
    # CL is the probability of at least 1 success in n trials.
    # CL stands for confidence level.   
    CL = 1 - f^n
    println("The 63% CL for n trials = 1 - f^n = ", CL)
    println("typeof(CL) = ", typeof(CL), "\n")   
    
    # Here is the 95% conf. level using 3n random trials.
    CL95 = 1 - f^(3n)
    println("The 95% CL for 3n trials = ", CL95)
    println("typeof(CL95) = ", typeof(CL95), "\n")
    
    # Here is the 99% conf. level using 5n random trials.
    CL99 = 1 - f^(5n)
    println("The 99% CL for 5n trials = ", CL99)
    println("typeof(CL99) = ", typeof(CL99), "\n")
    
    """ ============================= Output ==============================
    n = factorial(big"26") = 403291461126605635584000000
    Note that we never have to use "big" again in the following code.
    typeof(n) = BigInt
    
    p = 1/n = 2.4795962632247974600749435458479566174226555415e-27
    typeof(p) = BigFloat
    
    f = 1 - p = 9.9999999999999999999999999752040373677520254001e-01
    typeof(f) = BigFloat
    
    The 63% CL for n trials = 1 - f^n = 6.3212055882855767839219205578358958187929158048e-01
    typeof(CL) = BigFloat
    
    The 95% CL for 3n trials = 9.5021293163213605701567013782477488392169554992e-01
    typeof(CL95) = BigFloat
    
    The 99% CL for 5n trials = 9.9326205300091453290223898909666750856240017783e-01
    typeof(CL99) = BigFloat
    """