Probability Distributions

Summary: this tutorial shows how to use the various types of probability distributions in FSharp.Stats.

Table of contents

• Continuous
  • Normal distribution
  • Multivariate normal distribution
  • F distribution
• Discrete
  • Bernoulli distribution
  • Binomial distribution
  • Multinomial distribution
  • Hypergerometric distribution
  • Poisson distribution
  • Gamma distribution
  • Negative binomial distribution
• Empirical
• Density estimation
• Distance

FSharp.Stats provides a wide range of probability distributions. Given the distribution parameters they can be used to investigate their statistical properties or to sample non-uniform random numbers.

For every distribution the probability density function (PDF) and cumulative probability function (CDF) can be accessed. By using the PDF you can access the probability for exactly X=k success states. The CDF is used when the cumulative probabilities of X<=k is required.

Continuous

Normal distribution

The normal or Gaussian distribution is a very common continuous probability distribution. Due to the central limit theorem, randomly sampled means of a random and independent distribution tend to approximate a normal distribution It describes the probability, that under a given mean and a given dispersion (standard deviation) an event occurs exactly k times.

It is defined by two parameters N(mu,tau):

• mu = mean
• tau = standard deviation

Example: The distribution of bread weights of a local manufacturer follows a normal distribution with mean 500 g and a standard deviation of 20 g.

NormA: What is the probability of bread weights to be lower than 470 g?

NormB: What is the probability of bread weights to be higher than 505 g?

NormC: Sample independently 10 values from the normal distribution and calculate their mean.

open FSharp.Stats
open FSharp.Stats.Distributions

// Creates a normal distribution with µ = 500 and tau = 20 
let normal = Continuous.Normal.Init 500. 20.

// NormA: What is the probability of bread weights to be equal or lower than 470 g?
let normA = normal.CDF 470.
// Output: 0.06681 = 6.68 %

// NormB: What is the probability of bread weights to be higher than 505 g?
let normB = 1. - (normal.CDF 505.)
// Output: 0.401294 = 40.13 %

// NormC: Sample independently 10 values from the normal distribution and calculate their mean.
let normC = 
    Seq.init 10 (fun _ -> normal.Sample())
    |> Seq.mean

500.2734598

// Set a seed so that sampling is reproducible
let seed = 1

Random.SetSampleGenerator(Random.RandThreadSafe(seed))   
List.init 3 (fun _ -> normal.Sample())

[494.0130671; 496.2549818; 538.0814408]

Random.SetSampleGenerator(Random.RandThreadSafe(seed))   
List.init 3 (fun _ -> normal.Sample())

[494.0130671; 496.2549818; 538.0814408]

// Get back to unseeded sampling
Random.SetSampleGenerator(Random.RandThreadSafe())   
List.init 3 (fun _ -> normal.Sample())

[486.9647716; 499.1231583; 467.8681566]

open Plotly.NET

let plotNormal =
    [400. .. 600.]
    |> List.map (fun x -> x,normal.PDF x)
    |> Chart.Area
    |> Chart.withTemplate ChartTemplates.lightMirrored
    |> Chart.withTitle "N(500,20) PDF"

let plotNormalCDF =

    [400. .. 600.]
    |> List.map (fun x -> x,normal.CDF x)
    |> Chart.Area
    |> Chart.withTemplate ChartTemplates.lightMirrored
    |> Chart.withTitle "N(500,20) CDF"

Multivariate normal distribution

Multivariate normal distributions are initialized with a mean vector and a covariance matrix.

let mvn = Continuous.MultivariateNormal.Init (vector [-1.;5.]) (matrix [[0.5;1.];[0.25;1.2]])
let axisXRange = [-5. .. 0.2 .. 5.]
let axisYRange = [ 0. .. 0.2 .. 10.]

// probability density function 
let mvnPdfs =
    axisYRange |> List.map (fun y -> 
        axisXRange
        |> List.map (fun x -> 
            mvn.PDF (vector [x;y])
            )
        )

