Testing.fs

module TestingTests
open Expecto
open System

open FsMath
open FSharp.Stats
open FSharp.Stats.Testing

open TestExtensions
open FSharp.Stats.Testing.SAM
open Deedle

[<Tests>]
let testPostHocTests =
    //Tests taken from:
    //https://www.icalcu.com/stat/anova-tukey-hsd-calculator.html
    testList "Testing.PostHoc" [
        (*
        // Test ommitted due to extremely long runtime of CodeCov.
        testCase "tukeyHSD" <| fun () ->
            let dataA = [|3.;3.;4.;5.;2.;5.;5.;4.;4.;2.;2.;2.;4.;3.;5.;3.;4.;5.;3.;5.;                   |]
            let dataB = [|10.;7.;9.;6.;7.;7.;6.;7.;10.;7.;8.;8.;8.;6.;10.;9.;9.;6.;9.;8.;                |]
            let dataC = [|6.;5.;6.;4.;4.;6.;1.;4.;6.;5.;4.;7.;4.;2.;1.;1.;3.;4.;5.;3.;                   |]
            let dataD = [|10.;5.;6.;5.;8.;5.;6.;9.;3.;10.;5.;9.;5.;5.;6.;10.;9.;6.;9.;10.;               |]
            let dataE = [|14.;17.;14.;13.;18.;12.;17.;11.;12.;11.;12.;10.;17.;19.;18.;18.;15.;14.;18.;16.|]

            let data = [|dataA;dataB;dataC;dataD;dataE|]
                
            let contrastMatrix = 
                [|                
                    //[|-1.;1.;0.;0.;0.;|] pvalue = zero
                    [|-1.;0.;1.;0.;0.;|]
                    [|-1.;0.;0.;1.;0.;|]
                    //[|-1.;0.;0.;0.;1.;|] pvalue = zero
                    [|0.;-1.;1.;0.;0.;|]
                    [|0.;-1.;0.;1.;0.;|]
                    //[|0.;-1.;0.;0.;1.;|] pvalue = zero
                    [|0.;0.;-1.;1.;0.;|]
                    //[|0.;0.;-1.;0.;1.;|] pvalue = zero
                    //[|0.;0.;0.;-1.;1.;|] pvalue = zero
                |]

            let pValues = 
                PostHoc.tukeyHSD contrastMatrix data 
                |> Array.map (fun x -> x.Significance)

            //pvalues from R: TUKEY <- TukeyHSD(x=ANOVA, 'data$treatment', conf.level=0.95)
            let rpval = [0.9685630;0.0000045;0.0000003;0.7072882;0.0000618]
                
            Expect.floatClose Accuracy.low rpval.[0] pValues.[0] "p values should be equal."
            Expect.floatClose Accuracy.low rpval.[1] pValues.[1] "p values should be equal."
            Expect.floatClose Accuracy.low rpval.[2] pValues.[2] "p values should be equal."
            Expect.floatClose Accuracy.low rpval.[3] pValues.[3] "p values should be equal."
            Expect.floatClose Accuracy.low rpval.[4] pValues.[4] "p values should be equal."
        *)
        testCase "dunnett" <| fun () ->
            let data = 
                [|
                    [|1.84;2.49;1.50;2.42;|]
                    [|2.43;1.85;2.42;2.73;|]
                    [|3.95;3.67;3.23;2.31;|]
                    [|3.21;3.20;2.32;3.30;|]
                    [|3.21;3.13;2.32;3.30;3.20;2.42;|]
                |]

            //first sample is control
            let contrastMatrix = 
                [|                
                    [|-1.;1.;0.;0.;0.|]
                    [|-1.;0.;1.;0.;0.|]
                    [|-1.;0.;0.;1.;0.|]
                    [|-1.;0.;0.;0.;1.|]
                |]

            let dunnettResult = 
                PostHoc.dunnetts contrastMatrix data Tables.dunnettsTwoSided095

            //result from: SPSS Dunnett's test version 27
            let pval = [0.811;0.010;0.050;0.049]
            let dmean = [0.295;1.2275;0.945;0.8675]
                
            Expect.equal dunnettResult.[0].Significance (pval.[0]<0.05) "Significance should be equal."
            Expect.equal dunnettResult.[1].Significance (pval.[1]<0.05) "Significance should be equal."
            Expect.equal dunnettResult.[2].Significance (pval.[2]<0.05) "Significance should be equal."
            Expect.equal dunnettResult.[3].Significance (pval.[3]<0.05) "Significance should be equal."
            Expect.floatClose Accuracy.high dunnettResult.[0].L dmean.[0] "Mean differences should be equal."
            Expect.floatClose Accuracy.high dunnettResult.[1].L dmean.[1] "Mean differences should be equal."
            Expect.floatClose Accuracy.high dunnettResult.[2].L dmean.[2] "Mean differences should be equal."
            Expect.floatClose Accuracy.high dunnettResult.[3].L dmean.[3] "Mean differences should be equal."
    ]

[<Tests>]
let hTestTests = 
    // H-Test with ties tested against r implementation kruskal.test(weight ~ group, data = my_data)
    let groupA = [4.17; 5.18;  5.18;  6.11;  4.50;  4.61;  5.17;  4.53;  5.33;  5.18;] 
    let groupB = [4.81; 4.17;  4.41;  3.59;  5.87;  3.83;  6.03;  4.89;  4.32;  4.69;] 
    let groupC = [6.31; 5.12;  5.00;  5.00;  5.00;  5.29;  5.00;  6.15;  5.80;  5.26;]   
    let samples = [groupA;groupB;groupC]
        
    // calculation of the H test 
    let hResult = 
        HTest.createHTest samples 
        
    testList "Testing.HTest" [
        testCase "createHTest" <| fun () -> 
            Expect.isTrue (0.03781 = Math.Round(hResult.PValueRight,5)) "pValue should be equal."
            Expect.isTrue (6.5502  = Math.Round(hResult.Statistic,4)) "statistic should be equal."
                
    ]

[<Tests>]
let friedmanTestTests = 
    // Friedman-Test testes against dataset from https://www.methodenberatung.uzh.ch/de/datenanalyse_spss/unterschiede/zentral/friedman.html#3.2._Ergebnisse_des_Friedman-Tests and p-values obtained from distcalc and https://www.socscistatistics.com/pvalues/chidistribution.aspx 
    let A = [|275.;273.;288.;273.;244.|]
    let B = [|292.;283.;284.;285.;329.|]
    let C = [|281.;274.;298.;270.;252.|]
    let D = [|284.;275.;271.;272.;258.|]
    let E = [|285.;294.;307.;278.;275.|]
    let F = [|283.;279.;301.;276.;279.|]
    let G = [|290.;265.;298.;291.;295.|]
    let H = [|294.;277.;295.;290.;271.|]
    let I = [|300.;304.;293.;279.;271.|]
    let J = [|284.;297.;352.;292.;284.|]
    let samples = seq{A;B;C;D;E;F;G;H;I;J}

    // modified dataset from UZH for 3x equal ranks 
    let A2 = [|275.;273.;288.;273.;273.|]
    let B2 = [|292.;283.;284.;285.;329.|]
    let C2 = [|281.;274.;298.;270.;252.|]
    let D2 = [|284.;275.;271.;272.;258.|]
    let E2 = [|285.;294.;307.;278.;275.|]
    let F2 = [|283.;279.;301.;276.;279.|]
    let G2 = [|290.;265.;298.;291.;295.|]
    let H2 = [|294.;277.;295.;290.;271.|]
    let I2 = [|300.;304.;293.;279.;271.|]
    let J2 = [|284.;297.;284.;292.;284.|]
    let samples2 = seq{A2;B2;C2;D2;E2;F2;G2;H2;I2;J2}


    //calculation of friedman test
    let friedmanResult1 = 
        FriedmanTest.createFriedmanTest samples 
        
    let friedmanResult2 = 
        FriedmanTest.createFriedmanTest samples2 

    testList "Testing.FriedmanTest" [
        testCase "createFriedmanTest2equal" <| fun () -> 
            Expect.floatClose Accuracy.low friedmanResult1.Statistic 13.259 "statistics should be equal."
            Expect.floatClose Accuracy.low friedmanResult1.PValueRight 0.010077 "pValue should be equal."
        testCase "createFriedmanTest3equal" <| fun () -> 
            Expect.floatClose Accuracy.low friedmanResult2.Statistic 9.738 "statistics should be equal."
            Expect.floatClose Accuracy.low friedmanResult2.PValueRight 0.04508 "pValue should be equal."
        ]

[<Tests>]
let wilcoxonTestTests = 
    // tested against SciPy Version 1.7.1
    let before = seq{78.;24.;64.;45.;64.;52.;30.;50.;64.;50.;78.;22.;84.;40.;90.;72.}
    let after = seq{78.;24.;62.;48.;68.;56.;25.;44.;56.;40.;68.;36.;68.;20.;58.;32.}
    let differences = seq{0.;0.;2.;-3.;-4.;-4.;5.;6.;8.;10.;10.;-14.;16.;20.;32.;40.}
    // with continuity correction:
    let wilcoxon1 = WilcoxonTest.createWilcoxonTest before after true 
    let wilcoxon2 = WilcoxonTest.createWilcoxonTest before after false
    let wilcoxon3 = WilcoxonTest.createWilcoxonTestFromDifferences differences true 
    let wilcoxon4 = WilcoxonTest.createWilcoxonTestFromDifferences differences false
    // ties of size 3 (exposes the rank-multiplication bug): differences [1,1,1,2,-3]
    // expected values computed with SciPy: wilcoxon([1,1,1,2,-3], mode='approx', correction=False/True)
    let tieDiffs3 = seq{1.;1.;1.;2.;-3.}
    let wilcoxon5 = WilcoxonTest.createWilcoxonTestFromDifferences tieDiffs3 false
    let wilcoxon6 = WilcoxonTest.createWilcoxonTestFromDifferences tieDiffs3 true

    testList "Testing.WilcoxonTest" [
        testCase "wilcoxonWithCorrection" <| fun () -> 
            Expect.floatClose Accuracy.low wilcoxon1.PValueTwoTailed 0.0382 "pValue should be equal."
        testCase "wilcoxonWithoutCorrection" <| fun () -> 
            Expect.floatClose Accuracy.low wilcoxon2.PValueTwoTailed 0.03537 "pValue should be equal."
        testCase "wilcoxonDifferencesWithCorrection" <| fun () -> 
            Expect.floatClose Accuracy.low wilcoxon3.PValueTwoTailed 0.0382 "pValue should be equal."
        testCase "wilcoxonDifferencesWithoutCorrection" <| fun () -> 
            Expect.floatClose Accuracy.low wilcoxon4.PValueTwoTailed 0.03537 "pValue should be equal."
        testCase "wilcoxonOneSidedWithCorrection" <| fun () -> 
            Expect.floatClose Accuracy.low wilcoxon1.PValueLeft 0.019102 "pValue should be equal"
        testCase "wilcoxonOneSidedWithoutCorrection" <| fun () -> 
            Expect.floatClose Accuracy.low wilcoxon2.PValueRight 0.9823 "pValue should be equal"
        testCase "wilcoxonTieSize3WithoutCorrection" <| fun () ->
            // SciPy: wilcoxon([1,1,1,2,-3], mode='approx', correction=False) => pvalue ≈ 0.492207
            Expect.floatClose Accuracy.low wilcoxon5.PValueTwoTailed 0.4922 "pValue should match SciPy with ties of size 3"
        testCase "wilcoxonTieSize3WithCorrection" <| fun () ->
            // SciPy: wilcoxon([1,1,1,2,-3], mode='approx', correction=True) => pvalue ≈ 0.582702
            Expect.floatClose Accuracy.low wilcoxon6.PValueTwoTailed 0.5827 "pValue should match SciPy with ties of size 3 and correction"
            
        ]


