I want to generate six Random numbers $a,b,c,d,e,f$ (all are real and positive). They are related by $a^2+b^2+c^2+d^2+e^2=1$ and $f$ is between 0 and 3. How can I generate those numbers?
2 Answers
$\begingroup$
$\endgroup$
Letting the 5-vector $ x = \{a^2, b^2, c^2, d^2, e^2 \} = \{x_1, x_2, x_3, x_4, x_5\}$ have DirichletDistribution[{1,1,1,1,1}], RandomVariate[DirichletDistribution[{1,1,1,1,1}]] gives a random 4-vector that satisfies $ 0 \leq x_i \leq 1$ and $\sum_{i=1}^4 x_i <1$. The 5th component of $x$ is determined by the condition $x_5= 1 - x_1 - x_2 - x_3 - x_4$. So appending $1 - \sum_{i=1}^4 x_i$ to $\{x_1, x_2, x_3, x_4\}$ returned by RandomVariate[DirichletDistribution[{1,1,1,1,1}]] and taking Sqrt we get a random vector $\{a, b, c, d, e \}$ that satisfies the condition.
Combining the steps in a function:
ClearAll[rvF]
rvF[dim_Integer, ss_: 1] := Sqrt[Append[#, 1 - Total[#]] & /@
RandomVariate[DirichletDistribution[ConstantArray[1, dim]], ss]];
Examples:
rvF[5]
{{0.866565, 0.266605, 0.134535, 0.307307, 0.255831}}
rvF[5,10]
To get a sample of size 10 for the 6-vector, $\{a, b, c, d, e , f \}$ with f distributed uniformly on [0,3]:
Join @@@ Transpose[{rvF[5, 10], List /@ RandomReal[3, {10}]}]
Further examples:
cp = ContourPlot3D[ x^2 + y^2 + z^2 == 1, {x, 0, 1}, {y, 0, 1}, {z, 0, 1},
Mesh -> None, ContourStyle -> Directive[Yellow, Opacity[0.7]]];
Show[cp, ListPointPlot3D[rvF[3, 1000]]]
$\begingroup$
$\endgroup$
3
To settle this: Append[Normalize[RandomReal[1, 5]], RandomReal[3]] produces a set of numbers that satisfies the OP's request. What distribution this tuple follows is a different kettle of fish.
-
$\begingroup$ If I am reading OP correctly @bob, he only wants the first 5 to satisfy the normalization condition, and the sixth to be within $(0,3)$.
Normalize[]will map the sixth number back to $(0,1)$. $\endgroup$ May 14, 2016 at 15:06 -
2$\begingroup$ If we want to sample uniformly from the set of all tuples that satisfy the constraint, we could do
Normalize /@ Select[RandomReal[1, {100, 5}], 0.1 < Norm[#] < 1 &]. The upper bound on the norm ensures (five-dimensional) spherical symmetry. The lower bound avoids numerical instability whenNormalizeing. This has to be repeated until a sufficient number of tuples are generated. Since this is rejection sampling, we don't know in advance how many will pass the test. $\endgroup$– SzabolcsMay 14, 2016 at 15:12 -
$\begingroup$ @Szabolcs, yes, one alternative would be to use something like a
While[]loop to return the first admissible tuple from the rejection sampling. $\endgroup$ May 14, 2016 at 15:19
Append[Normalize[RandomReal[1, 5]], RandomReal[3]]? $\endgroup$