How do I assign random values in an array of n elements satisfying following constraints?

Question

Suppose a is my array of n elements.

Constraints:

1. alpha <= a[i] <= beta.
2. 0 <= alpha <= beta <= 1.
3. Sum of all elements in a should be exactly 1.

Assuming such a array is always possible for given value of alpha and beta

This is what I am thinking:

Randomly assign values between alpha and beta. Find sum of array elements. If sum>1, subtract all elements with x such that constraints are satisfied. If sum<1, add all elements with x such that constraints are satisfied.

But I am finding it difficult to find out the newly mapped values.

Answer 1

This can be satisfied iff alpha × n ≤ 1 ≤ beta × n.

You can do this by picking n random numbers. Let r_i be the i^th random number, let R be the sum of the random numbers and S = (beta - alpha)/(1-n×alpha) × R. Set a[i] = alpha + r_i/S * (beta - alpha). As long as S is not less than the maximum random number then all elements of a are between alpha and beta and their sum is 1.

#include <iostream>
#include <cassert>
#include <random>
#include <vector>
#include <algorithm>
#include <iterator>

int main() {
    const double alpha = 0.05, beta = 0.15;
    const int n = 10;

    if (!(n * alpha <= 1.0 && 1.0 <= n * beta)) {
        std::cout << "unable to satisfy costraints.\n";
        return 0;
    }

    if (alpha == beta || std::fabs(1.0 - n * alpha) < 1e-6 || std::fabs(1.0 - n * beta) < 1e-6) {
        std::cout << "Solution for alpha = " << alpha << ", beta = " << beta << " n = " << n << ":\n";
        std::fill_n(std::ostream_iterator<double>(std::cout, " "), n, 1.0/n);
        std::cout << "\nSum: 1.0\n";
        return 0;
    }

    std::vector<int> r;
    double S;

    {
        std::random_device rd;
        std::seed_seq seed{rd(), rd(), rd(), rd(), rd(), rd(), rd(), rd()};
        std::mt19937 eng(seed);
        std::uniform_int_distribution<> dist(0, 1000);

        do {
            r.clear();
            std::generate_n(back_inserter(r), n, std::bind(dist, std::ref(eng)));

            int sum = std::accumulate(begin(r), end(r), 0);
            S = (beta - alpha) / (1 - n * alpha) * sum;
        } while (S < *std::max_element(begin(r), end(r)));

    }

    std::vector<double> a(n);
    std::transform(begin(r), end(r), begin(a), [&] (int ri) { return alpha + ri/S * (beta - alpha); });

    for (double ai : a) {
        assert(alpha <= ai && ai <= beta);
    }

    std::cout << "Solution for alpha = " << alpha << ", beta = " << beta << " n = " << n << ":\n";
    std::copy(begin(a), end(a), std::ostream_iterator<double>(std::cout, " "));
    std::cout << '\n';
    std::cout << "Sum: " << std::accumulate(begin(a), end(a), 0.0) << '\n';
}

Solution for alpha = 0.05, beta = 0.15, n = 10:
0.073923 0.117644 0.0834555 0.139368 0.101696 0.0846471 0.115261 0.0759395 0.0918882 0.116178
Sum: 1

---

You didn't specify the particular distribution you wanted the algorithm to provide, but I thought I'd point out that this algorithm does not necessarily produce every possible solution with equal probability; I believe some solutions are more likely to be produced than others.