let mvnSamples = 
    Array.init 1000 (fun _ -> mvn.Sample())

let surface = Chart.Surface(mvnPdfs,axisXRange,axisYRange)

let samples = 
    mvnSamples
    |> Array.map (fun t -> t.[0],t.[1])
    |> Array.unzip
    |> fun (x,y) -> Chart.Scatter3D(x,y,Array.init x.Length (fun _ -> Random.rndgen.NextFloat() / 3.),StyleParam.Mode.Markers)

let mvnChart = 
    [surface;samples]
    |> Chart.combine
    |> Chart.withTitle "Bivariate normal distribution with sampling"

Students t distribution

let studentTParams = [(0.,1.,1.);(0.,1.,2.);(0.,1.,5.);]
let xStudentT = [-10. ..0.1.. 10.]

let pdfStudentT mu tau dof = 
    xStudentT 
    |> List.map (Continuous.StudentT.PDF mu tau dof)
    |> List.zip xStudentT


let cdfStudentT mu tau dof = 
    xStudentT 
    |> List.map (Continuous.StudentT.CDF  mu tau dof)
    |> List.zip xStudentT


let v =
    studentTParams
    |> List.map (fun (mu,tau,dof) -> Chart.Spline(pdfStudentT mu tau dof,Name=sprintf "mu=%.1f tau=%.1f dof=%.1f" mu tau dof,ShowMarkers=false))
    |> Chart.combine
    |> Chart.withTemplate ChartTemplates.lightMirrored

F distribution

The F distribution or Fisher distribution, also known as Fisher-Snedecor distribution, is a continuous probability distribution. An F-distributed random variable results from the quotient of two Chi-square-distributed random variables each divided by the associated number of degrees of freedom. The F-distribution has two independent degrees of freedom(dof) as parameters, and thus forms a two-parameter distribution family.

Generally speaking, the F-tests and the resulting F-Distribution is utilized for comparing multiple levels of independent variables with multiple groups. In practice, it is most commonly used to compare the variances within a group to the variance between different groups, as seen in the Analysis of varaince.

let fParams = [(2.,1.);(5.,2.);(10.,1.);(100.,100.)]
let xF = [0. .. 1. .. 5.]

let pdfF a b = 
    xF 
    |> List.map (Continuous.F.PDF a b)
    |> List.zip xF

let fPDFs =
    fParams
    |> List.map (fun (a,b) -> Chart.Line(pdfF a b,Name=sprintf "dof1=%.1f dof2=%.1f" a b,LineWidth=3.) )
    |> Chart.combine
    |> Chart.withTemplate ChartTemplates.lightMirrored
    |> Chart.withTitle "Different F-Distributions PDFs, x=[0,5]"

let cdfF a b = 
    xF 
    |> List.map (Continuous.F.CDF a b)
    |> List.zip xF

let fCDFs =
    fParams
    |> List.map (fun (a,b) -> Chart.Line(cdfF a b,Name=sprintf "dof1=%.1f dof2=%.1f" a b,LineWidth=3.) )
    |> Chart.combine
    |> Chart.withTemplate ChartTemplates.lightMirrored
    |> Chart.withTitle "Different F-Distributions CDFs, x=[0,5]"

Discrete

Bernoulli distribution

A Bernoulli distribution is a (..) random experiment that has only two outcomes (usually called a "Success" or a "Failure"). For example, the probability of getting a heads (a "success") while flipping a coin is 0.5. The probability of "failure" is 1 – P (1 minus the probability of success, which also equals 0.5 for a coin toss). It is a special case of the binomial distribution for n = 1. In other words, it is a binomial distribution with a single trial (e.g. a single coin toss).
~ by statisticshowto

Mathematically, "success" and "failure" are represented as 1.0 and 0.0, respectively.

It is defined by one parameter B(p):

• p = probability of success

Example: A weighted coin with a probability of 0.6 to land on tails. Most bernoulli distribution calculations are rather intuitive:

BernA: What is the mean of a bernoulli distribution with the weighted coin?