[<Tests>]
let tTestTests = 
    // tested in SPSS version 27
    let groupA = [|-5.;-3.;-3.;-4.;-5.;|] 
    let groupB = [|-2.;-4.;-4.;-6.;-6.;-6.;-5.;|] 
    let groupC = [|-3.;-7.;-8.;-4.;-2.; 1.;-1.;|]   
    let groupD = [|1.;-1.;0.;2.;2.;|]   
        
    let meanA = Seq.mean groupA
    let meanB = Seq.mean groupB
    let varA = Seq.var groupA
    let varB = Seq.var groupB
    let nA = float (Seq.length groupA)
    let nB = float (Seq.length groupB)

    // calculation of the H test 
    let tTest1 = TTest.twoSample true groupA groupB
    let tTest2 = TTest.twoSampleFromMeanAndVar true (meanA,varA,nA) (meanB,varB,nB) 
    let tTest3 = TTest.twoSample false groupA groupB
    let tTest4 = TTest.oneSample groupD 0.5

    testList "Testing.TTest" [
        testCase "twoSample" <| fun () -> 
            Expect.floatClose Accuracy.low tTest1.PValue 0.377 "pValue should be equal."
            Expect.floatClose Accuracy.low tTest1.DegreesOfFreedom 10. "df should be equal."
            Expect.floatClose Accuracy.low tTest1.Statistic 0.924 "t statistic should be equal."
            Expect.floatClose Accuracy.low tTest3.PValue 0.345 "pValue should be equal."
            Expect.floatClose Accuracy.low tTest3.DegreesOfFreedom 9.990 "df should be equal."
        
        testCase "twoSampleFromMeanAndVar" <| fun () -> 
            Expect.equal tTest1 tTest2 "results should be equal."
            
           // tested with R function (t.test(batcha, batchb, var.equal=TRUE))
            let sample1 = [1.;2.;3;]
            let sample2 = [1.;3.;5.;7;]
            let mean = Seq.mean sample1
            let mean2 = Seq.mean sample2
            let var1 = Seq.var sample1
            let var2 = Seq.var sample2
            let ttestTwoSample = Testing.TTest.twoSampleFromMeanAndVar true (mean,var1,3) (mean2,var2,4) 
            let expectedPval = 0.26716219523142071
            let expectedStatistic = -1.24837556786471859
            Expect.floatClose Accuracy.high ttestTwoSample.PValue expectedPval "pValue should be equal."
            Expect.floatClose Accuracy.high ttestTwoSample.Statistic expectedStatistic "t statistic should be equal."

            let sample3 = [-1.;3.;-5.;7;]
            let mean3 = Seq.mean sample2
            let var3 = Seq.var sample2
            let ttestTwoSample2 = Testing.TTest.twoSampleFromMeanAndVar true (mean,var1,3) (mean3,var3,4) 
            let expectedPval2 = 0.75954440793496059
            let expectedStatistic2 = -0.323310403056781825
            Expect.floatClose Accuracy.high ttestTwoSample2.PValue expectedPval "pValue should be equal."
            Expect.floatClose Accuracy.high ttestTwoSample2.Statistic expectedStatistic "t statistic should be equal."

            let ttestTwoSample3 = Testing.TTest.twoSampleFromMeanAndVar true (nan,var2,3) (mean2,var2,4) 
            Expect.isTrue (nan.Equals(ttestTwoSample3.PValue)) "pValue should be nan."
            Expect.isTrue (nan.Equals(ttestTwoSample3.Statistic)) "t statistic should be nan."
        
        
        testCase "oneSample" <| fun () -> 
            Expect.floatClose Accuracy.low tTest4.PValue 0.634 "pValue should be equal."
            Expect.equal tTest4.DegreesOfFreedom 4. "df should be equal."
            Expect.floatClose Accuracy.low tTest4.Statistic 0.514 "t statistic should be equal."
        
        //tested with R function (t.test(c(-1,-2,-3), mu = -3, alternative = "two.sided"))
        testCase "oneSampleFromMeanandStDev" <| fun () -> 
            let sample = [1.;2.;3;]
            let mean = Seq.mean sample
            let stdev = Seq.stDev sample
            let ttest = Testing.TTest.oneSampleFromMeanAndStDev(mean,stdev,3) -3.
            let expectedPval = 0.013072457560346513
            let expectedStatistic = 8.6602540378443873
            Expect.floatClose Accuracy.high ttest.PValue expectedPval "pValue should be equal."
            Expect.floatClose Accuracy.high ttest.Statistic expectedStatistic "t statistic should be equal."

            let sample = [-1.;-2.;-3;]
            let mean3 = Seq.mean sample
            let stdev1 = Seq.stDev sample
            let ttest2 = Testing.TTest.oneSampleFromMeanAndStDev(mean3,stdev1,3) 0.
            let expectedPval1 = 0.074179900227448525
            let expectedStatistic2 = -3.46410161513775483
            Expect.floatClose Accuracy.high ttest2.PValue expectedPval1 "pValue should be equal."
            Expect.floatClose Accuracy.high ttest2.Statistic expectedStatistic2 "t statistic should be equal."

            let mean2 = nan
            let ttest3 = Testing.TTest.oneSampleFromMeanAndStDev(mean2,stdev,3) 0.
            Expect.isTrue (nan.Equals(ttest3.PValue)) "pValue should be nan."
            Expect.isTrue (nan.Equals(ttest3.Statistic)) "t statistic should be nan."

        testCase "twoSamplePaired" <| fun () -> 
          
           // tested with R function t.test(x, y, paired = TRUE, alternative = "two.sided")
            let vectorX = [|1.;2.;4.;8.|]
            let vectorY = [|10.;23.;11;9.|]
            let expectedPval = 0.10836944173355316
            let expectedStatistic = 2.26554660552391818
            let twoSamplePaired = Testing.TTest.twoSamplePaired vectorX vectorY
            Expect.floatClose Accuracy.high twoSamplePaired.PValue expectedPval "pValue should be equal."
            Expect.floatClose Accuracy.high twoSamplePaired.Statistic expectedStatistic "t statistic should be equal."
            Expect.equal twoSamplePaired.DegreesOfFreedom 3. "df should be equal."
            
            let vectorZ = [|-5.;-9.;0.;-8.|]
            let twoSamplePaired2 = Testing.TTest.twoSamplePaired vectorX vectorZ
            let expectedPval1 = 0.041226646225439562
            let expectedStatistic1 = -3.44031028692427698
            Expect.floatClose Accuracy.high twoSamplePaired2.PValue expectedPval1 "pValue should be equal."
            Expect.floatClose Accuracy.high twoSamplePaired2.Statistic expectedStatistic1 "t statistic should be equal."

            let vectorNan = [|nan;10.;23.;11.|]
            let twoSamplePaired3 = Testing.TTest.twoSamplePaired vectorX vectorNan
            Expect.isTrue (nan.Equals(twoSamplePaired3.PValue)) "pValue should be nan."
            Expect.isTrue (nan.Equals(twoSamplePaired3.Statistic)) "t statistic should be nan."
            
            let vectorBefore = [|10.;4.;15.;12.;6.|]
            let twoSamplePaired4() = Testing.TTest.twoSamplePaired  vectorBefore vectorX |>ignore
            Expect.throws twoSamplePaired4 "Vectors of different length"

           // test if exception Test works
           // Expect.throws (fun _ -> Testing.TTest.twoSamplePaired vectorX vectorY|>ignore) "Vetors should have equal length"
           
            let vectorWithInfinity = [|infinity;4.;15.;12.|]
            let twoSamplePairedInfinity = Testing.TTest.twoSamplePaired vectorWithInfinity vectorX
            Expect.isTrue (nan.Equals(twoSamplePairedInfinity.PValue)) "pValue should be nan."
            Expect.isTrue (nan.Equals(twoSamplePairedInfinity.Statistic)) "t statistic should be nan."

            let vectorWithNegativeInfinity =  [|infinity;4.;15.;12.|]
            let twoSampleNegativeInfinity = Testing.TTest.twoSamplePaired vectorWithNegativeInfinity vectorX
            Expect.isTrue (nan.Equals(twoSampleNegativeInfinity.PValue)) "pValue should be nan."
            Expect.isTrue (nan.Equals(twoSampleNegativeInfinity.Statistic)) "t statistic should be nan."

            let vectorNull = [|0.;0.;0.;0.|]
            let twoSamplePairedWithNullVector = Testing.TTest.twoSamplePaired vectorNull vectorY
            let expectedPval2 = 0.0272
            let expectedStatistic2 = 4.0451
            Expect.floatClose Accuracy.low twoSamplePairedWithNullVector.PValue expectedPval2 "pValue should be equal."
            Expect.floatClose Accuracy.low twoSamplePairedWithNullVector.Statistic expectedStatistic2 "t statistic should be equal."
        ]
        
     
[<Tests>]
let fTestTests = 
    // F-Test validated against res.ftest <- var.test(samplea, sampleb, alternative = "two.sided") RStudio 2022.02.3+492 "Prairie Trillium" Release (1db809b8323ba0a87c148d16eb84efe39a8e7785, 2022-05-20) for Windows

    let sampleFA =  [|5.0; 6.0; 5.8; 5.7|] 
    let sampleFB =  [|3.5; 3.7; 4.0; 3.3; 3.6|]
    let sampleNaN =  [|5.0; 6.0; 5.8; nan|]
    let sampleInf =  [|5.0; 6.0; 5.8; infinity|]
    let sampleNegInf =  [|5.0; 6.0; 5.8; -infinity|]
    let sampleties = [|5.0; 5.0; 5.8; 5.3|]

    // calculation of the F test 
    let fResult = FTest.testVariances sampleFA sampleFB
    let fResultTies = FTest.testVariances sampleFA sampleties 



    testList "Testing.FTest" [
        testCase "createFTest" <| fun () -> 
            Expect.floatClose Accuracy.low fResult.Statistic 2.82338 "statistics should be equal."
            Expect.floatClose Accuracy.low fResult.PValueTwoTailed 0.34172 "pValue should be equal."
        testCase "FTest NaN" <| fun () -> 
            Expect.throws (fun () -> ((FTest.testVariances sampleNaN sampleFB).Statistic) |> printf "%A") "FTest works with NaN"
        testCase "FTest infinities" <| fun () -> 
            Expect.throws (fun () -> ((FTest.testVariances sampleInf sampleFB).Statistic)|> printf "%A") "FTest works with divisions by infitity -> nan"
            Expect.throws (fun () -> ((FTest.testVariances sampleNegInf sampleFB).Statistic)|> printf "%A") "FTest works with divisions by -infitity -> nan"    
        testCase "FTest 2 ties" <| fun () -> 
            Expect.floatClose Accuracy.low fResultTies.Statistic 1.32748538 "statistics should be equal."
            Expect.floatClose Accuracy.low fResultTies.PValueTwoTailed 0.8214 "pValue should be equal."
        ]
[<Tests>]
let chiSquaredTests = 
    // ChiSquared https://www.graphpad.com/quickcalcs/chisquared2/
    // example from R
    // obs <- c(315, 101, 108, 32)
    // exp <- c(0.5625, 0.1875, 0.1875, 0.0625) 
    // chisq.test(obs, p = exp)
    let testCase1 =
        let expected = [312.75;104.25;104.25;34.75]
        let observed = [315.;101.;108.;32.]
        let df = expected.Length - 1
        ChiSquareTest.compute df expected observed

    //obs <- c(315, 101, 80, 32, 50)
    //exp <- c(0.5625, 0.1875, 0.0875, 0.0625,0.1) 
    //chisq.test(obs, p = exp)
    let testCase2 =
        let expected = [325.125;108.375;50.575;36.125;57.8]
        let observed = [315.;101.;80.;32.;50.] 
        let df = expected.Length - 1
        ChiSquareTest.compute df expected observed
        
    testList "Testing.ChiSquaredTest" [
        testCase "compute" <| fun () -> 
            Expect.isTrue (0.9254 = Math.Round(testCase1.PValueRight,4)) "pValue should be equal."
            Expect.isTrue (0.4700 = Math.Round(testCase1.Statistic,4)) "statistic should be equal."
            Expect.isTrue (0.000638 = Math.Round(testCase2.PValueRight,6)) "pValue should be equal."
            Expect.isTrue (19.461 = Math.Round(testCase2.Statistic,3)) "statistic should be equal."
            
    ]

[<Tests>]
let pearsonTests = 
    // examples from R
    // cor.test(x,y)
    let testCase1 =
        let seq1 = [44.4; 45.9; 41.9; 53.3; 44.7; 44.1; 50.7; 45.2; 60.1;]
        let seq2 = [ 2.6;  3.1;  2.5;  5.0;  3.6;  4.0;  5.2;  2.8;  3.8;]
        Correlation.testPearson seq1 seq2

    let testCase2 =
        let seq1 = [312.7; 104.2; 104.; 34.7]
        let seq2 = [315.5; 101.3; 108.; 32.2]
        Correlation.testPearson seq1 seq2
        
    testList "Testing.Correlation" [
        testCase "testPearson" <| fun () -> 
            Expect.isTrue (0.108173054 = Math.Round(testCase1.PValue,9)) "pValue should be equal"
            Expect.isTrue (0.000294627 = Math.Round(testCase2.PValue,9)) "pValue should be equal"
    ]


[<Tests>]
let holmTests =

    let largeSetnan        = 
        Frame.ReadCsv(location = @"data/holmHochberg_Input_nan.csv",hasHeaders = true,separators = ",").GetColumn<float>("pValues")
        |> Series.valuesAll
        |> Array.ofSeq
        |> Array.map (fun x -> if x.IsSome then x.Value else nan ) 
 
    let largeSet        = 
        largeSetnan |> Array.filter (fun x -> not (nan.Equals x))
        
    // calculated using stats "The R Stats Package" v.4.5.1, using the p.adjust with method parameter set to  "holm" 
    let largeSet_Expectednan        = 
        Frame.ReadCsv(location = @"data/fwer_holm_results.csv",hasHeaders = true,separators = ",").GetColumn<float>("pValues")
        |> Series.valuesAll
        |> Array.ofSeq
        |> Array.map (fun x -> if x.IsSome then x.Value else nan ) 

    let largeSet_Expected        = 
        largeSet_Expectednan 
        |> Array.filter (fun x -> not (nan.Equals x))

    testList "Testing.MultipleTesting.holmFWER" [
        
        testCase "testHolmLarge" (fun () -> 
            Expect.sequenceEqual 
                (largeSet |> MultipleTesting.holmFWER |> Seq.map (fun x -> Math.Round(x,9))) 
                (largeSet_Expected |> Seq.map (fun x -> Math.Round(x,9)))
                "adjusted pValues should be equal to the reference implementation."
        )

        testCase "testHolmLargeNaN" (fun () -> 
            TestExtensions.sequenceEqualRoundedNaN 9
                (largeSetnan |> MultipleTesting.holmFWER |> Seq.ofArray) 
                (largeSet_Expectednan |> Seq.ofArray)
                "adjusted pValues should be equal to the reference implementation."
        )

    ]

[<Tests>]
let hochbergTests =


    let largeSetnan        = 
        Frame.ReadCsv(location = @"data/holmHochberg_Input_nan.csv",hasHeaders = true,separators = ",").GetColumn<float>("pValues")
        |> Series.valuesAll
        |> Array.ofSeq
        |> Array.map (fun x -> if x.IsSome then x.Value else nan ) 
 
    let largeSet        = 
        largeSetnan |> Array.filter (fun x -> not (nan.Equals x))

    // calculated using stats "The R Stats Package" v.4.5.1, using the p.adjust with method parameter set to  "hochberg"
    let largeSet_Expectednan        = 
        Frame.ReadCsv(location = @"data/fwer_hochberg_results.csv",hasHeaders = true,separators = ",").GetColumn<float>("pValues")
        |> Series.valuesAll
        |> Array.ofSeq
        |> Array.map (fun x -> if x.IsSome then x.Value else nan ) 

    let largeSet_Expected        = 
        largeSet_Expectednan 
        |> Array.filter (fun x -> not (nan.Equals x))

    testList "Testing.MultipleTesting.hochbergFWER" [
        
        testCase "testHochbergLarge" (fun () -> 
            Expect.sequenceEqual 
                (largeSet |> MultipleTesting.hochbergFWER |> Seq.map (fun x -> Math.Round(x,9))) 
                (largeSet_Expected |> Seq.map (fun x -> Math.Round(x,9)))
                "adjusted pValues should be equal to the reference implementation."
        )

        testCase "testHochbergLargeNaN" (fun () -> 
            TestExtensions.sequenceEqualRoundedNaN 9
                (largeSetnan |> MultipleTesting.hochbergFWER |> Seq.ofArray) 
                (largeSet_Expectednan |> Seq.ofArray)
                "adjusted pValues should be equal to the reference implementation."
        )

    ]

