<div dir="ltr"><div dir="ltr"><div class="gmail_default" style="font-family:trebuchet ms,sans-serif;font-size:small">Ok. I think I have a general solution of the standard n-dimensional simplex.</div><div class="gmail_default" style="font-family:trebuchet ms,sans-serif;font-size:small"><br></div><div class="gmail_default" style="font-family:trebuchet ms,sans-serif;font-size:small">Let {x1, x2, ... , xn} be the coordinates.</div><div class="gmail_default" style="font-family:trebuchet ms,sans-serif;font-size:small"><br></div><div class="gmail_default" style="font-family:trebuchet ms,sans-serif;font-size:small">For any given value of x1 the remaining coordinates can go from 0 to (1-x1) and span an (n-1)-simplex with volume V ~ (1-x1)^(n-1). To distribute the points uniformly, the probability of choosing a value x1 should be proportional to that volume:</div><div class="gmail_default" style="font-family:trebuchet ms,sans-serif;font-size:small"><br></div><div class="gmail_default" style="font-family:trebuchet ms,sans-serif;font-size:small">Prob: f1(x1) ~ (1-x1)^(n-1)</div><div class="gmail_default" style="font-family:trebuchet ms,sans-serif;font-size:small"><br></div><div class="gmail_default" style="font-family:trebuchet ms,sans-serif;font-size:small">The cumulative distribution is:</div><div class="gmail_default" style="font-family:trebuchet ms,sans-serif;font-size:small"><br></div><div class="gmail_default" style="font-family:trebuchet ms,sans-serif;font-size:small">F1(x1) = Integral from 0 to x1 of f1(x1)</div><div class="gmail_default" style="font-family:trebuchet ms,sans-serif;font-size:small"><br></div><div class="gmail_default" style="font-family:trebuchet ms,sans-serif;font-size:small">If we integrate this and normalize the integral so that F1(1) = 1 we get</div><div class="gmail_default" style="font-family:trebuchet ms,sans-serif;font-size:small"><br></div><div class="gmail_default" style="font-family:trebuchet ms,sans-serif;font-size:small">F1(x1) = 1 - (1 - x1)^n</div><div class="gmail_default" style="font-family:trebuchet ms,sans-serif;font-size:small"><br></div><div class="gmail_default" style="font-family:trebuchet ms,sans-serif;font-size:small">Let us now draw n values from the Sobol sequence. The way the Sobol sequence works, you specify the number of dimensions and Sobol_nD returns a single n-point from the unit N-cube.</div><div class="gmail_default" style="font-family:trebuchet ms,sans-serif;font-size:small"><br></div><div class="gmail_default" style="font-family:trebuchet ms,sans-serif;font-size:small">{s1, s2, ..., sn} <-- Sobol_nD()</div><div class="gmail_default" style="font-family:trebuchet ms,sans-serif;font-size:small"><br></div><div class="gmail_default" style="font-family:trebuchet ms,sans-serif;font-size:small">Since the cumulative distribution F1() goes from 0 to 1, I can set F1 = s1 or F1 = 1-s1. As it turns out, the latter gives a slightly nicer formula.</div><div class="gmail_default" style="font-family:trebuchet ms,sans-serif;font-size:small"><br></div><div class="gmail_default" style="font-family:trebuchet ms,sans-serif;font-size:small">Let F1(x1) = 1 - s1</div><div class="gmail_default" style="font-family:trebuchet ms,sans-serif;font-size:small"><br></div><div class="gmail_default" style="font-family:trebuchet ms,sans-serif;font-size:small">=> x1 = 1 - s1^(1/n)</div><div class="gmail_default" style="font-family:trebuchet ms,sans-serif;font-size:small"><br></div><div class="gmail_default" style="font-family:trebuchet ms,sans-serif;font-size:small">Now that I have generated x1, I can generate other values recursively.</div><div class="gmail_default" style="font-family:trebuchet ms,sans-serif;font-size:small"><br></div><div class="gmail_default" style="font-family:trebuchet ms,sans-serif;font-size:small">Suppose that I have generated values {x1, x2, ..., x{k-1}} and now I want to generate xk. We have that 0 < xk < 1 - x{k-1}. For any choice of xk the remaining coordinates span an (N-k)-simplex with volume (1-xk)^(N-k) so the probability of choosing xk has to be proportional to that. So we setup the probability distribution, integrate to get the cumulative distribution Fk(xk), and normalize so that Fk(xk = 1 - x{k-1}) = 1. When I got to this point the formula was slightly less ugly if I chose Fk(xk) = sk. Solve for xk and I obtain a formula for xk based on the kth Sobol number and the previous point x{k-1}.</div><div class="gmail_default" style="font-family:trebuchet ms,sans-serif;font-size:small"><br></div><div class="gmail_default" style="font-family:trebuchet ms,sans-serif;font-size:small">The formula is slightly less ugly if I write it as a formula for x{k+1} instead of xk. Long story short:</div><div class="gmail_default" style="font-family:trebuchet ms,sans-serif;font-size:small"><br></div><div class="gmail_default" style="font-family:trebuchet ms,sans-serif;font-size:small">x1 = 1 - s1^(1/n)<br></div><div class="gmail_default" style="font-family:trebuchet ms,sans-serif;font-size:small"><br></div><div class="gmail_default" style="font-family:trebuchet ms,sans-serif;font-size:small">x{k+1} = 1 - [ 1 - s{k+1} * ( 1 - xk^(n-k) ) ]^(1 / (n-k) )</div><div class="gmail_default" style="font-family:trebuchet ms,sans-serif;font-size:small"><br></div><div class="gmail_default" style="font-family:trebuchet ms,sans-serif;font-size:small">As an example, let's compute the last coordinate xn. In this instance there is a lot of cancellation and we get:</div><div class="gmail_default" style="font-family:trebuchet ms,sans-serif;font-size:small"><br></div><div class="gmail_default" style="font-family:trebuchet ms,sans-serif;font-size:small">xn = sn * (1 - x{n-1})</div></div><div><br></div><div><div class="gmail_default" style="font-family:"trebuchet ms",sans-serif;font-size:small">I haven't tested this, but this should produce uniformly distributed points {x1,x2,...,xn} inside the standard n-simplex.</div></div><div class="gmail_default" style="font-family:"trebuchet ms",sans-serif;font-size:small"><br></div><div class="gmail_default" style="font-family:"trebuchet ms",sans-serif;font-size:small">Cheers,</div><div class="gmail_default" style="font-family:"trebuchet ms",sans-serif;font-size:small">Daniel</div><div class="gmail_default" style="font-family:"trebuchet ms",sans-serif;font-size:small"><br></div><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Sat, Feb 5, 2022 at 10:32 PM Daniel Carrera <<a href="mailto:dcarrera@gmail.com">dcarrera@gmail.com</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div dir="ltr"><div style="font-family:"trebuchet ms",sans-serif;font-size:small"><span style="font-family:Arial,Helvetica,sans-serif">On Sat, Feb 5, 2022 at 3:55 PM Daniel Carrera <<a href="mailto:dcarrera@gmail.com" target="_blank">dcarrera@gmail.com</a>> wrote:</span><br></div></div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div dir="ltr"><div style="font-family:"trebuchet ms",sans-serif;font-size:small"><a href="https://postimg.cc/zV34pv0F" target="_blank">https://postimg.cc/zV34pv0F</a><br></div></div><div class="gmail_quote"><div style="font-family:"trebuchet ms",sans-serif;font-size:small"><br></div><div style="font-family:"trebuchet ms",sans-serif;font-size:small">Unfortunately, it doesn't seem to be space-filling. It looks like you have reinvented Sierpiński's gasket.</div></div></div>
</blockquote></div><br clear="all"><div><div style="font-family:"trebuchet ms",sans-serif;font-size:small">Incidentally, my own idea of using multiple barycenters doesn't work either. It gives points very clustered in the center of the simplex.</div></div><div style="font-family:"trebuchet ms",sans-serif;font-size:small"><br></div><div style="font-family:"trebuchet ms",sans-serif;font-size:small">I did find a solution that works great for a 2D simplex, but I don't see an obvious way to generalize it. We begin by generating a pair of points (u,v) from the 2D Sobol sequence. We want to map these to (x,y) inside the simplex in a way that retains at least some of the uniformity guarantees of the Sobol sequence. The simplex is defined by the boundary:</div><div style="font-family:"trebuchet ms",sans-serif;font-size:small"><br></div><div style="font-family:"trebuchet ms",sans-serif;font-size:small">x + y < 1</div><div style="font-family:"trebuchet ms",sans-serif;font-size:small"><br></div><div style="font-family:"trebuchet ms",sans-serif;font-size:small">That means for x = 0 there are more possible values of y than for x = 0.9, so we should produce more values close to x=0. This is a familiar problem in the context of drawing random numbers and it has a known solution.</div><div style="font-family:"trebuchet ms",sans-serif;font-size:small"><br></div><div style="font-family:"trebuchet ms",sans-serif;font-size:small">Let `y = f(x) = 1 - x` be (not normalized) the probability of generating the value `x`. The cumulative distribution is:</div><div style="font-family:"trebuchet ms",sans-serif;font-size:small"><br></div><div style="font-family:"trebuchet ms",sans-serif;font-size:small">F(x) = x - x^2/2</div><div style="font-family:"trebuchet ms",sans-serif;font-size:small"><br></div><div style="font-family:"trebuchet ms",sans-serif;font-size:small">Normalize it to `F(x) = 2x - x^2` so that we have the property that `F(0) = 0` and `F(1) = 1`. So we can let F(x) be equal to the first Sobol number:</div><div style="font-family:"trebuchet ms",sans-serif;font-size:small"><br></div><div style="font-family:"trebuchet ms",sans-serif;font-size:small">F(x) = u</div><div style="font-family:"trebuchet ms",sans-serif;font-size:small"><br></div><div style="font-family:"trebuchet ms",sans-serif;font-size:small">=> x = 1 - sqrt(1 - u)</div><div style="font-family:"trebuchet ms",sans-serif;font-size:small"><br></div><div style="font-family:"trebuchet ms",sans-serif;font-size:small">This gives an `x` with the correct probability distribution and we are free to draw `y` randomly between 0 and (1-x). We accomplish that with the second Sobol number:</div><div style="font-family:"trebuchet ms",sans-serif;font-size:small"><br></div><div style="font-family:"trebuchet ms",sans-serif;font-size:small">=> y = (1 - x)*v</div><div style="font-family:"trebuchet ms",sans-serif;font-size:small"><br></div><div style="font-family:"trebuchet ms",sans-serif;font-size:small">The result is a nicely uniform distribution of points without the clustering that you'd expect from random points:</div><div style="font-family:"trebuchet ms",sans-serif;font-size:small"><br></div><div style="font-family:"trebuchet ms",sans-serif;font-size:small"><a href="https://postimg.cc/XBbnD7fQ" target="_blank">https://postimg.cc/XBbnD7fQ</a></div><div style="font-family:"trebuchet ms",sans-serif;font-size:small"><br></div><div style="font-family:"trebuchet ms",sans-serif;font-size:small">There's probably a way to generalize these ideas to higher dimensions.</div><div style="font-family:"trebuchet ms",sans-serif;font-size:small"><br></div><div style="font-family:"trebuchet ms",sans-serif;font-size:small">Cheers,</div>-- <br><div dir="ltr"><div dir="ltr"><div><div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr"><font face="trebuchet ms, sans-serif">Dr. Daniel Carrera</font></div><div dir="ltr"><font face="trebuchet ms, sans-serif">Postdoctoral Research Associate</font></div><div><font face="trebuchet ms, sans-serif">Iowa State University</font></div></div></div></div></div></div></div></div></div></div></div>
</blockquote></div><br clear="all"><div><br></div>-- <br><div dir="ltr" class="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr"><font face="trebuchet ms, sans-serif">Dr. Daniel Carrera</font></div><div dir="ltr"><font face="trebuchet ms, sans-serif">Postdoctoral Research Associate</font></div><div><font face="trebuchet ms, sans-serif">Iowa State University</font></div></div></div></div></div></div></div></div></div></div></div>