BernB: What is the probability to land on heads?

open FSharp.Stats
open FSharp.Stats.Distributions

// Assumes "tails" to be success
let bernoulli = Discrete.Bernoulli.Init 0.6

// BernA: What is the mean of a bernoulli distribution with the weighted coin?
let bernA = bernoulli.Mean
// Output: 0.6
// Altough the bernoulli distribution can never return 0.6 (only 0.0 or 1.0) on average it will return heads at the same probability it has to land on heads.

// BernB: What is the probability to land on heads?
let bernB = bernoulli.PMF 0
// Output: 0.4
// Again: Heads = 0.0 = failure and tails = 1.0 = success. 

let plotBernoulli =
    [0; 1]
    |> List.map (fun x -> x, bernoulli.PMF x)
    |> Chart.Column
    |> Chart.withTemplate ChartTemplates.lightMirrored
    |> Chart.withTitle "B(0.6)"

Binomial distribution

The binomial distribution describes the probability, that under a given success probability and a given number of draws an event occurs exactly k times (with replacement).

It is defined by two parameters B(n,p):

• n = number of draws
• p = probability of success

Example: The school bus is late with a probability of 0.10.

BinoA: What is the probability of running late exactly 5 times during a 30 day month?

BinoB: What is the probability of running late for a maximum of 5 times?

BinoC: What is the probability of running late for at least 5 times?

open FSharp.Stats
open FSharp.Stats.Distributions

// Creates a binomial distribution with n=30 and p=0.90 
let binomial = Discrete.Binomial.Init 0.1 30

// BinoA: What is the probability of running late exactly 5 times during a 30 day month?
let binoA = binomial.PMF 5
// Output: 0.1023 = 10.23 %

// BinoB: What is the probability of running late for a maximum of 5 times?
let binoB = binomial.CDF 5.
// Output: 0.9268 = 92.68 %

// BinoC: What is the probability of running late for at least 5 times?
let binoC = 1. - (binomial.CDF 4.)
// Output: 0.1755 = 17.55 %

let plotBinomial =
    [0..30]
    |> List.map (fun x -> x,binomial.PMF x)
    |> Chart.Column
    |> Chart.withTemplate ChartTemplates.lightMirrored
    |> Chart.withTitle "B(30,0.1)"

Hypergerometric distribution

The hypergeometric distribution describes the probability, that under a given number of success and failure events and a given number of draws an event occurs exactly k times (without replacement).

It is defined by three parameters:

• N = finite population representing the total number of events.
• K = number of success events in this population.
• n = number of draws from the population.

Example: You participate in a lottery, where you have to choose 6 numbers out of 49. The lottery queen draws 6 numbers randomly, where the order does not matter.

HypA: What is the probability that your 6 numbers are right?

HypB: How to simulate artificial draws from the distribution?

HypC: What is the probability that you have at least 3 right ones?

HypD: What is the probability that you have a maximum of 3 right ones?

// Creates a hypergeometric distribution with N=49, K=6, n=6.
let hyper = Discrete.Hypergeometric.Init 49 6 6

// HypA: What is the probability that your 6 numbers are right?
let hypA = hyper.PMF 6
// Output: 7.15-08

// HypB: How to simulate artificial draws from the distribution?
let hypB = hyper.Sample()
// Output: Number of success events randomly sampled from the distribution.

// HypC: What is the probability that you have at least 3 right ones?
// CDF is implemented to calculate P(X <= k)
let hypC = 1. - hyper.CDF 2
// Output: 0.01864 = 1.86 %

// HypD: What is the probability that you have a maximum of 3 right ones? 
let hypD = hyper.CDF 3
// Output: 0.99901 = 99.90 %

Poisson distribution

The poisson distribution describes what the probability for a number of events is, occurring in a certain time, area, or volume.

It is defined by just one parameters Po(lambda) where lambda is the average rate.

Example: During a year, a forest is struck by a lightning 5.5 times.

PoA: What is the probability that the lightning strikes exactly 3 times?