[<Tests>]
let dunnSidakTests =
    // p = [0.01; 0.04; 0.03; 0.1; 0.5], m=5
    // expected values verified with Python: 1 - (1-p)^m
    let pValues = [| 0.01; 0.04; 0.03; 0.1; 0.5 |]
    let pNaN    = [| 0.01; nan;  0.03; 0.1; 0.5 |]

    testList "Testing.MultipleTesting.DunnSidak" [
        testCase "singleStepBasic" <| fun () ->
            let result = MultipleTesting.dunnSidakFWER pValues
            // 1 - (1-p)^5 for each p
            let expected = [| 0.049010; 0.184627; 0.141266; 0.40951; 0.96875 |]
            Array.iter2 (fun r e ->
                Expect.floatClose Accuracy.low r e "Single-step Šidák adjusted p-values should match."
            ) result expected
        testCase "singleStepNaN" <| fun () ->
            let result = MultipleTesting.dunnSidakFWER pNaN
            // m=4 valid values; 1 - (1-p)^4
            Expect.floatClose Accuracy.low result.[0] 0.039404 "p[0] with NaN should match."
            Expect.isTrue (Double.IsNaN result.[1]) "NaN position should remain NaN."
            Expect.floatClose Accuracy.low result.[2] 0.114707 "p[2] with NaN should match."
            Expect.floatClose Accuracy.low result.[3] 0.34390  "p[3] with NaN should match."
            Expect.floatClose Accuracy.low result.[4] 0.93750  "p[4] with NaN should match."
        testCase "holmSidakBasic" <| fun () ->
            let result = MultipleTesting.holmSidakFWER pValues
            // sorted p: [0.01, 0.03, 0.04, 0.1, 0.5]
            // raw: [1-0.99^5, 1-0.97^4, 1-0.96^3, 1-0.9^2, 1-0.5^1]
            // running max (already monotone): [0.04901, 0.11471, 0.11526, 0.19, 0.5]
            // back to original order: [0.04901, 0.11526, 0.11471, 0.19, 0.5]
            let expected = [| 0.049010; 0.115264; 0.114707; 0.19; 0.5 |]
            Array.iter2 (fun r e ->
                Expect.floatClose Accuracy.low r e "Holm–Šidák adjusted p-values should match."
            ) result expected
        testCase "holmSidakNaN" <| fun () ->
            let result = MultipleTesting.holmSidakFWER pNaN
            // m=4 valid values; sorted: [0.01, 0.03, 0.1, 0.5]
            // raw: [1-0.99^4, 1-0.97^3, 1-0.9^2, 1-0.5^1] = [0.039404, 0.087327, 0.19, 0.5]
            Expect.floatClose Accuracy.low result.[0] 0.039404 "p[0] Holm-Šidák NaN"
            Expect.isTrue (Double.IsNaN result.[1]) "NaN position should remain NaN."
            Expect.floatClose Accuracy.low result.[2] 0.087327 "p[2] Holm-Šidák NaN"
            Expect.floatClose Accuracy.low result.[3] 0.19000  "p[3] Holm-Šidák NaN"
            Expect.floatClose Accuracy.low result.[4] 0.50000  "p[4] Holm-Šidák NaN"
    ]

[<Tests>]
let benjaminiHochbergTests =


    let largeSetWithIdsnan = readCsv @"benjaminiHochberg_Input_nan.csv"
    let largeSetnan        = 
        largeSetWithIdsnan |> Array.map snd
        
    let largeSetWithIds        = 
        largeSetWithIdsnan |> Array.filter (fun (_,x) -> not (nan.Equals x))

    let largeSet        = 
        largeSetnan |> Array.filter (fun x -> not (nan.Equals x))

    let largeSetWithIds_Expectednan = readCsv @"benjaminiHochberg_AdjustedWithR_nan.csv"
    let largeSet_Expectednan        = largeSetWithIds_Expectednan |> Array.map snd
    
    let largeSetWithIds_Expected = 
        largeSetWithIds_Expectednan 
        |> Array.filter (fun (_,x) -> not (nan.Equals x))

    let largeSet_Expected        = 
        largeSet_Expectednan 
        |> Array.filter (fun x -> not (nan.Equals x))

    testList "Testing.MultipleTesting.BenjaminiHochberg" [
        
        testCase "testBHLarge" (fun () -> 
            Expect.sequenceEqual 
                (largeSet |> MultipleTesting.benjaminiHochbergFDR |> Seq.map (fun x -> Math.Round(x,9))) 
                (largeSet_Expected |> Seq.map (fun x -> Math.Round(x,9)))
                "adjusted pValues should be equal to the reference implementation."
        )

        testCase "testBHLargeNaN" (fun () -> 
            TestExtensions.sequenceEqualRoundedNaN 9
                (largeSetnan |> MultipleTesting.benjaminiHochbergFDR) 
                (largeSet_Expectednan |> Seq.ofArray)
                "adjusted pValues should be equal to the reference implementation."
        )

        testCase "testBHLargeBy" (fun () -> 
            Expect.sequenceEqual 
                (
                    largeSetWithIds 
                    |> MultipleTesting.benjaminiHochbergFDRBy id 
                    |> Seq.sortBy fst
                    |> Seq.map (fun (x,y) -> x, Math.Round(y,9))
                ) 
                (
                    largeSetWithIds_Expected 
                    |> Seq.sortBy fst
                    |> Seq.map (fun (x,y) -> x, Math.Round(y,9))
                )
                "adjusted pValues with keys should be equal to the reference implementation."
        )

        testCase "testBHLargeNaNBy" (fun () -> 
            Expect.sequenceEqual 
                (
                    [("A0",nan); ("A0",nan); yield! largeSetWithIds] 
                    |> MultipleTesting.benjaminiHochbergFDRBy id 
                    |> Seq.sortBy fst 
                    |> Seq.skip 2 
                    |> Seq.map (fun (x,y) -> x, Math.Round(y,9))
                ) 
                (
                    largeSetWithIds_Expected
                    |> Seq.sortBy fst
                    |> Seq.map (fun (x,y) -> x, Math.Round(y,9))
                )
                "adjusted pValues with keys should be equal to the reference implementation, ignoring nan."
        )
            
    ]



//[<Tests>] suspended because of long codecov run time
let qValuesTest =
  
    let largeSetWithIdsnan = readCsv @"benjaminiHochberg_Input_nan.csv"
    let largeSetnan = largeSetWithIdsnan |> Array.map snd
    
    let largeSetWithIds_Expectednan = readCsv @"qvaluesWithR_nan.csv"
    let largeSet_Expectednan = largeSetWithIds_Expectednan |> Array.map snd

    let largeSetWithIds_ExpectedRobustnan = readCsv @"qvaluesRobustWithR_nan.csv"
    let largeSet_ExpectedRobustnan = largeSetWithIds_ExpectedRobustnan |> Array.map snd

    
    let largeSet = 
        largeSetnan |> Array.filter (fun x -> not (nan.Equals x))
    
    let largeSet_Expected = 
        largeSet_Expectednan |> Array.filter (fun x -> not (nan.Equals x))

    let largeSet_ExpectedRobust = 
        largeSet_ExpectedRobustnan |> Array.filter (fun x -> not (nan.Equals x))


    testList "Testing.MultipleTesting.Qvalues" [
      
        testCase "ofPValues" (fun () -> 
            //tested against r qvalue package 2.26.0
            //pi0 estimation is in closed form in r package and therefore cannot be tested 
            //qvalue::qvalue(pvals,pi0=0.48345)
            let pi0 = 0.48345
            Expect.sequenceEqual 
                (largeSet |> MultipleTesting.Qvalues.ofPValues pi0 |> Seq.map (fun x -> Math.Round(x,9))) 
                (largeSet_Expected |> Seq.map (fun x -> Math.Round(x,9)))
                "qValues should be equal to the reference implementation."
        )

        testCase "ofPValues_nan" (fun () -> 
            //tested against r qvalue package 2.26.0
            //pi0 estimation is in closed form in r package and therefore cannot be tested 
            //qvalue::qvalue(pvals,pi0=0.48345)
            let pi0 = 0.48345
            TestExtensions.sequenceEqualRoundedNaN 9
                (largeSetnan |> MultipleTesting.Qvalues.ofPValues pi0 |> Seq.ofArray) 
                (largeSet_Expectednan |> Seq.ofArray)
                "qValues should be equal to the reference implementation."
        )
        
        testCase "ofPValuesRobust" (fun () -> 
            //tested against r qvalue package 2.26.0
            //pi0 estimation is in closed form in r package and therefore cannot be tested 
            //qvalue::qvalue(pvals,pi0=0.48345,pfdr=TRUE)
            let pi0 = 0.48345
            Expect.sequenceEqual 
                (largeSet |> MultipleTesting.Qvalues.ofPValuesRobust pi0 |> Seq.map (fun x -> Math.Round(x,9))) 
                (largeSet_ExpectedRobust |> Seq.map (fun x -> Math.Round(x,9)))
                "qValues Robust should be equal to the reference implementation."
        )

        testCase "ofPValuesRobust_nan" (fun () -> 
            //tested against r qvalue package 2.26.0
            //pi0 estimation is in closed form in r package and therefore cannot be tested 
            //qvalue::qvalue(pvals,pi0=0.48345,pfdr=TRUE)
            let pi0 = 0.48345
            TestExtensions.sequenceEqualRoundedNaN 9
                (largeSetnan |> MultipleTesting.Qvalues.ofPValuesRobust pi0 |> Seq.ofArray) 
                (largeSet_ExpectedRobustnan |> Seq.ofArray)
                "qValues Robust should be equal to the reference implementation."
        )

    ]


let createMetricTestInt metricName actual expected = testCase metricName (fun () -> Expect.equal actual expected (sprintf "Metric %s was calculated incorrectly." metricName))
let createMetricTestFloat accuracy metricName actual expected = testCase metricName (fun () -> Expect.floatClose accuracy actual expected (sprintf "Metric %s was calculated incorrectly." metricName))

[<Tests>]
let binaryConfusionMatrixTests =

    // binary classification
    //              | Predicted |
    //              |  P  |  N  |
    // | Actual | P |  3  |  1  |
    //          | N |  1  |  2  |
    
    let tp = 3.
    let tn = 2.
    let fp = 1.
    let fn = 1.
    
    let binaryCM = BinaryConfusionMatrix.create(int tp,int tn,int fp,int fn)
    let ofPredictions1 = BinaryConfusionMatrix.ofPredictions(1,[1;1;1;1;0;0;0],[1;1;1;0;1;0;0])
    let ofPredictions2 = BinaryConfusionMatrix.ofPredictions([true;true;true;true;false;false;false],[true;true;true;false;true;false;false])

    let expectedCM = {
        TP = 3
        TN = 2
        FP = 1
        FN = 1
    }

    testList "Testing.BinaryConfusionMatrix" [

        testCase "create" (fun _ -> Expect.equal binaryCM expectedCM "binary confusion matrix incorrectly created")
        testCase "ofPredictions1" (fun _ -> Expect.equal ofPredictions1 expectedCM "binary confusion matrix created incorrectly from observations with positive label")
        testCase "ofPredictions2" (fun _ -> Expect.equal ofPredictions2 expectedCM "binary confusion matrix created incorrectly from boolean observations")

        createMetricTestInt "TruePositives" binaryCM.TP 3
        createMetricTestInt "TrueNegatives" binaryCM.TN 2
        createMetricTestInt "FalsePositives" binaryCM.FP 1
        createMetricTestInt "FalseNegatives" binaryCM.FN 1
        
        testCase "thresholdMap implicit thresholds 1" (fun _ ->
            let actual = BinaryConfusionMatrix.thresholdMap(
                [true;true;true;true;false;false;false],
                [0.9 ;0.6 ;0.7 ; 0.2 ; 0.7; 0.3 ; 0.1]
            )
            let expected = [
                1.9, BinaryConfusionMatrix.create(0,3,0,4)
                0.9, BinaryConfusionMatrix.create(1,3,0,3)
                0.7, BinaryConfusionMatrix.create(2,2,1,2)
                0.6, BinaryConfusionMatrix.create(3,2,1,1)
                0.3, BinaryConfusionMatrix.create(3,1,2,1)
                0.2, BinaryConfusionMatrix.create(4,1,2,0)
                0.1, BinaryConfusionMatrix.create(4,0,3,0)
            ]
            Expect.sequenceEqual actual expected "binary threshold map not correctly created from binary predictions"
        )

        testCase "thresholdMap explicit thresholds 1" (fun _ ->
            let actual = BinaryConfusionMatrix.thresholdMap(
                [true;true;true;true;false;false;false],
                [0.9 ;0.6 ;0.7 ; 0.2 ; 0.7; 0.3 ; 0.1],
                [1.; 0.9; 0.8; 0.7; 0.6; 0.5; 0.4; 0.3; 0.2; 0.1; 0.]
            )
            let expected = [|
                1.9, BinaryConfusionMatrix.create(0,3,0,4)
                1.0, BinaryConfusionMatrix.create(0,3,0,4)
                0.9, BinaryConfusionMatrix.create(1,3,0,3)
                0.8, BinaryConfusionMatrix.create(1,3,0,3)
                0.7, BinaryConfusionMatrix.create(2,2,1,2)
                0.6, BinaryConfusionMatrix.create(3,2,1,1)
                0.5, BinaryConfusionMatrix.create(3,2,1,1)
                0.4, BinaryConfusionMatrix.create(3,2,1,1)
                0.3, BinaryConfusionMatrix.create(3,1,2,1)
                0.2, BinaryConfusionMatrix.create(4,1,2,0)
                0.1, BinaryConfusionMatrix.create(4,0,3,0)
                0. , BinaryConfusionMatrix.create(4,0,3,0)
            |]
            Expect.sequenceEqual actual expected "binary threshold map not correctly created from binary predictions"
        )

        testCase "thresholdMap: floating point error affects custom thresholds" (fun _ ->
            let actual = BinaryConfusionMatrix.thresholdMap(
                [true;true;true;true;false;false;false],
                [0.9 ;0.6 ;0.7 ; 0.2 ; 0.7; 0.3 ; 0.1],
                // these values are not exact due to floating point errors in addition. For example, the 0.7 is actually 0.70000000000000006661338147750939 which is > 0.7 and therefore produces an unexpected result
                [0. .. 0.1 .. 1.] |> List.rev 
            )
            let expected = [|
                1.9, BinaryConfusionMatrix.create(0,3,0,4)
                1.0, BinaryConfusionMatrix.create(0,3,0,4)
                0.9, BinaryConfusionMatrix.create(1,3,0,3)
                0.8, BinaryConfusionMatrix.create(1,3,0,3)
                0.7, BinaryConfusionMatrix.create(2,2,1,2)
                0.6, BinaryConfusionMatrix.create(3,2,1,1)
                0.5, BinaryConfusionMatrix.create(3,2,1,1)
                0.4, BinaryConfusionMatrix.create(3,2,1,1)
                0.3, BinaryConfusionMatrix.create(3,1,2,1)
                0.2, BinaryConfusionMatrix.create(4,1,2,0)
                0.1, BinaryConfusionMatrix.create(4,0,3,0)
                0. , BinaryConfusionMatrix.create(4,0,3,0)
            |]
            Expect.isFalse (actual = expected) "expected list comprehension threshold to produce slightly incorrent thresholds"
        )
    ]


