Fit a quantity of items to a binomial distribution / bell curve

Question

This is a sort of discrete math integration question - I need to fit a fixed quantity of items to a binomial distribution or bell curve over a fixed period.

Say I a total of M boxes being shipped over T days, where n is the number boxes that arrive on day t. I need a way to calculate n(t) for each day t, so that

1. The Sum ( n(t) ) 0 -> t = M

2. t is an integer, and n(t) is an integer, and

3. the shape of n(t) matches as closely as possible to a bell curve.

---

Edit

In case anyone does think this an SO-worthy question, here's the Javascript I cobbled together from the pointers in Yves Daoust answer.

/*
See https://stackoverflow.com/questions/5259421/
*/
function normal(x, mean, stdDev) {
    return stdNormal(( x - mean ) / stdDev);
}

function stdNormal(z) {
    // Power series is not stable at these extreme tail scenarios
    if (z < -6) { return 0; }
    if (z >  6) { return 1; }

    let j, k ;
    
    let m = 1;        // m(k) == (2**k)/factorial(k)
    let b = z;        // b(k) == z ** (2*k + 1)
    let z2 = z * z;    // cache of z squared
    let z4 = z2 * z2;  // cache of z to the 4th
    let values = [];

    // Compute the power series in groups of two terms.
    // This reduces floating point errors because the series
    // alternates between positive and negative.
    for (k=0; k<100; k+=2) {
        const a = 2*k + 1;
        let item = b / (a*m);
        item *= (1 - (a*z2)/((a+1)*(a+2)));
        values.push(item);
        m *= (4*(k+1)*(k+2));
        b *= z4;
    }

    // Add the smallest terms to the total first that
    // way we minimize the floating point errors.
    let total = 0;
    for (k=49; k>=0; k--) {
        total += values[k];
    }

    // Multiply total by 1/sqrt(2*PI)
    // Then add 0.5 so that stdNormal(0) === 0.5
    return 0.5 + 0.3989422804014327 * total;
}

/*
Compute the cdf of the binomial distribution between 0 and T, times M.
Round all values to integers. Let m(t) the numbers so obtained.
Use n(t) = m(t+1) - m(t). In doing so, you ensure Σ n(t) = n(T) - n(0) = M.
Thanks to Yves Daoust
*/
function distributeItems(itemsToPlace = 100, steps = 7) {
    const mean = Math.floor((steps - 1) / 2);
    const stdDev = mean / 2.96; /// for 'standard' std dev ;-)
    const m = [0]; // cdf
    const n = [0]; // items
    for (var step = 1; step <= steps; step++) {
        m.push(Math.round(normal(step, mean, stdDev) * itemsToPlace,0));
        n.push( m[step] - m[step - 1] );
    }
    const interimCount = n.reduce(function (sum, elt) { return sum+elt; }, 0);
    const discrepancy = itemsToPlace - interimCount;
    if (discrepancy !==0) {
        n[n.length-1] += discrepancy;
    }
    return n;
}

const n = distributeItems(40,7);
console.log('Items: ',n, 'Total: ',n.reduce(function (sum, elt) { return sum+elt; }, 0))

// Result
// [ 0, 1, 5, 14, 14, 5, 1, 0 ] 40

Answer 1

Maybe not an optimal solution, but should not be far.

Compute the cdf of the binomial distribution between 0 and T, times M.

Round all values to integers. Let m(t) the numbers so obtained.

Use n(t) = m(t+1) - m(t). In doing so, you ensure Σ n(t) = n(T) - n(0) = M.

In the picture below, the blue curve is an exact binomial N=7, p=0.3, adjusted for M=20. The orange curve is obtained by the above procedure. As you can check, 2+5+6+4+2+1+0+0=20.