PoB: What is the probability that the lightning strikes less than 2 times?

PoC: What is the probability that the lightning strikes more than 7 times?

// Creates a poisson distribution with lambda=  .
let poisson = Discrete.Poisson.Init 5.5

// PoA: What is the probability that the lightning strikes exactly 3 times?
let poA = poisson.PMF 3
// Output: 0.11332 = 11.33 %

// PoB: What is the probability that the lightning strikes less or equal to 2 times?
let poB = 
    // CDF not implemented yet
    //poisson.CDF 2.
    [0 .. 2] |> List.sumBy poisson.PMF
    // Output: 0.088376 = 8.84 %
    
// PoC: What is the probability that the lightning strikes more than 7 times?
let poC = 1. -  poisson.CDF 6.
// Output: Not implemented yet. Manual addition necessary

let plotPo =
    [0..20]
    |> List.map (fun x -> x,poisson.PMF x)
    |> Chart.Column
    |> Chart.withTemplate ChartTemplates.lightMirrored
    //|> Chart.withSize(600.,400.)
    |> Chart.withTitle "Po(5.5)"

Gamma distribution

let gammaParams = [(1.,2.);(2.,2.);(3.,2.);(5.,1.);(9.0,0.5);(7.5,1.);(0.5,1.);]
let xgamma = [0. ..0.1.. 20.]

let pdfGamma a b = 
    xgamma 
    |> List.map (Continuous.Gamma.PDF a b)
    |> List.zip xgamma

let gammaPlot =
    gammaParams
    |> List.map (fun (a,b) -> Chart.Point(pdfGamma a b,Name=sprintf "a=%.1f b=%.1f" a b) )//,ShowMarkers=false))
    |> Chart.combine
    |> Chart.withTemplate ChartTemplates.lightMirrored

let cdfGamma a b = 
    xgamma 
    |> List.map (Continuous.Gamma.CDF a b)
    |> List.zip xgamma

let gammaCDFPlot=
    gammaParams
    |> List.map (fun (a,b) -> Chart.Spline(cdfGamma a b,Name=sprintf "a=%.1f b=%.1f" a b) )//,ShowMarkers=false))
    |> Chart.combine
    |> Chart.withTemplate ChartTemplates.lightMirrored

Negative binomial distribution

Note: There is no unique definition of the negative binomial distribution. In definition A the random variable x is the number of trials, while in definition B the random variable x is the number of failures. The definition of r (number of successes) and p (probability of each bernoulli trial) is the same for both definitions.

In FSharp.Stats both definitions are implemented. The number of trials (A) can be modelled with negativeBinomial_trials, while the number of failures can be modelled with negativeBinomial_failures. For further discussion about this issue check out #256.

A - Negative binomial distribution (trials)

The distribution models the number of trials needed (x) to get the r^th success in repeated independent Bernoulli trials. Until the (x-1)^th trial, (r-1) successes must be accomplished, which matches the standard binomial distribution. Therefore, to get the r^th success in the x^th trial, you have to multiply Binom(p,x-1,r-1) by p.

B - Negative binomial distribution (faiures)

The distribution models the number of failures (k) prior to the r^th success in repeated independent Bernoulli trials. Until the (k+r-1)^th trial, (r-1) successes must be accomplished, which matches the standard binomial distribution. Therefore, to get the r^th success after k^th failures, you have to multiply Binom(p,k+r-1,r-1) by p.

The PDF and CDF of both distributions can be converted into each other by:

(negativeBinomial_trials   r p).PMF x = (negativeBinomial_failures r p).PMF (x-r)
(negativeBinomial_failures r p).PMF k = (negativeBinomial_trials   r p).PMF (k+r)

(negativeBinomial_trials   r p).CDF x = (negativeBinomial_failures r p).CDF (x-r)
(negativeBinomial_failures r p).CDF k = (negativeBinomial_trials   r p).CDF (k+r)

Task1: What is the probability of obtaining the third success on the 10 trial. The individual success probability is 0.09.

