mir-random statistical artifact

app.d

import std.algorithm;
import std.math;
import std.random;
import std.range : array, iota, take;
import std.stdio;
import mir.ndslice.sorting;
import mir.random.algorithm;
import mir.random.engine.pcg;
import mir.random.engine.splitmix;
import mir.random.variable;

enum n = 30; // sample size
enum numTrials = 5_000;

auto orderStat(double quantile)
{
    auto k = cast(int) (quantile * (n+1))-1;
    
    if(k < 0)
        return 0;
    else if(k >= n)
        return n-1;
    else
        return k;
}

auto cdf(double x)
{
    // trivial for uniform(0, 1)
    return x;

    // for exponential:
    // return 1-(exp(lambda * -x));
}

version = UsingMir;

void main()
{
    //auto rng = pcg32(mir.random.unpredictableSeed);
    //auto rng = SplitMixEngine!(staffordMix13, true)(1);
    auto ev = range!()(UniformVariable!double(0, 1));

    auto quantileFile = File("quantile.txt", "w");
    auto quantileDiffFile = File("quantile-diff.txt", "w");

    foreach(trialNum; 0 .. numTrials)
    {
        version(UsingMir)
        {
            auto sample = ev
                .take(n)
                .array;

            mir.ndslice.sorting.sort(sample);
        }
        else
        {
            auto sample = iota(n).map!(i => uniform!"[]"(0.0, 1.0)).array;
            std.algorithm.sort(sample);
        }

        auto requestedQuantile = uniform!"[]"(0.0, 1.0);
        auto k = orderStat(requestedQuantile);
        auto estimatedQuantileValue = sample[k];
        auto trueQuantile = cdf(estimatedQuantileValue);

        quantileFile.writefln("%s %s", requestedQuantile, trueQuantile);
        quantileDiffFile.writefln("%s %s", requestedQuantile,
            trueQuantile - requestedQuantile);
    }
}