OK, experimenting with how to make blog entries, and also providing some information on how to generate points that are all Pareto optimal but are otherwise as random as possible. I created these methods back in 2004 while working on my PhD thesis. Since then, they’ve been recreated many times over, and one can google on them and get lots of hits. Here’s a useful page.
http://stackoverflow.com/questions/3010837/sample-uniformly-at-random-from-an-n-dimensional-unit-simplex
The algorithm to generate points on a simplex
Here’s a picture of 1000 points in 3 dimensions.
Here’s the code. If you cut and paste this into Mathematica, you may need to delete and retype the “d-1” term, as some sort of hidden formatting seems to carry over and mess up the evaluation when I try it.
genSimplex[n_, d_] := Table[Differences[Sort[Flatten[{0, RandomReal[1, d – 1], 1}]]], {n}];
The algorithm generates points that are randomly distributed on an outer face of a simplex. The way to generate them is, for a d-dimensional problem…
1. Generate d-1 uniformly distributed random values in the range [0,1]
2. Add a 0 and a 1 to the list
3. Sort the list
4. Extract a list of the differences between the elements
You now have a list of random values that sum to 1 (so they are on a plane that is defined by points that sum to one) and that are otherwise independent of each other, so their dispersion is uniform.
Generating Points on Hypersphere
This algorithm uses the fact that multivariate normal points distributed with \[Mu] = {0,…,0} and \[CapitalSigma] = Identity Matrix are uniformly sprayed out in all directions. So generate the points, take the absolute value of the coordinate values to move them all to the positive quadrant, and then normalize them to project them onto the unit hypersphere.
genSphere[n_, d_] := Table[Normalize[Abs[RandomVariate[MultinormalDistribution[ConstantArray[0, d], IdentityMatrix[d]]]]], {n}]
Picture of 1000 pts in 3D, with the aspect ratio a little squashed.
Generating Points on a Plane And Then Rotating Them
This algorithmic approach first generates points that form a d-1 dimensional cloud embedded in d dimensions, so each point has value [0 x1 x2 …. xd]. Thus the points are on a plane orthogonal to a d-dimensional vector [1 0 …. 0]. We next create a rotation matrix that will rotate the vector to [1 1 1 …. 1], and orthogonalize it. We then just multiply each of the data points by this matrix.
First example is generating a could that is distributed uniformly over a d-1 hypercube.
genRotatedPlane[n_, d_] := Map[Transpose@Orthogonalize@ReplacePart[IdentityMatrix[d], {1, i_} -> 1].# &, Transpose@Prepend[RandomReal[1, {d – 1, n}], Array[0 &, n]]];
You may need to retype sections in case some odd hidden formatting carries through. Example data set, 1000 pts in 3D.
Similar to above, but rather than generating a d-1 hypercube, we generate a d-1 multi-variate normal cloud with μ = [0 0 … 0] and Σ as the identity matrix. Lots of other potential approaches to generating the initial data prior to rotation.
genRotatedPlaneNormal[n_, d_] :=
Map[Transpose@Orthogonalize@ReplacePart[IdentityMatrix[d], {1, i_} -> 1].# &, Transpose@Prepend[RandomVariate[NormalDistribution[], {d – 1, n}], Array[0 &, n]]];
Here’s a modification that removes a chunk of the data.
genRotatedPlaneNormalChunkMissing[n_, d_] :=
Map[Transpose@Orthogonalize@ReplacePart[IdentityMatrix[d], {1, i_} -> 1].# &,
Select[Transpose@Prepend[RandomVariate[NormalDistribution[], {d – 1, n}],
Array[0 &, n]], ! AllTrue[#, NonPositive] &]];