//number of trials
let x = 10

//probability of the independent bernoulli trials
let p = 0.09

//number of successes
let r = 3

// number of failures 
let k = x - 3

//Solution A:
let negB_A = Discrete.NegativeBinomial_trials.PMF r p x
//returns 0.01356187619

//Solution B:
let negB_B = Discrete.NegativeBinomial_failures.PMF r p k
//returns 0.01356187619

Comparison of PMF

let negBinom_trials =
    Distributions.Discrete.NegativeBinomial_trials.Init 3 0.3

let negBinom_failures = 
    Distributions.Discrete.NegativeBinomial_failures.Init 3 0.3

negBinom_trials.CDF 1


let pmfComparison = 
    let pmfs_trials   = [0..30] |> List.map (fun x -> x,negBinom_trials.PMF x)
    let pmfs_failures = [0..30] |> List.map (fun x -> x,negBinom_failures.PMF x)
    [
         Chart.Column (pmfs_trials, Opacity=0.5, Name= "trials (x)")
         Chart.Column (pmfs_failures, Opacity=0.5, Name= "failures (k)")
    ]
    |> Chart.combine
    |> Chart.withLayoutStyle(BarMode=StyleParam.BarMode.Overlay)
    |> Chart.withTemplate ChartTemplates.lightMirrored
    |> Chart.withTitle "PMF"
    |> Chart.withXAxisStyle "x or k"
    |> Chart.withYAxisStyle "P(X=x) or P(X=k)"

You can clearly see, that the PMF distributions are shifted according to the number of successes because: \(failures (k) = trials (x) - successes (r)\).

Comparison of CDF

let cdfComparison = 
    let cdfs_trials   = [0..30] |> List.map (fun x -> x,negBinom_trials.CDF x)
    let cdfs_failures = [0..30] |> List.map (fun x -> x,negBinom_failures.CDF x)
    [
         Chart.Column (cdfs_trials, Opacity=0.5, Name= "trials (x)")
         Chart.Column (cdfs_failures, Opacity=0.5, Name= "failures (k)")
    ]
    |> Chart.combine
    |> Chart.withLayoutStyle(BarMode=StyleParam.BarMode.Overlay)
    |> Chart.withTemplate ChartTemplates.lightMirrored
    |> Chart.withTitle "CDF"
    |> Chart.withXAxisStyle "x or k"
    |> Chart.withYAxisStyle "P(X=x) or P(X=k)"

You can clearly see, that the CDF distributions are shifted according to the number of successes because: \(failures (k) = trials (x) - successes (r)\).

Multinomial distribution

The multinomial distribution is a generic version of the binomial distribution. While for binomial, the probabilities for a single success state is of interest, the multinomial distribution deals with multiple exact success events.

Example: There are people from 3 different towns: 3 from town A, 7 from town B, and 20 from town C. When 5 people are randomly chosen, what is the probability, to obtain exactly 1 from town A, 1 from town B, and 3 from town C? The individual success probabilities can be easily accessed p(A)=5/30, p(B)=7/30, and p(C)=20/30.

// probabilities for all individual success states 
let multiNomProb = vector [(3./30.); (7./30.); (20./30.)]

// the success combination that is of interest
let multiNomKs   = Vector.Generic.ofList [1; 1; 3]

// gives the probability of obtaining exactly the pattern 1,1,3
let mNom = Discrete.Multinomial.PMF multiNomProb multiNomKs

Relation to binomial distribution

Input vectors of length 2 correspond to the binomial distribution. The following examples are identical. While for binomial it is required to input the total number (n), for the multinomial distribution you have to give the corresponding anto-probability:

let mNom_bin_A = (Discrete.Binomial.PMF 0.123 200 20)
let mNom_bin_B = Discrete.Multinomial.PMF (vector [|0.123; 0.877|]) (Vector.Generic.ofArray [|20; 180|])

mNom_bin_A //0.0556956956889893
mNom_bin_B //0.0556956956889898