[<Tests>]
let multiLabelConfusionMatrixTests =

    // multi label classification 
    //              |    Predicted  |
    //              |  A  |  B  | C |
    // | Actual | A |  3  |  1  | 1 |
    //          | B |  1  |  2  | 0 |
    //          | C |  2  |  0  | 4 |

    let c: Matrix<int> = 
        [
            [3; 1; 1]
            [1; 2; 0]
            [2; 0; 4]
        ]
        |> array2D
        |> Matrix.ofArray2D
    
    let expectedMLCM = 
        {
            Labels = [|"A"; "B"; "C"|]
            Confusion = 
                [
                    [3; 1; 1]
                    [1; 2; 0]
                    [2; 0; 4]
                ]
                |> array2D
                |> Matrix.ofArray2D
        }

    let multiLabelCM = MultiLabelConfusionMatrix.create([|"A";"B";"C"|], c)
    let ofPredictions = 
        MultiLabelConfusionMatrix.ofPredictions(
            [|"A"; "B"; "C"|],
            [|"A"; "A"; "A"; "A"; "A"; "B"; "B"; "B"; "C"; "C"; "C"; "C"; "C"; "C"|],
            [|"A"; "A"; "A"; "B"; "C"; "B"; "B"; "A"; "C"; "C"; "C"; "C"; "A"; "A"|]
        )
    let allVsAll = multiLabelCM |> MultiLabelConfusionMatrix.allVsAll
    let expectedAllVsAll =
        [
            "A", BinaryConfusionMatrix.create(3,6,3,2)
            "B", BinaryConfusionMatrix.create(2,10,1,1)
            "C", BinaryConfusionMatrix.create(4,7,1,2)
        ]

    testList "Testing.MultiLabelConfusionMatrix" [

        testCase "create" (fun _ -> Expect.equal multiLabelCM expectedMLCM "multi label confusion matrix incorrectly created")
        testCase "ofPredictions" (fun _ -> Expect.equal ofPredictions expectedMLCM "multi label confusion matrix created incorrectly from observations with positive label")

        testCase "oneVsAll1" (fun _ -> Expect.equal (snd expectedAllVsAll[0]) (multiLabelCM |> MultiLabelConfusionMatrix.oneVsRest "A") "all-vs-all binary confusion matrices incorrectly created from multi label confusion matrix")
        testCase "oneVsAll2" (fun _ -> Expect.equal (snd expectedAllVsAll[1]) (multiLabelCM |> MultiLabelConfusionMatrix.oneVsRest "B") "all-vs-all binary confusion matrices incorrectly created from multi label confusion matrix")
        testCase "oneVsAll3" (fun _ -> Expect.equal (snd expectedAllVsAll[2]) (multiLabelCM |> MultiLabelConfusionMatrix.oneVsRest "C") "all-vs-all binary confusion matrices incorrectly created from multi label confusion matrix")
        testCase "allVsAll" (fun _ -> Expect.sequenceEqual expectedAllVsAll allVsAll "all-vs-all binary confusion matrices incorrectly created from multi label confusion matrix")
        
    ]

