I'm currently working on implementing a simple genetic algorithm in Julia, and I need a way to generate n random numbers between -1 and 1 with a sum of 0. Ideally, the distribution would be uniformly random, but it would still work as long as it's a reasonable approximation. A search on Google returns these two posts (1, 2) which are related but not exactly what I'm looking for.
I've reduced the problem to generating n numbers in the range [0, 1] with sum n/2, but I'm not sure where to go from there.
Simulate u_1, ..., u_n between -1 and 1. Set x1 = (u_1-u_n)/2, x2 = (u_2-u_1)/2, x3 = (u_3-u_2)/2, ... xn = (u_n-u_{n-1})/2. Then each number x is between -1 and 1, and the numbers xs sum to 0.
n = 10
u = 2*rand(n) .- 1
x = [(u[1]-u[n])/2; diff(u) ./ 2]
sum(x)
Here is the DRS algorithm. The call DRS(n, s, a, b) generates n random numbers between a and b that sum to s, and they are uniform.
# adapted from Roger Stafford's Matlab implementation
# https://nl.mathworks.com/matlabcentral/fileexchange/9700-random-vectors-with-fixed-sum
# Roger Stafford (2023). Random Vectors with Fixed Sum (https://www.mathworks.com/matlabcentral/fileexchange/9700-random-vectors-with-fixed-sum), MATLAB Central File Exchange. Retrieved March 23, 2023.
function DRS(n, s, a, b)
if s < n*a || s > n*b || a >= b
error("invalid values")
end
# Rescale to a unit cube: 0 <= x(i) <= 1
s = (s - n*a) / (b-a)
# Construct the transition probability table, t.
# t(i,j) will be utilized only in the region where j <= i + 1.
k = max(min(Int(floor(s)), n-1), 0) # Must have 0 <= k <= n-1
s = max(min(s, k+1), k) # Must have k <= s <= k+1
s1 = [s - i for i in k:-1:k-n+1] # s1 & s2 will never be negative
s2 = [i - s for i in k+n:-1:k+1]
w = zeros(n, n+1)
w[1, 2] = prevfloat(typemax(Float64)) # Scale for full 'double' range
t = zeros(n-1, n)
tiny = eps(Float64) # The smallest positive 'double'
for i = 2:n
tmp1 = w[i-1, 2:i+1] .* s1[1:i]/i
tmp2 = w[i-1, 1:i] .* s2[n-i+1:n]/i
w[i, 2:i+1] = tmp1 + tmp2
tmp3 = w[i, 2:i+1] .+ tiny # In case tmp1 & tmp2 are both 0,
tmp4 = s2[n-i+1:n] .> s1[1:i] # then t is 0 on left & 1 on right
t[i-1, 1:i] = (tmp2 ./ tmp3) .* tmp4 + (1 .- tmp1 ./ tmp3) .* (.!tmp4)
end
# Derive the polytope volume v from the appropriate element in the bottom row of w.
v = n^(3/2) * (w[n, k+2] / prevfloat(typemax(Float64))) * (b - a)^(n - 1)
# Now compute the matrix x.
x = zeros(n)
rt = rand(n - 1) # For random selection of simplex type
rs = rand(n - 1) # For random location within a simplex
j = k + 1 # For indexing in the t table
sm = 0
pr = 1 # Start with sum zero & product 1
for i = (n-1):(-1):1 # Work backwards in the t table
e = (rt[n-i] <= t[i, j]) # Use rt to choose a transition
sx = rs[n-i] ^ (1/i) # Use rs to compute next simplex coord.
sm = sm + (1 - sx) * pr * s / (i + 1) # Update sum
pr = sx * pr # Update product
x[n-i] = sm + pr * e # Calculate x using simplex coords.
s = s - e
j = j - e # Transition adjustment
end
x[n] = sm + pr * s # Compute the last x
# Randomly permute the order in x and rescale.
rp = rand(n) # Use rp to carry out a random permutation
p = sortperm(rp) #
return (b - a) * x[p] .+ a
end