A cumulative density function (CDF) is not defined properly, as you do not know which success statement is of interest. If you want to investigate how probable it is to see at least 3 people of town C you would have to calculate the sum all possible combinations that result in this constellation:

Discrete.Multinomial.PMF multiNomProb [1;1;3]
Discrete.Multinomial.PMF multiNomProb [2;0;3]
Discrete.Multinomial.PMF multiNomProb [0;2;3]
Discrete.Multinomial.PMF multiNomProb [1;0;4]
Discrete.Multinomial.PMF multiNomProb [0;1;4]
Discrete.Multinomial.PMF multiNomProb [0;0;5]

The cumulative probability is 0.7901234.

It may be that in future, a dedicated CDF functionality is added that requests the index of the success state of interest and sums up all possible combination probabilities.

Empirical

You can create empirically derived distributions and sample randomly from these. In this example 100,000 random samples from a student t distribution are drawn and visualized

// Randomly taken samples; in this case from a gaussian normal distribution.
let sampleDistribution = 
    let template = Continuous.Normal.Init 5. 2.
    Seq.init 100000 (fun _ -> template.Sample())

//creates an empirical distribution with bandwidth 0.1
let empiricalDistribution = 
    EmpiricalDistribution.create 0.1 sampleDistribution

Categorical data

You also can create probability mass functions (PMF) from categorical (nominal) data. You can check for just the elements present in your input sequence, or include elements of a predefined set.

let myText = 
    "Hello World, I am a test Text with all kind of Characters!112"

// Define your set of characters that should be checked for
// Any character that is not present in these sets is neglected
let mySmallAlphabet = "abcdefghijklmnopqrstuvwxyz" |> Set.ofSeq
let myAlphabet      = "abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ" |> Set.ofSeq
let myAlphabetNum   = "abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ0123456789" |> Set.ofSeq

These alphabets can be used to create the probability maps.

// takes the characters and determines their probabilities without considering non-existing characters
let myFrequencies0 = EmpiricalDistribution.createNominal() myText

// takes upper and lower case characters and determines their probability
let myFrequencies1 = EmpiricalDistribution.createNominal(Template=myAlphabet) myText

// takes upper and lower case characters and numbers and determines their probability
let myFrequencies2 = EmpiricalDistribution.createNominal(Template=myAlphabetNum) myText

// takes only lower case characters and determines their probability
// The big M is neglected because it is not present in the template alphabet.
let myFrequencies3 = EmpiricalDistribution.createNominal(Template=mySmallAlphabet) myText

An additional field for transforming the input sequence may be beneficial if it does not matter if an character is lower case or upper case:

// converts all characters to lower case characters and determines their probability
let myFrequencies4 = EmpiricalDistribution.createNominal(Template=mySmallAlphabet,Transform=System.Char.ToLower) myText

With a template set, you can access the probability of 'z' even if it is not in your original input sequence.

myFrequencies3.['z'] //results in 0.0

Visualization of the PMF

let categoricalDistribution = 
    [
    Chart.Column(myFrequencies0 |> Map.toArray,"0_noTemplate")    |> Chart.withYAxisStyle "probability"
    Chart.Column(myFrequencies1 |> Map.toArray,"1_bigAlphabet")   |> Chart.withYAxisStyle "probability"
    Chart.Column(myFrequencies2 |> Map.toArray,"2_numAlphabet")   |> Chart.withYAxisStyle "probability"
    Chart.Column(myFrequencies3 |> Map.toArray,"3_smallAlphabet") |> Chart.withYAxisStyle "probability"
    Chart.Column(myFrequencies4 |> Map.toArray,"4_toLower + smallAlphabet") |> Chart.withYAxisStyle "probability"
    ]
    |> Chart.Grid(4,1)
    |> Chart.withTemplate ChartTemplates.lightMirrored

Distribution merging

You can merge two distributions by using Empirical.merge, subroutines like Empirical.add, or the generic function Empirical.mergeBy.

Merging two distributions leads to a combined distribution. If keys are present in both distributions the value at distA is superseded with the value at distB.