[<Tests>]
let comparisonMetricsTests =

    testList "Testing.ComparisonMetrics" [

        // values calculated by formulas at https://en.wikipedia.org/wiki/Confusion_matrix
        let sensitivity = 0.75
        let specificity = 0.6666666667
        let precision = 0.75
        let negativePredictiveValue = 0.6666666667
        let missrate = 0.25
        let fallOut = 0.3333333333
        let falseDiscoveryRate = 0.25
        let falseOmissionRate = 0.3333333333
        let positiveLikelihoodRatio = sensitivity / fallOut
        let negativeLikelihoodRatio = missrate / specificity
        let prevalenceThreshold = sqrt(fallOut) / (sqrt(sensitivity) + sqrt(fallOut))
        let threatScore = 0.6
        let prevalence = 0.5714285714
        let accuracy = 0.7142857143
        let balancedAccuracy = (sensitivity + specificity) / 2.
        let f1 = 0.75
        let phiCoefficient = 0.4166666667
        let fowlkesMallowsIndex = 0.75
        let informedness = 0.4166666667
        let markedness = 0.4166666667
        let diagnosticOddsRatio = positiveLikelihoodRatio / negativeLikelihoodRatio

        let tp = 3.
        let tn = 2.
        let fp = 1.
        let fn = 1.
        let p = 4.
        let n = 3.
        let samplesize = 7.
        let binaryCM = BinaryConfusionMatrix.create(3,2,1,1)
        let cm = ComparisonMetrics.create(binaryCM)

        testList "Metric calculation" [
            createMetricTestFloat Accuracy.veryHigh "Calculate Sensitivity" (ComparisonMetrics.calculateSensitivity tp p) sensitivity
            createMetricTestFloat Accuracy.veryHigh "Calculate Specificity" (ComparisonMetrics.calculateSpecificity tn n) specificity
            createMetricTestFloat Accuracy.veryHigh "Calculate Precision" (ComparisonMetrics.calculatePrecision tp fp) precision
            createMetricTestFloat Accuracy.veryHigh "Calculate NegativePredictiveValue" (ComparisonMetrics.calculateNegativePredictiveValue tn fn) negativePredictiveValue
            createMetricTestFloat Accuracy.veryHigh "Calculate Missrate" (ComparisonMetrics.calculateMissrate fn p) missrate 
            createMetricTestFloat Accuracy.veryHigh "Calculate FallOut" (ComparisonMetrics.calculateFallOut fp n) fallOut
            createMetricTestFloat Accuracy.veryHigh "Calculate FalseDiscoveryRate" (ComparisonMetrics.calculateFalseDiscoveryRate fp tp) falseDiscoveryRate
            createMetricTestFloat Accuracy.veryHigh "Calculate FalseOmissionRate" (ComparisonMetrics.calculateFalseOmissionRate fn tn) falseOmissionRate
            createMetricTestFloat Accuracy.veryHigh "Calculate PositiveLikelihoodRatio" (ComparisonMetrics.calculatePositiveLikelihoodRatio tp p fp n) positiveLikelihoodRatio
            createMetricTestFloat Accuracy.veryHigh "Calculate NegativeLikelihoodRatio" (ComparisonMetrics.calculateNegativeLikelihoodRatio fn p tn n) negativeLikelihoodRatio
            createMetricTestFloat Accuracy.veryHigh "Calculate PrevalenceThreshold" (ComparisonMetrics.calculatePrevalenceThreshold fp n tp p) prevalenceThreshold
            createMetricTestFloat Accuracy.veryHigh "Calculate ThreatScore" (ComparisonMetrics.calculateThreatScore tp fn fp) threatScore
            createMetricTestFloat Accuracy.veryHigh "Calculate Prevalence" (ComparisonMetrics.calculatePrevalence p samplesize) prevalence
            createMetricTestFloat Accuracy.veryHigh "Calculate Accuracy" (ComparisonMetrics.calculateAccuracy tp tn samplesize) accuracy
            createMetricTestFloat Accuracy.veryHigh "Calculate BalancedAccuracy" (ComparisonMetrics.calculateBalancedAccuracy tp p tn n) balancedAccuracy
            createMetricTestFloat Accuracy.veryHigh "Calculate F1" (ComparisonMetrics.calculateF1 tp fp fn) f1
            createMetricTestFloat Accuracy.veryHigh "Calculate PhiCoefficient" (ComparisonMetrics.calculatePhiCoefficient tp tn fp fn) phiCoefficient
            createMetricTestFloat Accuracy.veryHigh "Calculate FowlkesMallowsIndex" (ComparisonMetrics.calculateFowlkesMallowsIndex tp fp p) fowlkesMallowsIndex
            createMetricTestFloat Accuracy.veryHigh "Calculate Informedness" (ComparisonMetrics.calculateInformedness tp p tn n) informedness
            createMetricTestFloat Accuracy.veryHigh "Calculate Markedness" (ComparisonMetrics.calculateMarkedness tp fp tn fn) markedness
            createMetricTestFloat Accuracy.veryHigh "Calculate DiagnosticOddsRatio" (ComparisonMetrics.calculateDiagnosticOddsRatio tp tn fp fn p n) diagnosticOddsRatio
        ]
        testList "Binary predictions" [

            createMetricTestInt "TruePositives" cm.TP 3
            createMetricTestInt "TrueNegatives" cm.TN 2
            createMetricTestInt "FalsePositives" cm.FP 1
            createMetricTestInt "FalseNegatives" cm.FN 1
            createMetricTestInt "Positves" cm.P 4
            createMetricTestInt "Negatives" cm.N 3
            createMetricTestInt "Total" cm.SampleSize 7

            createMetricTestFloat Accuracy.veryHigh "Sensitivity" cm.Sensitivity sensitivity
            createMetricTestFloat Accuracy.veryHigh "Specificity" cm.Specificity specificity
            createMetricTestFloat Accuracy.veryHigh "Precision" cm.Precision precision
            createMetricTestFloat Accuracy.veryHigh "NegativePredictiveValue" cm.NegativePredictiveValue negativePredictiveValue
            createMetricTestFloat Accuracy.veryHigh "Missrate" cm.Missrate missrate 
            createMetricTestFloat Accuracy.veryHigh "FallOut" cm.FallOut fallOut
            createMetricTestFloat Accuracy.veryHigh "FalseDiscoveryRate" cm.FalseDiscoveryRate falseDiscoveryRate
            createMetricTestFloat Accuracy.veryHigh "FalseOmissionRate" cm.FalseOmissionRate falseOmissionRate
            createMetricTestFloat Accuracy.veryHigh "PositiveLikelihoodRatio" cm.PositiveLikelihoodRatio positiveLikelihoodRatio
            createMetricTestFloat Accuracy.veryHigh "NegativeLikelihoodRatio" cm.NegativeLikelihoodRatio negativeLikelihoodRatio
            createMetricTestFloat Accuracy.veryHigh "PrevalenceThreshold" cm.PrevalenceThreshold prevalenceThreshold
            createMetricTestFloat Accuracy.veryHigh "ThreatScore" cm.ThreatScore threatScore
            createMetricTestFloat Accuracy.veryHigh "Prevalence" cm.Prevalence prevalence
            createMetricTestFloat Accuracy.veryHigh "Accuracy" cm.Accuracy accuracy
            createMetricTestFloat Accuracy.veryHigh "BalancedAccuracy" cm.BalancedAccuracy balancedAccuracy
            createMetricTestFloat Accuracy.veryHigh "F1" cm.F1 f1
            createMetricTestFloat Accuracy.veryHigh "PhiCoefficient" cm.PhiCoefficient phiCoefficient
            createMetricTestFloat Accuracy.veryHigh "FowlkesMallowsIndex" cm.FowlkesMallowsIndex fowlkesMallowsIndex
            createMetricTestFloat Accuracy.veryHigh "Informedness" cm.Informedness informedness
            createMetricTestFloat Accuracy.veryHigh "Markedness" cm.Markedness markedness
            createMetricTestFloat Accuracy.veryHigh "DiagnosticOddsRatio" cm.DiagnosticOddsRatio diagnosticOddsRatio

        ]
        testList "Multi-label predictions" [
            let c: Matrix<int> = 
                [
                    [3; 1; 1]
                    [1; 2; 0]
                    [2; 0; 4]
                ]
                |> array2D
                |> Matrix.ofArray2D

            let multiLabelCM  = MultiLabelConfusionMatrix.create([|"A";"B";"C"|], c)

            let expectedAvsRest = BinaryConfusionMatrix.create(3,6,3,2)
            let expectedBvsRest = BinaryConfusionMatrix.create(2,10,1,1)
            let expectedCvsRest = BinaryConfusionMatrix.create(4,7,1,2)

            let expectedMicroAverage = ComparisonMetrics.create(BinaryConfusionMatrix.create(9,23,5,5))

            let expectedMacroAverage =
                [
                    expectedAvsRest
                    expectedBvsRest
                    expectedCvsRest
                ]
                |> List.map (fun x -> ComparisonMetrics.create(x))
                |> fun metrics ->
                    ComparisonMetrics.create(
                        ((metrics[0].P                       + metrics[1].P                      + metrics[2].P                      ) / 3.),
                        ((metrics[0].N                       + metrics[1].N                      + metrics[2].N                      ) / 3.),
                        ((metrics[0].SampleSize              + metrics[1].SampleSize             + metrics[2].SampleSize             ) / 3.),
                        ((metrics[0].TP                      + metrics[1].TP                     + metrics[2].TP                     ) / 3.),
                        ((metrics[0].TN                      + metrics[1].TN                     + metrics[2].TN                     ) / 3.),
                        ((metrics[0].FP                      + metrics[1].FP                     + metrics[2].FP                     ) / 3.),
                        ((metrics[0].FN                      + metrics[1].FN                     + metrics[2].FN                     ) / 3.),
                        ((metrics[0].Sensitivity             + metrics[1].Sensitivity            + metrics[2].Sensitivity            ) / 3.),
                        ((metrics[0].Specificity             + metrics[1].Specificity            + metrics[2].Specificity            ) / 3.),
                        ((metrics[0].Precision               + metrics[1].Precision              + metrics[2].Precision              ) / 3.),
                        ((metrics[0].NegativePredictiveValue + metrics[1].NegativePredictiveValue+ metrics[2].NegativePredictiveValue) / 3.),
                        ((metrics[0].Missrate                + metrics[1].Missrate               + metrics[2].Missrate               ) / 3.),
                        ((metrics[0].FallOut                 + metrics[1].FallOut                + metrics[2].FallOut                ) / 3.),
                        ((metrics[0].FalseDiscoveryRate      + metrics[1].FalseDiscoveryRate     + metrics[2].FalseDiscoveryRate     ) / 3.),
                        ((metrics[0].FalseOmissionRate       + metrics[1].FalseOmissionRate      + metrics[2].FalseOmissionRate      ) / 3.),
                        ((metrics[0].PositiveLikelihoodRatio + metrics[1].PositiveLikelihoodRatio+ metrics[2].PositiveLikelihoodRatio) / 3.),
                        ((metrics[0].NegativeLikelihoodRatio + metrics[1].NegativeLikelihoodRatio+ metrics[2].NegativeLikelihoodRatio) / 3.),
                        ((metrics[0].PrevalenceThreshold     + metrics[1].PrevalenceThreshold    + metrics[2].PrevalenceThreshold    ) / 3.),
                        ((metrics[0].ThreatScore             + metrics[1].ThreatScore            + metrics[2].ThreatScore            ) / 3.),
                        ((metrics[0].Prevalence              + metrics[1].Prevalence             + metrics[2].Prevalence             ) / 3.),
                        ((metrics[0].Accuracy                + metrics[1].Accuracy               + metrics[2].Accuracy               ) / 3.),
                        ((metrics[0].BalancedAccuracy        + metrics[1].BalancedAccuracy       + metrics[2].BalancedAccuracy       ) / 3.),
                        ((metrics[0].F1                      + metrics[1].F1                     + metrics[2].F1                     ) / 3.),
                        ((metrics[0].PhiCoefficient          + metrics[1].PhiCoefficient         + metrics[2].PhiCoefficient         ) / 3.),
                        ((metrics[0].FowlkesMallowsIndex     + metrics[1].FowlkesMallowsIndex    + metrics[2].FowlkesMallowsIndex    ) / 3.),
                        ((metrics[0].Informedness            + metrics[1].Informedness           + metrics[2].Informedness           ) / 3.),
                        ((metrics[0].Markedness              + metrics[1].Markedness             + metrics[2].Markedness             ) / 3.),
                        ((metrics[0].DiagnosticOddsRatio     + metrics[1].DiagnosticOddsRatio    + metrics[2].DiagnosticOddsRatio    ) / 3.)
                    )


            let cmMicroAverage1 = ComparisonMetrics.microAverage multiLabelCM
            let cmMacroAverage1 = ComparisonMetrics.macroAverage multiLabelCM
            let cmMicroAverage2 = multiLabelCM |> MultiLabelConfusionMatrix.allVsAll |> Array.map snd |> ComparisonMetrics.microAverage
            let cmMacroAverage2 = multiLabelCM |> MultiLabelConfusionMatrix.allVsAll |> Array.map (snd >> ComparisonMetrics.create) |> ComparisonMetrics.macroAverage
            
            createMetricTestFloat Accuracy.veryHigh "microAverage: Sensitivity 1" cmMicroAverage1.Sensitivity expectedMicroAverage.Sensitivity
            createMetricTestFloat Accuracy.veryHigh "microAverage: Specificity 1" cmMicroAverage1.Specificity expectedMicroAverage.Specificity
            createMetricTestFloat Accuracy.veryHigh "microAverage: Precision 1" cmMicroAverage1.Precision expectedMicroAverage.Precision
            createMetricTestFloat Accuracy.veryHigh "microAverage: NegativePredictiveValue 1" cmMicroAverage1.NegativePredictiveValue expectedMicroAverage.NegativePredictiveValue
            createMetricTestFloat Accuracy.veryHigh "microAverage: Missrate 1" cmMicroAverage1.Missrate expectedMicroAverage.Missrate 
            createMetricTestFloat Accuracy.veryHigh "microAverage: FallOut 1" cmMicroAverage1.FallOut expectedMicroAverage.FallOut
            createMetricTestFloat Accuracy.veryHigh "microAverage: FalseDiscoveryRate 1" cmMicroAverage1.FalseDiscoveryRate expectedMicroAverage.FalseDiscoveryRate
            createMetricTestFloat Accuracy.veryHigh "microAverage: FalseOmissionRate 1" cmMicroAverage1.FalseOmissionRate expectedMicroAverage.FalseOmissionRate
            createMetricTestFloat Accuracy.veryHigh "microAverage: PositiveLikelihoodRatio 1" cmMicroAverage1.PositiveLikelihoodRatio expectedMicroAverage.PositiveLikelihoodRatio
            createMetricTestFloat Accuracy.veryHigh "microAverage: NegativeLikelihoodRatio 1" cmMicroAverage1.NegativeLikelihoodRatio expectedMicroAverage.NegativeLikelihoodRatio
            createMetricTestFloat Accuracy.veryHigh "microAverage: PrevalenceThreshold 1" cmMicroAverage1.PrevalenceThreshold expectedMicroAverage.PrevalenceThreshold
            createMetricTestFloat Accuracy.veryHigh "microAverage: ThreatScore 1" cmMicroAverage1.ThreatScore expectedMicroAverage.ThreatScore
            createMetricTestFloat Accuracy.veryHigh "microAverage: Prevalence 1" cmMicroAverage1.Prevalence expectedMicroAverage.Prevalence
            createMetricTestFloat Accuracy.veryHigh "microAverage: Accuracy 1" cmMicroAverage1.Accuracy expectedMicroAverage.Accuracy
            createMetricTestFloat Accuracy.veryHigh "microAverage: BalancedAccuracy 1" cmMicroAverage1.BalancedAccuracy expectedMicroAverage.BalancedAccuracy
            createMetricTestFloat Accuracy.veryHigh "microAverage: F1 1" cmMicroAverage1.F1 expectedMicroAverage.F1
            createMetricTestFloat Accuracy.veryHigh "microAverage: PhiCoefficient 1" cmMicroAverage1.PhiCoefficient expectedMicroAverage.PhiCoefficient
            createMetricTestFloat Accuracy.veryHigh "microAverage: FowlkesMallowsIndex 1" cmMicroAverage1.FowlkesMallowsIndex expectedMicroAverage.FowlkesMallowsIndex
            createMetricTestFloat Accuracy.veryHigh "microAverage: Informedness 1" cmMicroAverage1.Informedness expectedMicroAverage.Informedness
            createMetricTestFloat Accuracy.veryHigh "microAverage: Markedness 1" cmMicroAverage1.Markedness expectedMicroAverage.Markedness
            createMetricTestFloat Accuracy.veryHigh "microAverage: DiagnosticOddsRatio 1" cmMicroAverage1.DiagnosticOddsRatio expectedMicroAverage.DiagnosticOddsRatio

            createMetricTestFloat Accuracy.veryHigh "microAverage: Sensitivity 2" cmMicroAverage2.Sensitivity expectedMicroAverage.Sensitivity
            createMetricTestFloat Accuracy.veryHigh "microAverage: Specificity 2" cmMicroAverage2.Specificity expectedMicroAverage.Specificity
            createMetricTestFloat Accuracy.veryHigh "microAverage: Precision 2" cmMicroAverage2.Precision expectedMicroAverage.Precision
            createMetricTestFloat Accuracy.veryHigh "microAverage: NegativePredictiveValue 2" cmMicroAverage2.NegativePredictiveValue expectedMicroAverage.NegativePredictiveValue
            createMetricTestFloat Accuracy.veryHigh "microAverage: Missrate 2" cmMicroAverage2.Missrate expectedMicroAverage.Missrate 
            createMetricTestFloat Accuracy.veryHigh "microAverage: FallOut 2" cmMicroAverage1.FallOut expectedMicroAverage.FallOut
            createMetricTestFloat Accuracy.veryHigh "microAverage: FalseDiscoveryRate 2" cmMicroAverage2.FalseDiscoveryRate expectedMicroAverage.FalseDiscoveryRate
            createMetricTestFloat Accuracy.veryHigh "microAverage: FalseOmissionRate 2" cmMicroAverage2.FalseOmissionRate expectedMicroAverage.FalseOmissionRate
            createMetricTestFloat Accuracy.veryHigh "microAverage: PositiveLikelihoodRatio 2" cmMicroAverage1.PositiveLikelihoodRatio expectedMicroAverage.PositiveLikelihoodRatio
            createMetricTestFloat Accuracy.veryHigh "microAverage: NegativeLikelihoodRatio 2" cmMicroAverage1.NegativeLikelihoodRatio expectedMicroAverage.NegativeLikelihoodRatio
            createMetricTestFloat Accuracy.veryHigh "microAverage: PrevalenceThreshold 2" cmMicroAverage1.PrevalenceThreshold expectedMicroAverage.PrevalenceThreshold
            createMetricTestFloat Accuracy.veryHigh "microAverage: ThreatScore 2" cmMicroAverage1.ThreatScore expectedMicroAverage.ThreatScore
            createMetricTestFloat Accuracy.veryHigh "microAverage: Prevalence 2" cmMicroAverage1.Prevalence expectedMicroAverage.Prevalence
            createMetricTestFloat Accuracy.veryHigh "microAverage: Accuracy 2" cmMicroAverage1.Accuracy expectedMicroAverage.Accuracy
            createMetricTestFloat Accuracy.veryHigh "microAverage: BalancedAccuracy 2" cmMicroAverage1.BalancedAccuracy expectedMicroAverage.BalancedAccuracy
            createMetricTestFloat Accuracy.veryHigh "microAverage: F1 2" cmMicroAverage1.F1 expectedMicroAverage.F1
            createMetricTestFloat Accuracy.veryHigh "microAverage: PhiCoefficient 2" cmMicroAverage1.PhiCoefficient expectedMicroAverage.PhiCoefficient
            createMetricTestFloat Accuracy.veryHigh "microAverage: FowlkesMallowsIndex 2" cmMicroAverage1.FowlkesMallowsIndex expectedMicroAverage.FowlkesMallowsIndex
            createMetricTestFloat Accuracy.veryHigh "microAverage: Informedness 2" cmMicroAverage1.Informedness expectedMicroAverage.Informedness
            createMetricTestFloat Accuracy.veryHigh "microAverage: Markedness 2" cmMicroAverage1.Markedness expectedMicroAverage.Markedness
            createMetricTestFloat Accuracy.veryHigh "microAverage: DiagnosticOddsRatio 2" cmMicroAverage1.DiagnosticOddsRatio expectedMicroAverage.DiagnosticOddsRatio
        
            createMetricTestFloat Accuracy.veryHigh "macroAverage: Sensitivity 1" cmMacroAverage1.Sensitivity expectedMacroAverage.Sensitivity
            createMetricTestFloat Accuracy.veryHigh "macroAverage: Specificity 1" cmMacroAverage1.Specificity expectedMacroAverage.Specificity
            createMetricTestFloat Accuracy.veryHigh "macroAverage: Precision 1" cmMacroAverage1.Precision expectedMacroAverage.Precision
            createMetricTestFloat Accuracy.veryHigh "macroAverage: NegativePredictiveValue 1" cmMacroAverage1.NegativePredictiveValue expectedMacroAverage.NegativePredictiveValue
            createMetricTestFloat Accuracy.veryHigh "macroAverage: Missrate 1" cmMacroAverage1.Missrate expectedMacroAverage.Missrate 
            createMetricTestFloat Accuracy.veryHigh "macroAverage: FallOut 1" cmMacroAverage1.FallOut expectedMacroAverage.FallOut
            createMetricTestFloat Accuracy.veryHigh "macroAverage: FalseDiscoveryRate 1" cmMacroAverage1.FalseDiscoveryRate expectedMacroAverage.FalseDiscoveryRate
            createMetricTestFloat Accuracy.veryHigh "macroAverage: FalseOmissionRate 1" cmMacroAverage1.FalseOmissionRate expectedMacroAverage.FalseOmissionRate
            createMetricTestFloat Accuracy.veryHigh "macroAverage: PositiveLikelihoodRatio 1" cmMacroAverage1.PositiveLikelihoodRatio expectedMacroAverage.PositiveLikelihoodRatio
            createMetricTestFloat Accuracy.veryHigh "macroAverage: NegativeLikelihoodRatio 1" cmMacroAverage1.NegativeLikelihoodRatio expectedMacroAverage.NegativeLikelihoodRatio
            createMetricTestFloat Accuracy.veryHigh "macroAverage: PrevalenceThreshold 1" cmMacroAverage1.PrevalenceThreshold expectedMacroAverage.PrevalenceThreshold
            createMetricTestFloat Accuracy.veryHigh "macroAverage: ThreatScore 1" cmMacroAverage1.ThreatScore expectedMacroAverage.ThreatScore
            createMetricTestFloat Accuracy.veryHigh "macroAverage: Prevalence 1" cmMacroAverage1.Prevalence expectedMacroAverage.Prevalence
            createMetricTestFloat Accuracy.veryHigh "macroAverage: Accuracy 1" cmMacroAverage1.Accuracy expectedMacroAverage.Accuracy
            createMetricTestFloat Accuracy.veryHigh "macroAverage: BalancedAccuracy 1" cmMacroAverage1.BalancedAccuracy expectedMacroAverage.BalancedAccuracy
            createMetricTestFloat Accuracy.veryHigh "macroAverage: F1 1" cmMacroAverage1.F1 expectedMacroAverage.F1
            createMetricTestFloat Accuracy.veryHigh "macroAverage: PhiCoefficient 1" cmMacroAverage1.PhiCoefficient expectedMacroAverage.PhiCoefficient
            createMetricTestFloat Accuracy.veryHigh "macroAverage: FowlkesMallowsIndex 1" cmMacroAverage1.FowlkesMallowsIndex expectedMacroAverage.FowlkesMallowsIndex
            createMetricTestFloat Accuracy.veryHigh "macroAverage: Informedness 1" cmMacroAverage1.Informedness expectedMacroAverage.Informedness
            createMetricTestFloat Accuracy.veryHigh "macroAverage: Markedness 1" cmMacroAverage1.Markedness expectedMacroAverage.Markedness
            createMetricTestFloat Accuracy.veryHigh "macroAverage: DiagnosticOddsRatio 1" cmMacroAverage1.DiagnosticOddsRatio expectedMacroAverage.DiagnosticOddsRatio

            createMetricTestFloat Accuracy.veryHigh "macroAverage: Sensitivity 2" cmMacroAverage2.Sensitivity expectedMacroAverage.Sensitivity
            createMetricTestFloat Accuracy.veryHigh "macroAverage: Specificity 2" cmMacroAverage2.Specificity expectedMacroAverage.Specificity
            createMetricTestFloat Accuracy.veryHigh "macroAverage: Precision 2" cmMacroAverage2.Precision expectedMacroAverage.Precision
            createMetricTestFloat Accuracy.veryHigh "macroAverage: NegativePredictiveValue 2" cmMacroAverage2.NegativePredictiveValue expectedMacroAverage.NegativePredictiveValue
            createMetricTestFloat Accuracy.veryHigh "macroAverage: Missrate 2" cmMacroAverage2.Missrate expectedMacroAverage.Missrate 
            createMetricTestFloat Accuracy.veryHigh "macroAverage: FallOut 2" cmMacroAverage2.FallOut expectedMacroAverage.FallOut
            createMetricTestFloat Accuracy.veryHigh "macroAverage: FalseDiscoveryRate 2" cmMacroAverage2.FalseDiscoveryRate expectedMacroAverage.FalseDiscoveryRate
            createMetricTestFloat Accuracy.veryHigh "macroAverage: FalseOmissionRate 2" cmMacroAverage2.FalseOmissionRate expectedMacroAverage.FalseOmissionRate
            createMetricTestFloat Accuracy.veryHigh "macroAverage: PositiveLikelihoodRatio 2" cmMacroAverage2.PositiveLikelihoodRatio expectedMacroAverage.PositiveLikelihoodRatio
            createMetricTestFloat Accuracy.veryHigh "macroAverage: NegativeLikelihoodRatio 2" cmMacroAverage2.NegativeLikelihoodRatio expectedMacroAverage.NegativeLikelihoodRatio
            createMetricTestFloat Accuracy.veryHigh "macroAverage: PrevalenceThreshold 2" cmMacroAverage2.PrevalenceThreshold expectedMacroAverage.PrevalenceThreshold
            createMetricTestFloat Accuracy.veryHigh "macroAverage: ThreatScore 2" cmMacroAverage2.ThreatScore expectedMacroAverage.ThreatScore
            createMetricTestFloat Accuracy.veryHigh "macroAverage: Prevalence 2" cmMacroAverage2.Prevalence expectedMacroAverage.Prevalence
            createMetricTestFloat Accuracy.veryHigh "macroAverage: Accuracy 2" cmMacroAverage2.Accuracy expectedMacroAverage.Accuracy
            createMetricTestFloat Accuracy.veryHigh "macroAverage: BalancedAccuracy 2" cmMacroAverage2.BalancedAccuracy expectedMacroAverage.BalancedAccuracy
            createMetricTestFloat Accuracy.veryHigh "macroAverage: F1 2" cmMacroAverage2.F1 expectedMacroAverage.F1
            createMetricTestFloat Accuracy.veryHigh "macroAverage: PhiCoefficient 2" cmMacroAverage2.PhiCoefficient expectedMacroAverage.PhiCoefficient
            createMetricTestFloat Accuracy.veryHigh "macroAverage: FowlkesMallowsIndex 2" cmMacroAverage2.FowlkesMallowsIndex expectedMacroAverage.FowlkesMallowsIndex
            createMetricTestFloat Accuracy.veryHigh "macroAverage: Informedness 2" cmMacroAverage2.Informedness expectedMacroAverage.Informedness
            createMetricTestFloat Accuracy.veryHigh "macroAverage: Markedness 2" cmMacroAverage2.Markedness expectedMacroAverage.Markedness
            createMetricTestFloat Accuracy.veryHigh "macroAverage: DiagnosticOddsRatio 2" cmMacroAverage2.DiagnosticOddsRatio expectedMacroAverage.DiagnosticOddsRatio
        ]
        testList "binary threshold map" [
            let actual = 
                ComparisonMetrics.binaryThresholdMap(
                    [true;true;true;true;false;false;false],
                    [0.9 ;0.6 ;0.7 ; 0.2 ; 0.7; 0.3 ; 0.1]
                )

            let expected = [|
                1.9, BinaryConfusionMatrix.create(0,3,0,4) |> ComparisonMetrics.create
                0.9, BinaryConfusionMatrix.create(1,3,0,3) |> ComparisonMetrics.create
                0.7, BinaryConfusionMatrix.create(2,2,1,2) |> ComparisonMetrics.create
                0.6, BinaryConfusionMatrix.create(3,2,1,1) |> ComparisonMetrics.create
                0.3, BinaryConfusionMatrix.create(3,1,2,1) |> ComparisonMetrics.create
                0.2, BinaryConfusionMatrix.create(4,1,2,0) |> ComparisonMetrics.create
                0.1, BinaryConfusionMatrix.create(4,0,3,0) |> ComparisonMetrics.create
            |]
            testCase "threshold 1-9" (fun _ -> TestExtensions.comparisonMetricsEqualRounded 3 (snd actual[0]) (snd expected[0]) "Incorrect metrics for threshold 1.9")
            testCase "threshold 0-9" (fun _ -> TestExtensions.comparisonMetricsEqualRounded 3 (snd actual[1]) (snd expected[1]) "Incorrect metrics for threshold 0.9")
            testCase "threshold 0-7" (fun _ -> TestExtensions.comparisonMetricsEqualRounded 3 (snd actual[2]) (snd expected[2]) "Incorrect metrics for threshold 0.7")
            testCase "threshold 0-6" (fun _ -> TestExtensions.comparisonMetricsEqualRounded 3 (snd actual[3]) (snd expected[3]) "Incorrect metrics for threshold 0.6")
            testCase "threshold 0-3" (fun _ -> TestExtensions.comparisonMetricsEqualRounded 3 (snd actual[4]) (snd expected[4]) "Incorrect metrics for threshold 0.3")
            testCase "threshold 0-2" (fun _ -> TestExtensions.comparisonMetricsEqualRounded 3 (snd actual[5]) (snd expected[5]) "Incorrect metrics for threshold 0.2")
            testCase "threshold 0-1" (fun _ -> TestExtensions.comparisonMetricsEqualRounded 3 (snd actual[6]) (snd expected[6]) "Incorrect metrics for threshold 0.1")
        ]
        testList "multi-label threshold map" [
            let expectedMetricsMap =
                Map.ofList [
                    "A", [|
                        1.9, BinaryConfusionMatrix.create(0,9,0,5) |> ComparisonMetrics.create
                        0.9, BinaryConfusionMatrix.create(1,9,0,4) |> ComparisonMetrics.create
                        0.8, BinaryConfusionMatrix.create(2,9,0,3) |> ComparisonMetrics.create
                        0.7, BinaryConfusionMatrix.create(3,9,0,2) |> ComparisonMetrics.create
                        0.6, BinaryConfusionMatrix.create(3,9,0,2) |> ComparisonMetrics.create
                        0.5, BinaryConfusionMatrix.create(3,7,2,2) |> ComparisonMetrics.create
                        0.4, BinaryConfusionMatrix.create(4,6,3,1) |> ComparisonMetrics.create
                        0.3, BinaryConfusionMatrix.create(5,5,4,0) |> ComparisonMetrics.create
                        0.2, BinaryConfusionMatrix.create(5,4,5,0) |> ComparisonMetrics.create
                        0.1, BinaryConfusionMatrix.create(5,0,9,0) |> ComparisonMetrics.create
                        0.0, BinaryConfusionMatrix.create(5,0,9,0) |> ComparisonMetrics.create
                    |]
                    "B", [|
                        1.9, BinaryConfusionMatrix.create(0,11,0,3) |> ComparisonMetrics.create
                        0.9, BinaryConfusionMatrix.create(0,11,0,3) |> ComparisonMetrics.create
                        0.8, BinaryConfusionMatrix.create(1,11,0,2) |> ComparisonMetrics.create
                        0.7, BinaryConfusionMatrix.create(2,11,0,1) |> ComparisonMetrics.create
                        0.6, BinaryConfusionMatrix.create(2,11,0,1) |> ComparisonMetrics.create
                        0.5, BinaryConfusionMatrix.create(2,10,1,1) |> ComparisonMetrics.create
                        0.4, BinaryConfusionMatrix.create(3,10,1,0) |> ComparisonMetrics.create
                        0.3, BinaryConfusionMatrix.create(3,9,2,0) |> ComparisonMetrics.create
                        0.2, BinaryConfusionMatrix.create(3,9,2,0) |> ComparisonMetrics.create
                        0.1, BinaryConfusionMatrix.create(3,4,7,0) |> ComparisonMetrics.create
                        0.0, BinaryConfusionMatrix.create(3,0,11,0) |> ComparisonMetrics.create
                    |]
                    "C", [|
                        1.9, BinaryConfusionMatrix.create(0,8,0,6) |> ComparisonMetrics.create
                        0.9, BinaryConfusionMatrix.create(1,8,0,5) |> ComparisonMetrics.create
                        0.8, BinaryConfusionMatrix.create(3,8,0,3) |> ComparisonMetrics.create
                        0.7, BinaryConfusionMatrix.create(4,8,0,2) |> ComparisonMetrics.create
                        0.6, BinaryConfusionMatrix.create(4,7,1,2) |> ComparisonMetrics.create
                        0.5, BinaryConfusionMatrix.create(4,7,1,2) |> ComparisonMetrics.create
                        0.4, BinaryConfusionMatrix.create(5,7,1,1) |> ComparisonMetrics.create
                        0.3, BinaryConfusionMatrix.create(6,7,1,0) |> ComparisonMetrics.create
                        0.2, BinaryConfusionMatrix.create(6,5,3,0) |> ComparisonMetrics.create
                        0.1, BinaryConfusionMatrix.create(6,0,8,0) |> ComparisonMetrics.create
                        0.0, BinaryConfusionMatrix.create(6,0,8,0) |> ComparisonMetrics.create
                    |]
                ]   
            
            let actual = 
                ComparisonMetrics.multiLabelThresholdMap(
                    actual = [|"A"; "A"; "A"; "A"; "A"; "B"; "B"; "B"; "C"; "C"; "C"; "C"; "C"; "C"|],
                    predictions = [|
                        "A", [|0.8; 0.7; 0.9; 0.4; 0.3; 0.1; 0.2; 0.5; 0.1; 0.1; 0.1; 0.3; 0.5; 0.4|]
                        "B", [|0.0; 0.1; 0.0; 0.5; 0.1; 0.8; 0.7; 0.4; 0.0; 0.1; 0.1; 0.0; 0.1; 0.3|]
                        "C", [|0.2; 0.2; 0.1; 0.1; 0.6; 0.1; 0.1; 0.1; 0.9; 0.8; 0.8; 0.7; 0.4; 0.3|]
                    |]
                )

            testCase "A: threshold 1-9" (fun _ -> TestExtensions.comparisonMetricsEqualRounded 3 (snd (actual["A"][0])) (snd (expectedMetricsMap["A"][0])) "Incorrect metrics for threshold 1.9")
            testCase "A: threshold 0-9" (fun _ -> TestExtensions.comparisonMetricsEqualRounded 3 (snd (actual["A"][1])) (snd (expectedMetricsMap["A"][1])) "Incorrect metrics for threshold 0.9")
            testCase "A: threshold 0-8" (fun _ -> TestExtensions.comparisonMetricsEqualRounded 3 (snd (actual["A"][2])) (snd (expectedMetricsMap["A"][2])) "Incorrect metrics for threshold 0.8")
            testCase "A: threshold 0-7" (fun _ -> TestExtensions.comparisonMetricsEqualRounded 3 (snd (actual["A"][3])) (snd (expectedMetricsMap["A"][3])) "Incorrect metrics for threshold 0.7")
            testCase "A: threshold 0-6" (fun _ -> TestExtensions.comparisonMetricsEqualRounded 3 (snd (actual["A"][4])) (snd (expectedMetricsMap["A"][4])) "Incorrect metrics for threshold 0.6")
            testCase "A: threshold 0-5" (fun _ -> TestExtensions.comparisonMetricsEqualRounded 3 (snd (actual["A"][5])) (snd (expectedMetricsMap["A"][5])) "Incorrect metrics for threshold 0.5")
            testCase "A: threshold 0-4" (fun _ -> TestExtensions.comparisonMetricsEqualRounded 3 (snd (actual["A"][6])) (snd (expectedMetricsMap["A"][6])) "Incorrect metrics for threshold 0.4")
            testCase "A: threshold 0-3" (fun _ -> TestExtensions.comparisonMetricsEqualRounded 3 (snd (actual["A"][7])) (snd (expectedMetricsMap["A"][7])) "Incorrect metrics for threshold 0.3")
            testCase "A: threshold 0-2" (fun _ -> TestExtensions.comparisonMetricsEqualRounded 3 (snd (actual["A"][8])) (snd (expectedMetricsMap["A"][8])) "Incorrect metrics for threshold 0.2")
            testCase "A: threshold 0-1" (fun _ -> TestExtensions.comparisonMetricsEqualRounded 3 (snd (actual["A"][9])) (snd (expectedMetricsMap["A"][9])) "Incorrect metrics for threshold 0.1")
            testCase "A: threshold 0-0" (fun _ -> TestExtensions.comparisonMetricsEqualRounded 3 (snd (actual["A"][10])) (snd (expectedMetricsMap["A"][10])) "Incorrect metrics for threshold 0.0")

            testCase "B: threshold 1-9" (fun _ -> TestExtensions.comparisonMetricsEqualRounded 3 (snd (actual["B"][0])) (snd (expectedMetricsMap["B"][0])) "Incorrect metrics for threshold 1.9")
            testCase "B: threshold 0-9" (fun _ -> TestExtensions.comparisonMetricsEqualRounded 3 (snd (actual["B"][1])) (snd (expectedMetricsMap["B"][1])) "Incorrect metrics for threshold 0.9")
            testCase "B: threshold 0-8" (fun _ -> TestExtensions.comparisonMetricsEqualRounded 3 (snd (actual["B"][2])) (snd (expectedMetricsMap["B"][2])) "Incorrect metrics for threshold 0.8")
            testCase "B: threshold 0-7" (fun _ -> TestExtensions.comparisonMetricsEqualRounded 3 (snd (actual["B"][3])) (snd (expectedMetricsMap["B"][3])) "Incorrect metrics for threshold 0.7")
            testCase "B: threshold 0-6" (fun _ -> TestExtensions.comparisonMetricsEqualRounded 3 (snd (actual["B"][4])) (snd (expectedMetricsMap["B"][4])) "Incorrect metrics for threshold 0.6")
            testCase "B: threshold 0-5" (fun _ -> TestExtensions.comparisonMetricsEqualRounded 3 (snd (actual["B"][5])) (snd (expectedMetricsMap["B"][5])) "Incorrect metrics for threshold 0.5")
            testCase "B: threshold 0-4" (fun _ -> TestExtensions.comparisonMetricsEqualRounded 3 (snd (actual["B"][6])) (snd (expectedMetricsMap["B"][6])) "Incorrect metrics for threshold 0.4")
            testCase "B: threshold 0-3" (fun _ -> TestExtensions.comparisonMetricsEqualRounded 3 (snd (actual["B"][7])) (snd (expectedMetricsMap["B"][7])) "Incorrect metrics for threshold 0.3")
            testCase "B: threshold 0-2" (fun _ -> TestExtensions.comparisonMetricsEqualRounded 3 (snd (actual["B"][8])) (snd (expectedMetricsMap["B"][8])) "Incorrect metrics for threshold 0.2")
            testCase "B: threshold 0-1" (fun _ -> TestExtensions.comparisonMetricsEqualRounded 3 (snd (actual["B"][9])) (snd (expectedMetricsMap["B"][9])) "Incorrect metrics for threshold 0.1")
            testCase "B: threshold 0-0" (fun _ -> TestExtensions.comparisonMetricsEqualRounded 3 (snd (actual["B"][10])) (snd (expectedMetricsMap["B"][10])) "Incorrect metrics for threshold 0.0")

            testCase "C: threshold 1-9" (fun _ -> TestExtensions.comparisonMetricsEqualRounded 3 (snd (actual["C"][0])) (snd (expectedMetricsMap["C"][0])) "Incorrect metrics for threshold 1.9")
            testCase "C: threshold 0-9" (fun _ -> TestExtensions.comparisonMetricsEqualRounded 3 (snd (actual["C"][1])) (snd (expectedMetricsMap["C"][1])) "Incorrect metrics for threshold 0.9")
            testCase "C: threshold 0-8" (fun _ -> TestExtensions.comparisonMetricsEqualRounded 3 (snd (actual["C"][2])) (snd (expectedMetricsMap["C"][2])) "Incorrect metrics for threshold 0.8")
            testCase "C: threshold 0-7" (fun _ -> TestExtensions.comparisonMetricsEqualRounded 3 (snd (actual["C"][3])) (snd (expectedMetricsMap["C"][3])) "Incorrect metrics for threshold 0.7")
            testCase "C: threshold 0-6" (fun _ -> TestExtensions.comparisonMetricsEqualRounded 3 (snd (actual["C"][4])) (snd (expectedMetricsMap["C"][4])) "Incorrect metrics for threshold 0.6")
            testCase "C: threshold 0-5" (fun _ -> TestExtensions.comparisonMetricsEqualRounded 3 (snd (actual["C"][5])) (snd (expectedMetricsMap["C"][5])) "Incorrect metrics for threshold 0.5")
            testCase "C: threshold 0-4" (fun _ -> TestExtensions.comparisonMetricsEqualRounded 3 (snd (actual["C"][6])) (snd (expectedMetricsMap["C"][6])) "Incorrect metrics for threshold 0.4")
            testCase "C: threshold 0-3" (fun _ -> TestExtensions.comparisonMetricsEqualRounded 3 (snd (actual["C"][7])) (snd (expectedMetricsMap["C"][7])) "Incorrect metrics for threshold 0.3")
            testCase "C: threshold 0-2" (fun _ -> TestExtensions.comparisonMetricsEqualRounded 3 (snd (actual["C"][8])) (snd (expectedMetricsMap["C"][8])) "Incorrect metrics for threshold 0.2")
            testCase "C: threshold 0-1" (fun _ -> TestExtensions.comparisonMetricsEqualRounded 3 (snd (actual["C"][9])) (snd (expectedMetricsMap["C"][9])) "Incorrect metrics for threshold 0.1")
            testCase "C: threshold 0-0" (fun _ -> TestExtensions.comparisonMetricsEqualRounded 3 (snd (actual["C"][10])) (snd (expectedMetricsMap["C"][10])) "Incorrect metrics for threshold 0.0")
        ]
    ]
    
    
