# Generating uniformly distributed random solutions of a linear equation

Given $$n+1$$ variables $$p_0, p_1, \ldots, p_n$$ defined over $$\mathbb{R}^{+}$$ so that $$\sum_{i=0}^np_i=1$$, and given a real number $$1, I want to generate random solutions of the equation so that every solution is equiprobable (or close enough to be equiprobable) $$0p_0 + 1p_1+\ldots+np_n=x$$

Since all variables are real numbers, the number of solutions is infinite, but I don't know how to algorithmically choose any of them randomly in an (at least close to) uniform way

NOTE: In case it's needed, suppose I have a method $$g(a, b)$$ that, when called, returns a random real number between $$a$$ and $$b$$.

Possible solution: (I strongly suspect this is a solution far from being uniform). Generate $$n$$ random real numbers $$r_i$$ defined in the interval 0 to 1. They will represent coefficients of $$x$$ so that $$ip_i=xr_i$$ (since $$0p_0=0$$ add nothing I don't need to generate $$r_0$$). That means $$\sum r_i$$ must be equal to $$1$$. To make that happen, generate all the $$r_i$$ randomly as said, but if $$\sum r_i=k\neq 1$$, multiply each $$r_i$$ by $$1/k$$.

Assign $$p_i=xr_i/i$$ and $$p_0=1-\sum_{i=1}^np_i$$. For the sum of all $$p_i$$ to be equal to $$1$$ is required that $$\frac{r_1}{1}+\frac{r_2}{2}+\ldots+\frac{r_n}{n}\leq\frac{1}{x}$$

but I don't know if that's guaranteed or not and, even if it's guaranteed, I suspect the generated solution might not be choosen uniformely because $$p_0$$ is deterministically calculated.

• This question might get you going in a useful direction? Feb 27, 2023 at 18:44
• @Eoin there's 2 equations and $n$ real variables (actually, in my case $n=7$ and $0<x<3$ but I wanted to present the problem in a more general form). As far as $n>2$ there should be infinite solutions right?
– ABu
Feb 27, 2023 at 18:49
• The notion of "uniformity" depends on the reference measure used on the resulting manifold. Feb 27, 2023 at 20:35
• You say "given $p_0,\dots,p_n$" and "given $x$", but if you are given all of those, then there is nothing to generate or sample. I suspect you mean "given $x$, you want to sample $p_0,\dots,p_n$" from a particular distribution. Is that correct? If so, I encourage you to edit the question accordingly.
– D.W.
May 9, 2023 at 3:37

The intersection of the $$n+1$$ simplex $$\{\mathbf{p}\in\mathbb R_+^{n+1};\ \mathbf{p}^\top\mathbf 1_{n+1}=1\}$$ when $$\mathbf 1_{n+1}=(1,\ldots1)^\top$$ and of the constrained hyperplane $$\{\mathbf{p}\in\mathbb R_+^{n+1};\ \mathbf{p}^\top\mathbf \iota_{n+1}=x\}$$ when $$\iota_{n+1}=(0,1,\ldots n)^\top$$ is within a $$(n−1)$$-dimensional affine space $$\{\mathbf{p}\in\mathbb R_+^{n+1};\ A\mathbf p=b\}\tag{1}$$ where $$A$$ is a $$(2,n+1)$$ matrix whose rows are orthonormal (wlog).

If $$\mathbf p^0$$ is a particular solution, i.e., a particular element of (1), the other members will be of the form $$\mathbf p^0+\eta$$ with $$A\eta=0$$, which can be expressed via an orthonormal basis of vectors satisfying $$A\eta=0$$. It is then sufficient to find an hypercube containing (1) by finding upper and lower bounds on the components of $$\eta$$ in (1), to generate uniformly points in that hypercube and accept simulations such that $$\mathbf p^0+\eta\in\mathbb R_+^{n+1}$$

• As $n$ grows, this approach becomes exponentially less efficient.
– whuber
May 9, 2023 at 17:03

Here's on possible approach, in R. It takes random samples from a simplex as starting points, and then uses optimisation to find nearby points that meet your criteria. I'm not sure how to verify that the resulting values are still "uniformly" distributed in the admissible parameter space.

find_solution = function(n, x){
# Set up our function
basis_w = t(as.matrix(0:n))
calculate_x = function(ps){
xhat = basis_w %*% as.matrix(ps)
xhat[1,1]
}

# Create a cost function to minimize.
# It should penalise values of p out of the range 0-1,
# ps that don't sum to 1, and values of x that don't match the target
cost_func = function(ps){
# Out of bounds loss
too_high = ifelse(ps > 1, ps - 1, 0) %>% sq() %>% sum()
too_low = ifelse(ps < 0,  -ps, 0) %>% sq() %>% sum()
# Doesn't sum to 1-loss
simplex_loss = sq(sum(ps) - 1)
# Wrong x value loss
xhat = calculate_x(ps)
err_loss = sq(xhat - x)
# Manually chosing tuning parameters
100*too_high + 100*too_low + simplex_loss + err_loss
}

# Randomly sample from a simplex for starting values
# (as array rather than matrix for compatibility with optim() )
starting_values = hitandrun::simplex.sample(n + 1, 1)$$samples %>% c() # Optimise result = optim(cost_func, par = starting_values) if(abs(result$$value > .001)) warning('Solution not found')
return(result$par) } find_solution(7, 1.5) # [1] 0.3474378316 0.3221905230 0.0008495489 0.2713962006 0.0002233641 0.0003346271 0.0436698092 0.0138942988 find_solution(7, 1.5) # [1] 3.889560e-01 2.736482e-01 1.178040e-01 7.962354e-02 3.027082e-05 8.923485e-02 4.931859e-02 1.381217e-03 find_solution(7, 1.5) # [1] 2.184026e-01 3.564897e-01 2.410103e-01 1.512842e-01 9.407019e-04 -1.540388e-05 7.533937e-06 2.939822e-02 find_solution(7, 1.5) # [1] 1.931718e-01 4.007192e-01 1.561489e-01 2.386094e-01 2.693653e-03 1.691173e-05 6.235019e-05 8.569269e-03 find_solution(7, 1.5) # [1] 0.4523700298 0.0139939226 0.3043223032 0.1165318696 0.0502990295 0.0478577440 0.0138911475 0.0005681597 find_solution(7, 1.5) # [1] 0.527091876 0.035665984 0.006260825 0.372083535 0.011015178 0.014363521 0.007900288 0.024702057  • If indeed you are projecting to nearest solutions -- I haven't checked your code -- then your results will not have the intended uniform distribution. The larger$n\$ is, the more they will deviate from uniformity. (This is another manifestation of the "curse of dimensionality.")
– whuber
Feb 27, 2023 at 20:24