Please note, that when handling continuous data, the binning of both input distributions must be identical! When using categorical data, the binning does not matter and the parameter can be set to true.

let a =
    [("k1",1);("k2",3)]
    |> Map.ofList

let b =
    [("k2",2);("k3",4)]
    |> Map.ofList

let mergedDist = Empirical.merge true a b

Adding two distributions leads to a combined distribution. If keys are present in both distributions the values at distA and distB are added.

let addedDist = Empirical.add true a b

A custom merging function can be defined:

let customDist = Empirical.mergeBy true (fun valueA valueB -> valueA * valueB) a b

Density estimation

let nv = Array.init 1000 (fun _ -> Distributions.Continuous.Normal.Sample 5. 2.)

let xy = KernelDensity.estimate KernelDensity.Kernel.gaussian 1.0 nv

let kernelChart = 
    Chart.SplineArea xy
    |> Chart.withTemplate ChartTemplates.lightMirrored
 

Distance

In this example we will calculate the Wasserstein Metric between 3 distributions. One can imagine this metric as the amount of work needed to move the area (pile of dirt) of one distribution to the area of the other. That's why it's also called Earth Movers Distance.

Be aware that this distance measure is dependent on the actual absolute values of the distributions.

let distribution1 = 
    let normal = Continuous.Normal.Init 300. 15.
    Array.init 1000 (fun _ -> normal.Sample())
let distribution2 =
    let normal = Continuous.Normal.Init 350. 20.
    Array.init 500 (fun _ -> normal.Sample())
let distribution3 =
    let normal = Continuous.Normal.Init 500. 20.
    Array.init 1000 (fun _ -> normal.Sample())

let pilesOfDirt =
    [
    Chart.Histogram(distribution1,Name = "Distribution1")
    Chart.Histogram(distribution2,Name = "Distribution2")
    Chart.Histogram(distribution3,Name = "Distribution3")
    ]
    |> Chart.combine

let distance1and2 = Distance.OneDimensional.wassersteinDistance distribution1 distribution2

49.56468837

let distance1and3 = Distance.OneDimensional.wassersteinDistance distribution1 distribution3

199.0156561

As expected the distance between Distribution 1 and 2 is the lowest:

let distributions = 
    [|distribution1;distribution2;distribution3|]

let mapColor min max value = 
    let proportionR = 
        ((255. - 200.) * (value - min) / (max - min))
        |> int
        |> fun x -> 255 - x
    let proportionG = 
        ((255. - 200.) * (value - min) / (max - min))
        |> int
        |> fun x -> 255 - x
    let proportionB = 
        ((255. - 200.) * (value - min) / (max - min))
        |> int
        |> fun x -> 255 - x
    Color.fromARGB 1 proportionR proportionG proportionB

let distancesTable =
    let headerColors = ["white";"#1f77b4";"#ff7f0e";"#2ca02c"] |> List.map Color.fromString |> Color.fromColors
    let distances = 
        distributions
        |> Array.map (fun x ->
            distributions
            |> Array.map (fun y ->
                Distance.OneDimensional.wassersteinDistance x y
                |> (sprintf "%.2f")
            )
        )
        |> Array.transpose 
        |> Array.append [|[|"Distribution1";"Distribution2";"Distribution3"|]|]
        |> Array.transpose 
    let cellColors = 
        distances
        |> Array.mapi (fun i a ->
            a
            |> Array.mapi (fun j v -> 
                if j = 0 then 
                    if i = 0 then Color.fromHex "#1f77b4"
                    elif i = 1 then Color.fromHex "#ff7f0e"
                    else Color.fromHex "#2ca02c"
                else 
                    mapColor 0. 200. (float v)
            )
        )
        |> Array.transpose
        |> Seq.map Color.fromColors
        |> Color.fromColors
 
    Chart.Table(
        ["";"Distribution1";"Distribution2";"Distribution3"],         
        distances,
        HeaderFillColor = headerColors,
        CellsFillColor = cellColors)