// declaration of a custom ID type for SAM tests      
type MyId = MyId of string

let inline isSameNumber a b =
    if Ops.isNan a && Ops.isNan b then true
    elif Ops.isNegInf a && Ops.isNegInf b then true
    elif Ops.isPosInf a && Ops.isPosInf b then true
    else a = b
    
type SAM.SAM<'id> when 'id :> IComparable with
    // make it convenient to swap the id from one type to another for SAM.SAM<'id>
    static member mapId (f: 'a -> 'b) (x: SAM.SAM<'a>) : SAM.SAM<'b> = SAM.SAM<_>.Create (f x.ID) x.Ri x.Si x.Statistics x.QValue x.Foldchange x.MeanA x.StDevA x.MeanB x.StDevB
    // make it convenient to compare the contents
    static member areSame a b isSameId =
        if not (isSameId a.ID b.ID) then false else
        [
            a.Ri        , b.Ri
            a.Si        , b.Si
            a.Statistics, b.Statistics
            a.QValue    , b.QValue
            a.Foldchange, b.Foldchange
            a.MeanA     , b.MeanA
            a.StDevA    , b.StDevA
            a.MeanB     , b.MeanB
            a.StDevB    , b.StDevB
        ]
        |> List.forall (fun (a,b) -> isSameNumber a b)
  
type SAM.SAMResult<'id> when 'id :> IComparable with
    // make it convenient to swap the id from one type to another for SAM.SAMResult<'id>
    static member mapId (f: 'a -> 'b) (x: SAM.SAMResult<'a>) : SAM.SAMResult<'b> =
        SAM.SAMResult<_>.Create
          x.S0
          x.Pi0
          x.Delta
          x.UpperCut
          x.LowerCut
          (x.PosSigBioitem       |> Array.map (SAM.SAM<_>.mapId f))
          (x.NegSigBioitem       |> Array.map (SAM.SAM<_>.mapId f))
          (x.NonSigBioitem       |> Array.map (SAM.SAM<_>.mapId f))
          (x.AveragePermutations |> Array.map (SAM.SAM<_>.mapId f))
          x.FDR
          x.SigCalledCount
          x.MedianFalsePositivesCount

    static member areSame a b isSameId =
        let individualValuesAreSame =
            [
                a.S0                       , b.S0
                a.Pi0                      , b.Pi0
                a.Delta                    , b.Delta
                a.UpperCut                 , b.UpperCut
                a.LowerCut                 , b.LowerCut
                a.FDR                      , b.FDR
                a.SigCalledCount           , b.SigCalledCount
                a.MedianFalsePositivesCount, b.MedianFalsePositivesCount
                
            ]
            |> List.forall (fun (a,b) -> isSameNumber a b)
        if not individualValuesAreSame then false else
        let permsAreSame =
            Array.zip a.AveragePermutations b.AveragePermutations
            |> Array.forall (fun (a,b) -> isSameId a.ID b.ID && SAM.SAM.areSame a b isSameId)
        if not permsAreSame then false else
        let nonSigAreSame =
            Array.zip a.NonSigBioitem b.NonSigBioitem
            |> Array.forall (fun (a,b) -> isSameId a.ID b.ID && SAM.SAM.areSame a b isSameId)
        if not nonSigAreSame then false else
        let negSigAreSame =
            Array.zip a.NegSigBioitem b.NegSigBioitem
            |> Array.forall (fun (a,b) -> isSameId a.ID b.ID && SAM.SAM.areSame a b isSameId)
        if not negSigAreSame then false else
        let posSigAreSame =
            Array.zip a.PosSigBioitem b.PosSigBioitem
            |> Array.forall (fun (a,b) -> isSameId a.ID b.ID && SAM.SAM.areSame a b isSameId)
        if not posSigAreSame then false else
        true
        
[<Tests>]
let SAMTests = 
     // data preparation 
     let df:Frame<string,string> = 
         Frame.ReadCsv(@"data/TestDataSAM.txt",hasHeaders=true,separators = "\t")
         |> Frame.indexRows "gene"
     // get Rowkeys as Array
     let rowheader :string[] = df.RowKeys |> Seq.toArray
     // chunk data into Arrays for sample1 and sample2 to compare them 
     let (sample1,sample2) :float[][] * float [][]=  
         df
         |> Frame.getRows
         |> Series.values
         //|> Seq.map (Series.values >> Seq.toArray >> Array.chunkBySize 4 >> fun x -> x.[0],x.[1])
         |> Seq.map (Series.values >> Seq.toArray >> Array.chunkBySize 3 >> fun x -> x.[0],x.[1])
         |> Array.ofSeq
         |> Array.unzip
     // map rowheader to samples, so its (GeneName,[float;float;float])
     
     let data1 = Array.zip rowheader sample1 
     let data2 = Array.zip  rowheader sample2
     let corrected1 = 
         let medCorrect = 
             medianCentering data1 
         Array.zip rowheader medCorrect
     let corrected2 = 
         let medCorrect = 
             medianCentering data2
         Array.zip rowheader medCorrect


     let result1 = FSharp.Stats.Testing.SAM.twoClassUnpaired 100 0.05 data1 data2 (System.Random(27))
     let result1NonStringId =
        let data1 = data1 |> Array.map (fun (id,rest) -> MyId id, rest)
        let data2 = data2 |> Array.map (fun (id,rest) -> MyId id, rest)
        FSharp.Stats.Testing.SAM.twoClassUnpaired 100 0.05 data1 data2 (System.Random(27))
         
     let result2 = twoClassUnpaired 100 0.05 data1 data2 (System.Random(1337))
     let result3 = twoClassUnpaired 100 0.05 corrected1 corrected2 (System.Random(1337))
     let result4 = twoClassUnpaired 100 0.05 data1 data1 (System.Random(27))

     testList "SAM Tests" [
         testCase "twoClassUnpaired Seed 27" <| fun () -> 
             Expect.floatClose Accuracy.low result1.S0 0.041419 "S0 should be equal."
             Expect.floatClose Accuracy.low result1.Pi0 0.388060 "Pi0 should be equal."   
             Expect.floatClose Accuracy.low result1.Delta 0.986042 "Delta should be equal."   
             Expect.floatClose Accuracy.low result1.UpperCut 1.301745 "Upper Cut should be equal."   
             Expect.floatClose Accuracy.low result1.LowerCut -1.632438 "Lower Cut should be equal."   
             Expect.floatClose Accuracy.low (result1.PosSigBioitem |> Array.length |> float ) 80. "PosSigBioitems should be equal."   
             Expect.floatClose Accuracy.low (result1.NegSigBioitem |> Array.length |> float ) 61. "NegSigBioitems should be equal."   
             Expect.floatClose Accuracy.low (result1.MedianFalsePositivesCount |> float ) 14. "medFP should be equal."   

         testCase "twoClassUnpaired Seed 1337" <| fun () -> 
             Expect.floatClose Accuracy.low result2.S0 0.041419 "S0 should be equal."
             Expect.floatClose Accuracy.low result2.Pi0 0.388060 "Pi0 should be equal."   
             Expect.floatClose Accuracy.low result2.Delta 0.997591 "Delta should be equal."   
             Expect.floatClose Accuracy.low result2.UpperCut 1.301745 "Upper Cut should be equal."   
             Expect.floatClose Accuracy.low result2.LowerCut -1.632438 "Lower Cut should be equal."   
             Expect.floatClose Accuracy.low (result2.PosSigBioitem |> Array.length |> float ) 80. "PosSigBioitems should be equal."   
             Expect.floatClose Accuracy.low (result2.NegSigBioitem |> Array.length |> float ) 61. "NegSigBioitems should be equal."   
             Expect.floatClose Accuracy.low (result2.MedianFalsePositivesCount |> float ) 17. "medFP should be equal."   

         testCase "twoClassUnpaired median centered Seed 1337" <| fun () -> 
             Expect.floatClose Accuracy.low result3.S0 0.026303 "S0 should be equal."
             Expect.floatClose Accuracy.low result3.Pi0 0.407960 "Pi0 should be equal."   
             Expect.floatClose Accuracy.low result3.Delta 1.036888 "Delta should be equal."   
             Expect.floatClose Accuracy.low result3.UpperCut 1.381420 "Upper Cut should be equal."   
             Expect.floatClose Accuracy.low result3.LowerCut -1.740646 "Lower Cut should be equal."   
             Expect.floatClose Accuracy.low (result3.PosSigBioitem |> Array.length |> float ) 78. "PosSigBioitems should be equal."   
             Expect.floatClose Accuracy.low (result3.NegSigBioitem |> Array.length |> float ) 61. "NegSigBioitems should be equal."   
             Expect.floatClose Accuracy.low (result3.MedianFalsePositivesCount |> float ) 17. "medFP should be equal."   
         testCase "twoClassUnpaired Seed 27 data similar " <| fun () -> 
             Expect.floatClose Accuracy.low result4.S0 0.005275473477 "S0 should be equal."
             Expect.floatClose Accuracy.low result4.Pi0 2.0 "Pi0 should be equal."   
             Expect.floatClose Accuracy.low result4.Delta 100000.0 "Delta should be equal."   
             Expect.floatClose Accuracy.low result4.UpperCut 0.0 "Upper Cut should be equal."   
             Expect.floatClose Accuracy.low result4.LowerCut 0.0 "Lower Cut should be equal."   
             Expect.floatClose Accuracy.low (result4.PosSigBioitem |> Array.length |> float ) 0. "PosSigBioitems should be equal."   
             Expect.floatClose Accuracy.low (result4.NegSigBioitem |> Array.length |> float ) 0. "NegSigBioitems should be equal."   
             Expect.floatClose Accuracy.low (result4.MedianFalsePositivesCount |> float ) 0. "medFP should be equal."   

         testCase "non string ID type doesn't affect results" <| fun () ->
             let result1 : SAMResult<MyId> = SAMResult.mapId MyId result1
             Expect.isTrue (SAMResult.areSame result1 result1NonStringId (=)) "the ID type being other than string shouldn't affect the results"
            
     ]


[<Tests>]
let anovaTests = 
    // tested against R and Python.statsmodels
    let a0 = [|9.105820655954515; 9.03375793269564;  8.915088867146936|]
    let b0 = [|9.218533797307847; 9.374165278011928; 9.09902059371791|]
    let a1 = [|9.054802062519045; 8.980015323427445; 9.066383039929153|]
    let b1 = [|8.819119236212753; 8.96115712165076;  9.042933513927682|]
    let a2 = [|8.80043362511518;  8.938413526335156; 8.886644395458491|]
    let b2 = [|8.614386401606065; 8.53829813235829;  8.611071068887409|]

    let oneWaySetup = [|a0;a1;a2|]
    let twoWaySetup = [|[|a0;b0|];[|a1;b1|];[|a2;b2|]|]
    
    let oneWayAnova         = Testing.Anova.oneWayAnova oneWaySetup
    let twoWayANOVAFixed    = twoWaySetup|>Anova.twoWayANOVA Anova.TwoWayAnovaModel.Fixed 
    let twoWayANOVAMixed    = twoWaySetup|>Anova.twoWayANOVA Anova.TwoWayAnovaModel.Mixed 
    let twoWayANOVARandom   = twoWaySetup|>Anova.twoWayANOVA Anova.TwoWayAnovaModel.Random 

    testList "Testing.Anova" [
        testCase "onewayAnova" <| fun () -> 
            Expect.equal                                oneWayAnova.Factor.DegreesOfFreedom                     2                       "Factor.DegreesOfFreedom deviates from expected value"       
            Expect.floatClose Accuracy.low              (Math.Round (oneWayAnova.Factor.MeanSquares,8))         0.022925                "Factor.MeanSquares deviates from expected value"
            Expect.floatClose Accuracy.low              (Math.Round (oneWayAnova.Factor.Significance,8))        0.072                   "Factor.Significance deviates from expected value"
            Expect.equal                                oneWayAnova.Factor.Source                               Anova.BetweenGroups     "Factor.Source deviates from expected value" 
            Expect.floatClose Accuracy.low              (Math.Round (oneWayAnova.Factor.Statistic,8))           4.211                   "Factor.Statistic deviates from expected value"  
            Expect.floatClose Accuracy.low              (Math.Round (oneWayAnova.Factor.SumOfSquares,8))        0.04585                 "Factor.SumOfSquares deviates from expected value"
            Expect.equal                                oneWayAnova.Error.DegreesOfFreedom                      6                       "Error.DegreesOfFreedom deviates from expected value"      
            Expect.floatClose Accuracy.low              (Math.Round (oneWayAnova.Error.MeanSquares,8))          0.005444                "Error.MeanSquares deviates from expected value" 
            Expect.equal                                oneWayAnova.Error.Source                                Anova.WithinGroups      "Error.Source deviates from expected value"  
            Expect.floatClose Accuracy.low              (Math.Round (oneWayAnova.Error.SumOfSquares,8))         0.03266                 "Error.SumOfSquares deviates from expected value"
            Expect.equal                                oneWayAnova.Total.DegreesOfFreedom                      8                       "Total.DegreesOfFreedom deviates from expected value"      
            Expect.floatClose Accuracy.low              (Math.Round (oneWayAnova.Total.MeanSquares,8))          0.009814                "Total.MeanSquares deviates from expected value" 
            Expect.floatClose Accuracy.low              (Math.Round (oneWayAnova.Total.SumOfSquares,8))         0.07851                 "Total.SumOfSquares deviates from expected value"
            // Expect.floatClose Accuracy.high 0.    oneWayAnova.Error.Significance        "Error.Significance deviates from expected value"
            // Expect.floatClose Accuracy.high 0.    oneWayAnova.Error.Statistic           "Error.Statistic deviates from expected value"
            // Expect.floatClose Accuracy.high 0.    oneWayAnova.Total.Significance        "Total.Significance deviates from expected value"
            // Expect.floatClose Accuracy.high 0.   oneWayAnova.Total.Source              "Total.Source deviates from expected value"  
            // Expect.floatClose Accuracy.high 0.    oneWayAnova.Total.Statistic           "Total.Statistic deviates from expected value"

        testCase "twoWayANOVAFixed" <| fun () -> 
            Expect.equal                        twoWayANOVAFixed.FactorFst.DegreesOfFreedom                     2                       "FactorFst.DegreesOfFreed deviated from expected value"                                   
            Expect.floatClose Accuracy.low      (Math.Round (twoWayANOVAFixed.FactorFst.MeanSquares,8))         0.2386                  "FactorFst.MeanSquares deviated from expected value"                                  
            Expect.floatClose Accuracy.low      (Math.Round (twoWayANOVAFixed.FactorFst.Significance,8))        0.00002717              "FactorFst.Significance deviated from expected value"                                 
            Expect.equal                        twoWayANOVAFixed.FactorFst.Source                               Anova.Residual          "FactorFst.Source deviated from expected value"                       
            Expect.floatClose Accuracy.low      (Math.Round (twoWayANOVAFixed.FactorFst.Statistic,8))           28.6043                 "FactorFst.Statistic deviated from expected value"                            
            Expect.floatClose Accuracy.low      (Math.Round (twoWayANOVAFixed.FactorFst.SumOfSquares,8))        0.4771                  "FactorFst.SumOfSquares deviated from expected value"                                 
            Expect.equal                        twoWayANOVAFixed.FactorSnd.DegreesOfFreedom                     1                       "FactorSnd.DegreesOfFreed deviated from expected value"                                   
            Expect.floatClose Accuracy.low      (Math.Round (twoWayANOVAFixed.FactorSnd.MeanSquares,8))         0.01404                 "FactorSnd.MeanSquares deviated from expected value"                                  
            Expect.floatClose Accuracy.low      (Math.Round (twoWayANOVAFixed.FactorSnd.Significance,8))        0.21898                 "FactorSnd.Significance deviated from expected value"                                 
            Expect.equal                        twoWayANOVAFixed.FactorSnd.Source                               Anova.Residual          "FactorSnd.Source deviated from expected value"                       
            Expect.floatClose Accuracy.low      (Math.Round (twoWayANOVAFixed.FactorSnd.Statistic,8))           1.6831                  "FactorSnd.Statistic deviated from expected value"                            
            Expect.floatClose Accuracy.low      (Math.Round (twoWayANOVAFixed.FactorSnd.SumOfSquares,8))        0.01404                 "FactorSnd.SumOfSquares deviated from expected value"                                 
            Expect.equal                        twoWayANOVAFixed.Interaction.DegreesOfFreedom                   2                       "Interaction.DegreesOfFre deviated from expected value"                                   
            Expect.floatClose Accuracy.low      (Math.Round (twoWayANOVAFixed.Interaction.MeanSquares,8))       0.09512                 "Interaction.MeanSquares deviated from expected value"                                    
            Expect.floatClose Accuracy.low      (Math.Round (twoWayANOVAFixed.Interaction.Significance,8))      0.001678                "Interaction.Significance deviated from expected value"                                   
            Expect.equal                        twoWayANOVAFixed.Interaction.Source                             Anova.BetweenGroups     "Interaction.Source deviated from expected value"                             
            Expect.floatClose Accuracy.low      (Math.Round (twoWayANOVAFixed.Interaction.Statistic,8))         11.4052                 "Interaction.Statistic deviated from expected value"                                  
            Expect.floatClose Accuracy.low      (Math.Round (twoWayANOVAFixed.Interaction.SumOfSquares,8))      0.1902                  "Interaction.SumOfSquares deviated from expected value"                                   
            Expect.equal                        twoWayANOVAFixed.Error.DegreesOfFreedom                         12                      "Error.DegreesOfFreedom deviated from expected value"                                 
            Expect.floatClose Accuracy.low      (Math.Round (twoWayANOVAFixed.Error.MeanSquares,8))             0.00834                 "Error.MeanSquares deviated from expected value"                              
            Expect.equal                        twoWayANOVAFixed.Total.DegreesOfFreedom                         17                      "Total.DegreesOfFreedom deviated from expected value"               
            // Expect.floatClose Accuracy.low      (Math.Round (twoWayANOVAFixed.Total.SumOfSquares,8))            0.7815                  "Total.SumOfSquares deviated from expected value"               
            // Expect.floatClose Accuracy.low      (Math.Round (twoWayANOVAFixed.Total.MeanSquares,8))             0.04597                 "Total.MeanSquares deviated from expected value"               

        testCase "twoWayANOVAMixed" <| fun () -> 
            Expect.equal                        twoWayANOVAMixed.FactorFst.DegreesOfFreedom                    2                       "FactorFst.DegreesOfFreed deviated from expected value"  
            Expect.floatClose Accuracy.low      (Math.Round (twoWayANOVAMixed.FactorFst.MeanSquares,8))        0.2386                  "FactorFst.MeanSquares deviated from expected value"     
            Expect.floatClose Accuracy.low      (Math.Round (twoWayANOVAMixed.FactorFst.Significance,8))       0.00002717              "FactorFst.Significance deviated from expected value"    
            Expect.equal                        twoWayANOVAMixed.FactorFst.Source                              Anova.Residual          "FactorFst.Source deviated from expected value"          
            Expect.floatClose Accuracy.low      (Math.Round (twoWayANOVAMixed.FactorFst.Statistic,8))          28.6043                 "FactorFst.Statistic deviated from expected value"       
            Expect.floatClose Accuracy.low      (Math.Round (twoWayANOVAMixed.FactorFst.SumOfSquares,8))       0.4771                  "FactorFst.SumOfSquares deviated from expected value"    
            Expect.equal                        twoWayANOVAMixed.FactorSnd.DegreesOfFreedom                     1                       "FactorSnd.DegreesOfFreed deviated from expected value"  
            Expect.floatClose Accuracy.low      (Math.Round (twoWayANOVAMixed.FactorSnd.MeanSquares,8))        0.01404                 "FactorSnd.MeanSquares deviated from expected value"     
            Expect.floatClose Accuracy.low      (Math.Round (twoWayANOVAMixed.FactorSnd.Significance,8))       0.7379                  "FactorSnd.Significance deviated from expected value"    
            Expect.equal                        twoWayANOVAMixed.FactorSnd.Source                              Anova.Residual          "FactorSnd.Source deviated from expected value"          
            Expect.floatClose Accuracy.low      (Math.Round (twoWayANOVAMixed.FactorSnd.Statistic,8))          0.1476                  "FactorSnd.Statistic deviated from expected value"       
            Expect.floatClose Accuracy.low      (Math.Round (twoWayANOVAMixed.FactorSnd.SumOfSquares,8))       0.01404                 "FactorSnd.SumOfSquares deviated from expected value"    
            Expect.equal                        twoWayANOVAMixed.Interaction.DegreesOfFreedom                  2                       "Interaction.DegreesOfFre deviated from expected value"  
            Expect.floatClose Accuracy.low      (Math.Round (twoWayANOVAMixed.Interaction.MeanSquares,8))      0.09512                 "Interaction.MeanSquares deviated from expected value"   
            Expect.floatClose Accuracy.low      (Math.Round (twoWayANOVAMixed.Interaction.Significance,8))     0.001678                "Interaction.Significance deviated from expected value"  
            Expect.equal                        twoWayANOVAMixed.Interaction.Source                            Anova.BetweenGroups     "Interaction.Source deviated from expected value"        
            Expect.floatClose Accuracy.low      (Math.Round (twoWayANOVAMixed.Interaction.Statistic,8))        11.4052                 "Interaction.Statistic deviated from expected value"     
            Expect.floatClose Accuracy.low      (Math.Round (twoWayANOVAMixed.Interaction.SumOfSquares,8))     0.1902                  "Interaction.SumOfSquares deviated from expected value"  
            Expect.equal                        twoWayANOVAMixed.Error.DegreesOfFreedom                        12                      "Error.DegreesOfFreedom deviated from expected value"    
            Expect.floatClose Accuracy.low      (Math.Round (twoWayANOVAMixed.Error.MeanSquares,8))            0.00834                 "Error.MeanSquares deviated from expected value"         
            Expect.floatClose Accuracy.low      (Math.Round (twoWayANOVAMixed.Error.SumOfSquares,8))           0.1001                  "Error.MeanSquares deviated from expected value"         
            Expect.equal                        twoWayANOVAMixed.Total.DegreesOfFreedom                        17                      "Total.DegreesOfFreedom deviated from expected value"               
            // Expect.floatClose Accuracy.low      (Math.Round (twoWayANOVAMixed.Total.SumOfSquares,8))           0.7815                  "Total.SumOfSquares deviated from expected value"               
            // Expect.floatClose Accuracy.low      (Math.Round (twoWayANOVAMixed.Total.MeanSquares,8))            0.04597                 "Total.MeanSquares deviated from expected value"               

        testCase "twoWayANOVARandom" <| fun () -> 
            Expect.equal                        twoWayANOVARandom.FactorFst.DegreesOfFreedom                    2                      "FactorFst.DegreesOfFreed deviated from expected value"  
            Expect.floatClose Accuracy.low      (Math.Round (twoWayANOVARandom.FactorFst.MeanSquares,8))        0.2386                 "FactorFst.MeanSquares deviated from expected value"     
            Expect.floatClose Accuracy.low      (Math.Round (twoWayANOVARandom.FactorFst.Significance,8))       0.2851                 "FactorFst.Significance deviated from expected value"    
            Expect.equal                        twoWayANOVARandom.FactorFst.Source                              Anova.Residual         "FactorFst.Source deviated from expected value"          
            Expect.floatClose Accuracy.low      (Math.Round (twoWayANOVARandom.FactorFst.Statistic,8))          2.508                  "FactorFst.Statistic deviated from expected value"       
            Expect.floatClose Accuracy.low      (Math.Round (twoWayANOVARandom.FactorFst.SumOfSquares,8))       0.4771                 "FactorFst.SumOfSquares deviated from expected value"    
            Expect.equal                        twoWayANOVARandom.FactorSnd.DegreesOfFreedom                    1                      "FactorSnd.DegreesOfFreed deviated from expected value"  
            Expect.floatClose Accuracy.low      (Math.Round (twoWayANOVARandom.FactorSnd.MeanSquares,8))        0.01404                "FactorSnd.MeanSquares deviated from expected value"     
            Expect.floatClose Accuracy.low      (Math.Round (twoWayANOVARandom.FactorSnd.Significance,8))       0.7379                 "FactorSnd.Significance deviated from expected value"    
            Expect.equal                        twoWayANOVARandom.FactorSnd.Source                              Anova.Residual         "FactorSnd.Source deviated from expected value"          
            Expect.floatClose Accuracy.low      (Math.Round (twoWayANOVARandom.FactorSnd.Statistic,8))          0.1476                 "FactorSnd.Statistic deviated from expected value"       
            Expect.floatClose Accuracy.low      (Math.Round (twoWayANOVARandom.FactorSnd.SumOfSquares,8))       0.01404                "FactorSnd.SumOfSquares deviated from expected value"    
            Expect.equal                        twoWayANOVARandom.Interaction.DegreesOfFreedom                  2                      "Interaction.DegreesOfFre deviated from expected value"  
            Expect.floatClose Accuracy.low      (Math.Round (twoWayANOVARandom.Interaction.MeanSquares,8))      0.09512                "Interaction.MeanSquares deviated from expected value"   
            Expect.floatClose Accuracy.low      (Math.Round (twoWayANOVARandom.Interaction.Significance,8))     0.001678               "Interaction.Significance deviated from expected value"  
            Expect.equal                        twoWayANOVARandom.Interaction.Source                            Anova.BetweenGroups    "Interaction.Source deviated from expected value"        
            Expect.floatClose Accuracy.low      (Math.Round (twoWayANOVARandom.Interaction.Statistic,8))        11.4052                "Interaction.Statistic deviated from expected value"     
            Expect.floatClose Accuracy.low      (Math.Round (twoWayANOVARandom.Interaction.SumOfSquares,8))     0.1902                 "Interaction.SumOfSquares deviated from expected value"  
            Expect.equal                        twoWayANOVARandom.Error.DegreesOfFreedom                        12                     "Error.DegreesOfFreedom deviated from expected value"    
            Expect.floatClose Accuracy.low      (Math.Round (twoWayANOVARandom.Error.MeanSquares,8))            0.00834                "Error.MeanSquares deviated from expected value"         
            Expect.floatClose Accuracy.low      (Math.Round (twoWayANOVARandom.Error.SumOfSquares,8))           0.1001                 "Error.SumOfSquares deviated from expected value"         
            Expect.equal                        twoWayANOVARandom.Total.DegreesOfFreedom                        17                     "Total.DegreesOfFreedom deviated from expected value"               
            // Expect.floatClose Accuracy.low      (Math.Round (twoWayANOVARandom.Total.SumOfSquares,8))           0.7815                 "Total.SumOfSquares deviated from expected value"               
            // Expect.floatClose Accuracy.low      (Math.Round (twoWayANOVARandom.Total.MeanSquares,8))            0.04597                "Total.MeanSquares deviated from expected value"               
            
        ]