[EM] High Order Numerical Cubature2.0

[Forest Simmons]

Good work!

That's much better than what I was thinking!

El dom., 6 de feb. de 2022 1:22 a. m., Daniel Carrera <dcarrera at gmail.com>
escribió:

> Ok. I think I have a general solution of the standard n-dimensional
> simplex.
>
> Let {x1, x2, ... , xn} be the coordinates.
>
> 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:
>
> Prob: f1(x1) ~ (1-x1)^(n-1)
>
> The cumulative distribution is:
>
> F1(x1) = Integral from 0 to x1 of f1(x1)
>
> If we integrate this and normalize the integral so that F1(1) = 1 we get
>
> F1(x1) = 1 - (1 - x1)^n
>
> 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.
>
> {s1, s2, ..., sn} <-- Sobol_nD()
>
> 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.
>
> Let F1(x1) = 1 - s1
>
> => x1 = 1 - s1^(1/n)
>
> Now that I have generated x1, I can generate other values recursively.
>
> 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}.
>
> The formula is slightly less ugly if I write it as a formula for x{k+1}
> instead of xk. Long story short:
>
> x1 = 1 - s1^(1/n)
>
> x{k+1} = 1 - [ 1 - s{k+1} * ( 1 - xk^(n-k) ) ]^(1 / (n-k) )
>
> As an example, let's compute the last coordinate xn. In this instance
> there is a lot of cancellation and we get:
>
> xn = sn * (1 - x{n-1})
>
> I haven't tested this, but this should produce uniformly distributed
> points {x1,x2,...,xn} inside the standard n-simplex.
>
> Cheers,
> Daniel
>
> On Sat, Feb 5, 2022 at 10:32 PM Daniel Carrera <dcarrera at gmail.com> wrote:
>
>> On Sat, Feb 5, 2022 at 3:55 PM Daniel Carrera <dcarrera at gmail.com> wrote:
>>
>>> https://postimg.cc/zV34pv0F
>>>
>>> Unfortunately, it doesn't seem to be space-filling. It looks like you
>>> have reinvented Sierpiński's gasket.
>>>
>>
>> Incidentally, my own idea of using multiple barycenters doesn't work
>> either. It gives points very clustered in the center of the simplex.
>>
>> 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:
>>
>> x + y < 1
>>
>> 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.
>>
>> Let `y = f(x) = 1 - x` be (not normalized) the probability of generating
>> the value `x`. The cumulative distribution is:
>>
>> F(x) = x - x^2/2
>>
>> 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:
>>
>> F(x) = u
>>
>> => x = 1 - sqrt(1 - u)
>>
>> 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:
>>
>> => y = (1 - x)*v
>>
>> The result is a nicely uniform distribution of points without the
>> clustering that you'd expect from random points:
>>
>> https://postimg.cc/XBbnD7fQ
>>
>> There's probably a way to generalize these ideas to higher dimensions.
>>
>> Cheers,
>> --
>> Dr. Daniel Carrera
>> Postdoctoral Research Associate
>> Iowa State University
>>
>
>
> --
> Dr. Daniel Carrera
> Postdoctoral Research Associate
> Iowa State University
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.electorama.com/pipermail/election-methods-electorama.com/attachments/20220206/e6be5706/attachment-0001.html>