Mixing
![Find the minimum of $MA + MB$ when $M$ is on the circle [closed]](https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhZCTB4-bb-s-32SAm6ytjzFOU06cH9o8q-a_wq_UnH6lcd8hvws23CmBWHM6g256LzTX9yK1EiCCzTC1NSgtxlNk_J9hdr9bA3tD0Nqntbl521j4bLAm_uXPANuWLBhcCKTKzns2gaUodx/s72-c/1.jpg)
Gaussian measure of the standard simplex
Clash Royale CLAN TAG#URR8PPP
up vote
2
down vote
favorite
What is
$$int_Delta_n e^-(x_1^2 + x_2^2 + cdots + x_n^2)/2 d x$$
where $Delta_n = (x_1, x_2, dots, x_n) in mathbbR^n mid x_i geq 0, sum_i=1^n x_i leq 1$
The motivation is that $(2pi)^-n/2e^-(x_1^2 + cdots + x_n^2)/2$ is the multivariate Gaussian PDF, so this number is the probability that a random multivariate Gaussian lies in the standard simplex $Delta_n$.
The physical interpretation looks a little nicer when you replicate $Delta_n$ in every orthant, e.g. for $n=2$ it's a rotated square, for $n=3$ it's the octahedron.
integration normal-distribution simplex gaussian-integral
asked Jul 27 at 1:23
Chris Jones
702414
add a comment |Â
up vote
2
down vote
favorite
What is
$$int_Delta_n e^-(x_1^2 + x_2^2 + cdots + x_n^2)/2 d x$$
where $Delta_n = (x_1, x_2, dots, x_n) in mathbbR^n mid x_i geq 0, sum_i=1^n x_i leq 1$
The motivation is that $(2pi)^-n/2e^-(x_1^2 + cdots + x_n^2)/2$ is the multivariate Gaussian PDF, so this number is the probability that a random multivariate Gaussian lies in the standard simplex $Delta_n$.
The physical interpretation looks a little nicer when you replicate $Delta_n$ in every orthant, e.g. for $n=2$ it's a rotated square, for $n=3$ it's the octahedron.
integration normal-distribution simplex gaussian-integral
asked Jul 27 at 1:23
Chris Jones
702414
1
$n=2$ is doable by rotating the square 45° and realizing the coordinates are independent Gaussians, which are in the square iff they are at most $sqrt2/2$ in magnitude. But there is no rotation of the octahedron which admits a similar simplification...
– Mike Earnest
Jul 27 at 1:47
add a comment |Â
up vote
2
down vote
favorite
up vote
2
down vote
favorite
What is
$$int_Delta_n e^-(x_1^2 + x_2^2 + cdots + x_n^2)/2 d x$$
where $Delta_n = (x_1, x_2, dots, x_n) in mathbbR^n mid x_i geq 0, sum_i=1^n x_i leq 1$
The motivation is that $(2pi)^-n/2e^-(x_1^2 + cdots + x_n^2)/2$ is the multivariate Gaussian PDF, so this number is the probability that a random multivariate Gaussian lies in the standard simplex $Delta_n$.
The physical interpretation looks a little nicer when you replicate $Delta_n$ in every orthant, e.g. for $n=2$ it's a rotated square, for $n=3$ it's the octahedron.
integration normal-distribution simplex gaussian-integral
asked Jul 27 at 1:23
Chris Jones
702414
What is
$$int_Delta_n e^-(x_1^2 + x_2^2 + cdots + x_n^2)/2 d x$$
where $Delta_n = (x_1, x_2, dots, x_n) in mathbbR^n mid x_i geq 0, sum_i=1^n x_i leq 1$
The motivation is that $(2pi)^-n/2e^-(x_1^2 + cdots + x_n^2)/2$ is the multivariate Gaussian PDF, so this number is the probability that a random multivariate Gaussian lies in the standard simplex $Delta_n$.
The physical interpretation looks a little nicer when you replicate $Delta_n$ in every orthant, e.g. for $n=2$ it's a rotated square, for $n=3$ it's the octahedron.
integration normal-distribution simplex gaussian-integral
asked Jul 27 at 1:23
Chris Jones
702414
asked Jul 27 at 1:23
Chris Jones
702414
asked Jul 27 at 1:23
Chris Jones
702414
asked Jul 27 at 1:23
Chris Jones
702414
702414
1
$n=2$ is doable by rotating the square 45° and realizing the coordinates are independent Gaussians, which are in the square iff they are at most $sqrt2/2$ in magnitude. But there is no rotation of the octahedron which admits a similar simplification...
– Mike Earnest
Jul 27 at 1:47
add a comment |Â
1
$n=2$ is doable by rotating the square 45° and realizing the coordinates are independent Gaussians, which are in the square iff they are at most $sqrt2/2$ in magnitude. But there is no rotation of the octahedron which admits a similar simplification...
– Mike Earnest
Jul 27 at 1:47
1
1
$n=2$ is doable by rotating the square 45° and realizing the coordinates are independent Gaussians, which are in the square iff they are at most $sqrt2/2$ in magnitude. But there is no rotation of the octahedron which admits a similar simplification...
– Mike Earnest
Jul 27 at 1:47
$n=2$ is doable by rotating the square 45° and realizing the coordinates are independent Gaussians, which are in the square iff they are at most $sqrt2/2$ in magnitude. But there is no rotation of the octahedron which admits a similar simplification...
– Mike Earnest
Jul 27 at 1:47
add a comment |Â
1 Answer
1
active
oldest
votes
up vote
1
down vote
accepted
I somewhat doubt that a closed form exists for general $n$ even admitting familiar special functions. Even so, there are some observations that might be useful.
First, I'll generalize your notation a bit. It will be useful to look at "standard" simplices of varying sizes, which I'll refer to as follows:
$$
Delta_n(r)=leftvec xinmathbb R^n: x_ige 0, sum_i x_ile rright
$$
Also, to keep prefactors to a minimum, I'll use a rescaling of the error function that I'll call $E$:
$$
E(x)=int_0^xe^-y^2/2dy=sqrtfracpi2texterf(x/sqrt2)
$$
We can write the family of integrals you're interested in as:
$$
I_n(r)=int_Delta_n(r)e^-(x_1^2+...+x_n^2)/2dx_1...dx_n
$$
This integral can be factorized somewhat when written in the typical simplical format.
$$
I_n(r)=int_0^re^-x_1^2/2dx_1int_0^r-x_1e^-x_2^2/2dx_2int_0^r-x_1-x_2e^-x_3^2/2dx_3...int_0^r-x_1-...-x_n-1e^-x_n^2/2dx_n
$$
Here we can observe a recurrence relation, Notice that if we remove the leftmost integral, we get the same expression with one less integral and $r-x_1$ in place of $r$.
$$
I_n+1(r)=int_0^re^-x^2/2I_n(r-x)dx=int_0^re^-(r-x)^2/2I_n(x)dx
$$
Now we can enumerate some cases. $n=0$ isn't necessarily meaningful, but this choice satisfies the recurrence relation.
$$
I_0(r)=1
$$
$n=0$ is a of course just a Gaussian integral.
$$
I_1(r)=E(r)
$$
$n=2$ can be solved using Mike Earnest's insightful observation that the 4 such simplical regions can be rotated into a square, yielding a separable integral. This can also be shown (tediously) by differentiating with respect or $r$ under the integral sign in the recurrence relation for $n=1$.
$$
I_2(r)=frac 14int_-r/sqrt2^r/sqrt2int_-r/sqrt2^r/sqrt2e^-(x^2+y^2)/2dy dx=E^2(r/sqrt2)
$$
We can use the same trick to derive a 2nd order recurrence relation with a single integral, but It's not clearly more useful than the first.
$$
I_n+2(r)=sqrt2int_0^r e^-(r-x)^2/4E((r-x)/2)I_n(x)dx
$$
If you want to compute these values numerically, it is best to transform the functional recurrence relation into a discrete one. We can do that by assuming that $I_n(r)$ is an analytic function of $r$, and writing it in terms of its Maclaurian expansion $sum_p a_p^(n)r^p$. As a shorthand, let $g_k$ be the Maclaurian coefficients of the Gaussian.
$$
sum_pa^(n+1)_pr^p=int_0^rsum_k=0^infty g_k(r-x)^ksum_l=0^infty a^(n)_lx^ldr, g_k=begincasesfrac(-1)^k/22^k/2(k/2)! & k texteven \ 0 & k textodd endcases
$$
$$
sum_pa^(n+1)_pr^p=int_0^rsum_k=0^infty sum_l=0^inftysum_m=0^kg_ka^(n)_l(-1)^mbinomkmr^mx^l+k-mdr
$$
$$
sum_pa^(n+1)_pr^p=sum_k=0^infty sum_l=0^inftysum_m=0^kg_ka^(n)_lfrac(-1)^mk+l-m+1binomkmr^k+l+1
$$
Now, we let $c(p,l)=sum_m=0^p-l-1g_p-l-1frac(-1)^mp-mbinomp-l-1m$
$$
sum_pa^(n+1)_pr^p=sum_k=0^infty sum_l=0^infty c(k+l+1,l)a^(n)_lr^k+l+1
$$
Reindexing to collect like powers:
$$
sum_pa^(n+1)_pr^p=sum_p=1^inftysum_l=0^p-1c(p, l)a^(n)_lr^p
$$
We obtain a somewhat serviceable recurrence relation. We can start with the initial condition $a^(0)_p=delta_p0$.
$$
a^(n+1)_p=sum_l=0^p-1c(p, l)a^(n)_l
$$
To evaluate $I_n(1)$, these coefficients can be generated reasonably quickly, and the resulting series converges quite quickly. Here are the results of a simple python script, which agrees with the exact results we already have:
$$
beginarray hline
n & I_n(1) \ hline
0 & 1 \ hline
1 & 0.855624392 \ hline
2 & 0.425560334 \ hline
3 & 0.1438817189 \ hline
4 & 0.0365313603 \ hline
5 & 7.40668651cdot 10^-3 \ hline
6 & 1.24877137cdot 10^-3 \ hline
7 & 1.80133459cdot 10^-4 \ hline
8 & 2.27017316cdot 10^-5 \ hline
9 & 2.54005636ecdot 10^-6 \ hline
10 & 2.55531558cdot 10^-7 \ hline
endarray
$$
If $n$ is very large, there is also a useful approximation. The integral $int_Delta_n(1)x_1^2$ is perfectly doable, which allows us to write the mean value of the radius squared $rho^2=x_1^2+...+x_n^2$ on the standard simplex.
$$
langlerho^2rangle=frac2n(n+1)(n+2)
$$
Using this, along with the fact that the distribution of $rho^2$ becomes increasingly sharp for large $n$, allows us to approximate the Gaussian as a constant.
$$
I_n(1)approxfrac1n!e^frac-n(n+1)(n+2), n>>1
$$
Numerically, it appears the relative error of this approximation tends toward zero as $n$ increases.
answered Aug 3 at 4:05
Kajelad
1,848619
Actually, your recurrence relation for I_n(r) looks like it might be solvable via Fourier analysis! You showed that I_n is the n-fold convolution of the function e^-x^2/2. Taking Fourier transforms, convolution becomes multiplication F[e^-x^2/2]^n. If I recall, the Gaussian is an eigenfunction of the Fourier transform to some constant, so we get C^n e^-nx^2/2. Then we have to inverse transform to recover I_n.
– Chris Jones
Aug 3 at 22:55
1
That was my first thought as well. Unfortunately, the bounds of the integral become somewhat problematic, since it isn't a "proper" convolution. If we let $I_0(r)=Theta(r)$, then it can be shown by induction that all subsequent $I_n(r)$ are zero for $r<0$. In that case, the integral becomes a convolution of $I_n$ with $e^-x^2/2Theta(x)$. It might be possible to attack this convolution analytically, but it gets a bit messy as far as I can tell.
– Kajelad
Aug 3 at 23:00
add a comment |Â
1 Answer
1
active
oldest
votes
1 Answer
1
active
oldest
votes
active
oldest
votes
active
oldest
votes
up vote
1
down vote
accepted
I somewhat doubt that a closed form exists for general $n$ even admitting familiar special functions. Even so, there are some observations that might be useful.
First, I'll generalize your notation a bit. It will be useful to look at "standard" simplices of varying sizes, which I'll refer to as follows:
$$
Delta_n(r)=leftvec xinmathbb R^n: x_ige 0, sum_i x_ile rright
$$
Also, to keep prefactors to a minimum, I'll use a rescaling of the error function that I'll call $E$:
$$
E(x)=int_0^xe^-y^2/2dy=sqrtfracpi2texterf(x/sqrt2)
$$
We can write the family of integrals you're interested in as:
$$
I_n(r)=int_Delta_n(r)e^-(x_1^2+...+x_n^2)/2dx_1...dx_n
$$
This integral can be factorized somewhat when written in the typical simplical format.
$$
I_n(r)=int_0^re^-x_1^2/2dx_1int_0^r-x_1e^-x_2^2/2dx_2int_0^r-x_1-x_2e^-x_3^2/2dx_3...int_0^r-x_1-...-x_n-1e^-x_n^2/2dx_n
$$
Here we can observe a recurrence relation, Notice that if we remove the leftmost integral, we get the same expression with one less integral and $r-x_1$ in place of $r$.
$$
I_n+1(r)=int_0^re^-x^2/2I_n(r-x)dx=int_0^re^-(r-x)^2/2I_n(x)dx
$$
Now we can enumerate some cases. $n=0$ isn't necessarily meaningful, but this choice satisfies the recurrence relation.
$$
I_0(r)=1
$$
$n=0$ is a of course just a Gaussian integral.
$$
I_1(r)=E(r)
$$
$n=2$ can be solved using Mike Earnest's insightful observation that the 4 such simplical regions can be rotated into a square, yielding a separable integral. This can also be shown (tediously) by differentiating with respect or $r$ under the integral sign in the recurrence relation for $n=1$.
$$
I_2(r)=frac 14int_-r/sqrt2^r/sqrt2int_-r/sqrt2^r/sqrt2e^-(x^2+y^2)/2dy dx=E^2(r/sqrt2)
$$
We can use the same trick to derive a 2nd order recurrence relation with a single integral, but It's not clearly more useful than the first.
$$
I_n+2(r)=sqrt2int_0^r e^-(r-x)^2/4E((r-x)/2)I_n(x)dx
$$
If you want to compute these values numerically, it is best to transform the functional recurrence relation into a discrete one. We can do that by assuming that $I_n(r)$ is an analytic function of $r$, and writing it in terms of its Maclaurian expansion $sum_p a_p^(n)r^p$. As a shorthand, let $g_k$ be the Maclaurian coefficients of the Gaussian.
$$
sum_pa^(n+1)_pr^p=int_0^rsum_k=0^infty g_k(r-x)^ksum_l=0^infty a^(n)_lx^ldr, g_k=begincasesfrac(-1)^k/22^k/2(k/2)! & k texteven \ 0 & k textodd endcases
$$
$$
sum_pa^(n+1)_pr^p=int_0^rsum_k=0^infty sum_l=0^inftysum_m=0^kg_ka^(n)_l(-1)^mbinomkmr^mx^l+k-mdr
$$
$$
sum_pa^(n+1)_pr^p=sum_k=0^infty sum_l=0^inftysum_m=0^kg_ka^(n)_lfrac(-1)^mk+l-m+1binomkmr^k+l+1
$$
Now, we let $c(p,l)=sum_m=0^p-l-1g_p-l-1frac(-1)^mp-mbinomp-l-1m$
$$
sum_pa^(n+1)_pr^p=sum_k=0^infty sum_l=0^infty c(k+l+1,l)a^(n)_lr^k+l+1
$$
Reindexing to collect like powers:
$$
sum_pa^(n+1)_pr^p=sum_p=1^inftysum_l=0^p-1c(p, l)a^(n)_lr^p
$$
We obtain a somewhat serviceable recurrence relation. We can start with the initial condition $a^(0)_p=delta_p0$.
$$
a^(n+1)_p=sum_l=0^p-1c(p, l)a^(n)_l
$$
To evaluate $I_n(1)$, these coefficients can be generated reasonably quickly, and the resulting series converges quite quickly. Here are the results of a simple python script, which agrees with the exact results we already have:
$$
beginarray hline
n & I_n(1) \ hline
0 & 1 \ hline
1 & 0.855624392 \ hline
2 & 0.425560334 \ hline
3 & 0.1438817189 \ hline
4 & 0.0365313603 \ hline
5 & 7.40668651cdot 10^-3 \ hline
6 & 1.24877137cdot 10^-3 \ hline
7 & 1.80133459cdot 10^-4 \ hline
8 & 2.27017316cdot 10^-5 \ hline
9 & 2.54005636ecdot 10^-6 \ hline
10 & 2.55531558cdot 10^-7 \ hline
endarray
$$
If $n$ is very large, there is also a useful approximation. The integral $int_Delta_n(1)x_1^2$ is perfectly doable, which allows us to write the mean value of the radius squared $rho^2=x_1^2+...+x_n^2$ on the standard simplex.
$$
langlerho^2rangle=frac2n(n+1)(n+2)
$$
Using this, along with the fact that the distribution of $rho^2$ becomes increasingly sharp for large $n$, allows us to approximate the Gaussian as a constant.
$$
I_n(1)approxfrac1n!e^frac-n(n+1)(n+2), n>>1
$$
Numerically, it appears the relative error of this approximation tends toward zero as $n$ increases.
answered Aug 3 at 4:05
Kajelad
1,848619
Actually, your recurrence relation for I_n(r) looks like it might be solvable via Fourier analysis! You showed that I_n is the n-fold convolution of the function e^-x^2/2. Taking Fourier transforms, convolution becomes multiplication F[e^-x^2/2]^n. If I recall, the Gaussian is an eigenfunction of the Fourier transform to some constant, so we get C^n e^-nx^2/2. Then we have to inverse transform to recover I_n.
– Chris Jones
Aug 3 at 22:55
1
That was my first thought as well. Unfortunately, the bounds of the integral become somewhat problematic, since it isn't a "proper" convolution. If we let $I_0(r)=Theta(r)$, then it can be shown by induction that all subsequent $I_n(r)$ are zero for $r<0$. In that case, the integral becomes a convolution of $I_n$ with $e^-x^2/2Theta(x)$. It might be possible to attack this convolution analytically, but it gets a bit messy as far as I can tell.
– Kajelad
Aug 3 at 23:00
add a comment |Â
up vote
1
down vote
accepted
I somewhat doubt that a closed form exists for general $n$ even admitting familiar special functions. Even so, there are some observations that might be useful.
First, I'll generalize your notation a bit. It will be useful to look at "standard" simplices of varying sizes, which I'll refer to as follows:
$$
Delta_n(r)=leftvec xinmathbb R^n: x_ige 0, sum_i x_ile rright
$$
Also, to keep prefactors to a minimum, I'll use a rescaling of the error function that I'll call $E$:
$$
E(x)=int_0^xe^-y^2/2dy=sqrtfracpi2texterf(x/sqrt2)
$$
We can write the family of integrals you're interested in as:
$$
I_n(r)=int_Delta_n(r)e^-(x_1^2+...+x_n^2)/2dx_1...dx_n
$$
This integral can be factorized somewhat when written in the typical simplical format.
$$
I_n(r)=int_0^re^-x_1^2/2dx_1int_0^r-x_1e^-x_2^2/2dx_2int_0^r-x_1-x_2e^-x_3^2/2dx_3...int_0^r-x_1-...-x_n-1e^-x_n^2/2dx_n
$$
Here we can observe a recurrence relation, Notice that if we remove the leftmost integral, we get the same expression with one less integral and $r-x_1$ in place of $r$.
$$
I_n+1(r)=int_0^re^-x^2/2I_n(r-x)dx=int_0^re^-(r-x)^2/2I_n(x)dx
$$
Now we can enumerate some cases. $n=0$ isn't necessarily meaningful, but this choice satisfies the recurrence relation.
$$
I_0(r)=1
$$
$n=0$ is a of course just a Gaussian integral.
$$
I_1(r)=E(r)
$$
$n=2$ can be solved using Mike Earnest's insightful observation that the 4 such simplical regions can be rotated into a square, yielding a separable integral. This can also be shown (tediously) by differentiating with respect or $r$ under the integral sign in the recurrence relation for $n=1$.
$$
I_2(r)=frac 14int_-r/sqrt2^r/sqrt2int_-r/sqrt2^r/sqrt2e^-(x^2+y^2)/2dy dx=E^2(r/sqrt2)
$$
We can use the same trick to derive a 2nd order recurrence relation with a single integral, but It's not clearly more useful than the first.
$$
I_n+2(r)=sqrt2int_0^r e^-(r-x)^2/4E((r-x)/2)I_n(x)dx
$$
If you want to compute these values numerically, it is best to transform the functional recurrence relation into a discrete one. We can do that by assuming that $I_n(r)$ is an analytic function of $r$, and writing it in terms of its Maclaurian expansion $sum_p a_p^(n)r^p$. As a shorthand, let $g_k$ be the Maclaurian coefficients of the Gaussian.
$$
sum_pa^(n+1)_pr^p=int_0^rsum_k=0^infty g_k(r-x)^ksum_l=0^infty a^(n)_lx^ldr, g_k=begincasesfrac(-1)^k/22^k/2(k/2)! & k texteven \ 0 & k textodd endcases
$$
$$
sum_pa^(n+1)_pr^p=int_0^rsum_k=0^infty sum_l=0^inftysum_m=0^kg_ka^(n)_l(-1)^mbinomkmr^mx^l+k-mdr
$$
$$
sum_pa^(n+1)_pr^p=sum_k=0^infty sum_l=0^inftysum_m=0^kg_ka^(n)_lfrac(-1)^mk+l-m+1binomkmr^k+l+1
$$
Now, we let $c(p,l)=sum_m=0^p-l-1g_p-l-1frac(-1)^mp-mbinomp-l-1m$
$$
sum_pa^(n+1)_pr^p=sum_k=0^infty sum_l=0^infty c(k+l+1,l)a^(n)_lr^k+l+1
$$
Reindexing to collect like powers:
$$
sum_pa^(n+1)_pr^p=sum_p=1^inftysum_l=0^p-1c(p, l)a^(n)_lr^p
$$
We obtain a somewhat serviceable recurrence relation. We can start with the initial condition $a^(0)_p=delta_p0$.
$$
a^(n+1)_p=sum_l=0^p-1c(p, l)a^(n)_l
$$
To evaluate $I_n(1)$, these coefficients can be generated reasonably quickly, and the resulting series converges quite quickly. Here are the results of a simple python script, which agrees with the exact results we already have:
$$
beginarray hline
n & I_n(1) \ hline
0 & 1 \ hline
1 & 0.855624392 \ hline
2 & 0.425560334 \ hline
3 & 0.1438817189 \ hline
4 & 0.0365313603 \ hline
5 & 7.40668651cdot 10^-3 \ hline
6 & 1.24877137cdot 10^-3 \ hline
7 & 1.80133459cdot 10^-4 \ hline
8 & 2.27017316cdot 10^-5 \ hline
9 & 2.54005636ecdot 10^-6 \ hline
10 & 2.55531558cdot 10^-7 \ hline
endarray
$$
If $n$ is very large, there is also a useful approximation. The integral $int_Delta_n(1)x_1^2$ is perfectly doable, which allows us to write the mean value of the radius squared $rho^2=x_1^2+...+x_n^2$ on the standard simplex.
$$
langlerho^2rangle=frac2n(n+1)(n+2)
$$
Using this, along with the fact that the distribution of $rho^2$ becomes increasingly sharp for large $n$, allows us to approximate the Gaussian as a constant.
$$
I_n(1)approxfrac1n!e^frac-n(n+1)(n+2), n>>1
$$
Numerically, it appears the relative error of this approximation tends toward zero as $n$ increases.
answered Aug 3 at 4:05
Kajelad
1,848619
Actually, your recurrence relation for I_n(r) looks like it might be solvable via Fourier analysis! You showed that I_n is the n-fold convolution of the function e^-x^2/2. Taking Fourier transforms, convolution becomes multiplication F[e^-x^2/2]^n. If I recall, the Gaussian is an eigenfunction of the Fourier transform to some constant, so we get C^n e^-nx^2/2. Then we have to inverse transform to recover I_n.
– Chris Jones
Aug 3 at 22:55
1
That was my first thought as well. Unfortunately, the bounds of the integral become somewhat problematic, since it isn't a "proper" convolution. If we let $I_0(r)=Theta(r)$, then it can be shown by induction that all subsequent $I_n(r)$ are zero for $r<0$. In that case, the integral becomes a convolution of $I_n$ with $e^-x^2/2Theta(x)$. It might be possible to attack this convolution analytically, but it gets a bit messy as far as I can tell.
– Kajelad
Aug 3 at 23:00
add a comment |Â
up vote
1
down vote
accepted
up vote
1
down vote
accepted
I somewhat doubt that a closed form exists for general $n$ even admitting familiar special functions. Even so, there are some observations that might be useful.
First, I'll generalize your notation a bit. It will be useful to look at "standard" simplices of varying sizes, which I'll refer to as follows:
$$
Delta_n(r)=leftvec xinmathbb R^n: x_ige 0, sum_i x_ile rright
$$
Also, to keep prefactors to a minimum, I'll use a rescaling of the error function that I'll call $E$:
$$
E(x)=int_0^xe^-y^2/2dy=sqrtfracpi2texterf(x/sqrt2)
$$
We can write the family of integrals you're interested in as:
$$
I_n(r)=int_Delta_n(r)e^-(x_1^2+...+x_n^2)/2dx_1...dx_n
$$
This integral can be factorized somewhat when written in the typical simplical format.
$$
I_n(r)=int_0^re^-x_1^2/2dx_1int_0^r-x_1e^-x_2^2/2dx_2int_0^r-x_1-x_2e^-x_3^2/2dx_3...int_0^r-x_1-...-x_n-1e^-x_n^2/2dx_n
$$
Here we can observe a recurrence relation, Notice that if we remove the leftmost integral, we get the same expression with one less integral and $r-x_1$ in place of $r$.
$$
I_n+1(r)=int_0^re^-x^2/2I_n(r-x)dx=int_0^re^-(r-x)^2/2I_n(x)dx
$$
Now we can enumerate some cases. $n=0$ isn't necessarily meaningful, but this choice satisfies the recurrence relation.
$$
I_0(r)=1
$$
$n=0$ is a of course just a Gaussian integral.
$$
I_1(r)=E(r)
$$
$n=2$ can be solved using Mike Earnest's insightful observation that the 4 such simplical regions can be rotated into a square, yielding a separable integral. This can also be shown (tediously) by differentiating with respect or $r$ under the integral sign in the recurrence relation for $n=1$.
$$
I_2(r)=frac 14int_-r/sqrt2^r/sqrt2int_-r/sqrt2^r/sqrt2e^-(x^2+y^2)/2dy dx=E^2(r/sqrt2)
$$
We can use the same trick to derive a 2nd order recurrence relation with a single integral, but It's not clearly more useful than the first.
$$
I_n+2(r)=sqrt2int_0^r e^-(r-x)^2/4E((r-x)/2)I_n(x)dx
$$
If you want to compute these values numerically, it is best to transform the functional recurrence relation into a discrete one. We can do that by assuming that $I_n(r)$ is an analytic function of $r$, and writing it in terms of its Maclaurian expansion $sum_p a_p^(n)r^p$. As a shorthand, let $g_k$ be the Maclaurian coefficients of the Gaussian.
$$
sum_pa^(n+1)_pr^p=int_0^rsum_k=0^infty g_k(r-x)^ksum_l=0^infty a^(n)_lx^ldr, g_k=begincasesfrac(-1)^k/22^k/2(k/2)! & k texteven \ 0 & k textodd endcases
$$
$$
sum_pa^(n+1)_pr^p=int_0^rsum_k=0^infty sum_l=0^inftysum_m=0^kg_ka^(n)_l(-1)^mbinomkmr^mx^l+k-mdr
$$
$$
sum_pa^(n+1)_pr^p=sum_k=0^infty sum_l=0^inftysum_m=0^kg_ka^(n)_lfrac(-1)^mk+l-m+1binomkmr^k+l+1
$$
Now, we let $c(p,l)=sum_m=0^p-l-1g_p-l-1frac(-1)^mp-mbinomp-l-1m$
$$
sum_pa^(n+1)_pr^p=sum_k=0^infty sum_l=0^infty c(k+l+1,l)a^(n)_lr^k+l+1
$$
Reindexing to collect like powers:
$$
sum_pa^(n+1)_pr^p=sum_p=1^inftysum_l=0^p-1c(p, l)a^(n)_lr^p
$$
We obtain a somewhat serviceable recurrence relation. We can start with the initial condition $a^(0)_p=delta_p0$.
$$
a^(n+1)_p=sum_l=0^p-1c(p, l)a^(n)_l
$$
To evaluate $I_n(1)$, these coefficients can be generated reasonably quickly, and the resulting series converges quite quickly. Here are the results of a simple python script, which agrees with the exact results we already have:
$$
beginarray hline
n & I_n(1) \ hline
0 & 1 \ hline
1 & 0.855624392 \ hline
2 & 0.425560334 \ hline
3 & 0.1438817189 \ hline
4 & 0.0365313603 \ hline
5 & 7.40668651cdot 10^-3 \ hline
6 & 1.24877137cdot 10^-3 \ hline
7 & 1.80133459cdot 10^-4 \ hline
8 & 2.27017316cdot 10^-5 \ hline
9 & 2.54005636ecdot 10^-6 \ hline
10 & 2.55531558cdot 10^-7 \ hline
endarray
$$
If $n$ is very large, there is also a useful approximation. The integral $int_Delta_n(1)x_1^2$ is perfectly doable, which allows us to write the mean value of the radius squared $rho^2=x_1^2+...+x_n^2$ on the standard simplex.
$$
langlerho^2rangle=frac2n(n+1)(n+2)
$$
Using this, along with the fact that the distribution of $rho^2$ becomes increasingly sharp for large $n$, allows us to approximate the Gaussian as a constant.
$$
I_n(1)approxfrac1n!e^frac-n(n+1)(n+2), n>>1
$$
Numerically, it appears the relative error of this approximation tends toward zero as $n$ increases.
answered Aug 3 at 4:05
Kajelad
1,848619
I somewhat doubt that a closed form exists for general $n$ even admitting familiar special functions. Even so, there are some observations that might be useful.
First, I'll generalize your notation a bit. It will be useful to look at "standard" simplices of varying sizes, which I'll refer to as follows:
$$
Delta_n(r)=leftvec xinmathbb R^n: x_ige 0, sum_i x_ile rright
$$
Also, to keep prefactors to a minimum, I'll use a rescaling of the error function that I'll call $E$:
$$
E(x)=int_0^xe^-y^2/2dy=sqrtfracpi2texterf(x/sqrt2)
$$
We can write the family of integrals you're interested in as:
$$
I_n(r)=int_Delta_n(r)e^-(x_1^2+...+x_n^2)/2dx_1...dx_n
$$
This integral can be factorized somewhat when written in the typical simplical format.
$$
I_n(r)=int_0^re^-x_1^2/2dx_1int_0^r-x_1e^-x_2^2/2dx_2int_0^r-x_1-x_2e^-x_3^2/2dx_3...int_0^r-x_1-...-x_n-1e^-x_n^2/2dx_n
$$
Here we can observe a recurrence relation, Notice that if we remove the leftmost integral, we get the same expression with one less integral and $r-x_1$ in place of $r$.
$$
I_n+1(r)=int_0^re^-x^2/2I_n(r-x)dx=int_0^re^-(r-x)^2/2I_n(x)dx
$$
Now we can enumerate some cases. $n=0$ isn't necessarily meaningful, but this choice satisfies the recurrence relation.
$$
I_0(r)=1
$$
$n=0$ is a of course just a Gaussian integral.
$$
I_1(r)=E(r)
$$
$n=2$ can be solved using Mike Earnest's insightful observation that the 4 such simplical regions can be rotated into a square, yielding a separable integral. This can also be shown (tediously) by differentiating with respect or $r$ under the integral sign in the recurrence relation for $n=1$.
$$
I_2(r)=frac 14int_-r/sqrt2^r/sqrt2int_-r/sqrt2^r/sqrt2e^-(x^2+y^2)/2dy dx=E^2(r/sqrt2)
$$
We can use the same trick to derive a 2nd order recurrence relation with a single integral, but It's not clearly more useful than the first.
$$
I_n+2(r)=sqrt2int_0^r e^-(r-x)^2/4E((r-x)/2)I_n(x)dx
$$
If you want to compute these values numerically, it is best to transform the functional recurrence relation into a discrete one. We can do that by assuming that $I_n(r)$ is an analytic function of $r$, and writing it in terms of its Maclaurian expansion $sum_p a_p^(n)r^p$. As a shorthand, let $g_k$ be the Maclaurian coefficients of the Gaussian.
$$
sum_pa^(n+1)_pr^p=int_0^rsum_k=0^infty g_k(r-x)^ksum_l=0^infty a^(n)_lx^ldr, g_k=begincasesfrac(-1)^k/22^k/2(k/2)! & k texteven \ 0 & k textodd endcases
$$
$$
sum_pa^(n+1)_pr^p=int_0^rsum_k=0^infty sum_l=0^inftysum_m=0^kg_ka^(n)_l(-1)^mbinomkmr^mx^l+k-mdr
$$
$$
sum_pa^(n+1)_pr^p=sum_k=0^infty sum_l=0^inftysum_m=0^kg_ka^(n)_lfrac(-1)^mk+l-m+1binomkmr^k+l+1
$$
Now, we let $c(p,l)=sum_m=0^p-l-1g_p-l-1frac(-1)^mp-mbinomp-l-1m$
$$
sum_pa^(n+1)_pr^p=sum_k=0^infty sum_l=0^infty c(k+l+1,l)a^(n)_lr^k+l+1
$$
Reindexing to collect like powers:
$$
sum_pa^(n+1)_pr^p=sum_p=1^inftysum_l=0^p-1c(p, l)a^(n)_lr^p
$$
We obtain a somewhat serviceable recurrence relation. We can start with the initial condition $a^(0)_p=delta_p0$.
$$
a^(n+1)_p=sum_l=0^p-1c(p, l)a^(n)_l
$$
To evaluate $I_n(1)$, these coefficients can be generated reasonably quickly, and the resulting series converges quite quickly. Here are the results of a simple python script, which agrees with the exact results we already have:
$$
beginarray hline
n & I_n(1) \ hline
0 & 1 \ hline
1 & 0.855624392 \ hline
2 & 0.425560334 \ hline
3 & 0.1438817189 \ hline
4 & 0.0365313603 \ hline
5 & 7.40668651cdot 10^-3 \ hline
6 & 1.24877137cdot 10^-3 \ hline
7 & 1.80133459cdot 10^-4 \ hline
8 & 2.27017316cdot 10^-5 \ hline
9 & 2.54005636ecdot 10^-6 \ hline
10 & 2.55531558cdot 10^-7 \ hline
endarray
$$
If $n$ is very large, there is also a useful approximation. The integral $int_Delta_n(1)x_1^2$ is perfectly doable, which allows us to write the mean value of the radius squared $rho^2=x_1^2+...+x_n^2$ on the standard simplex.
$$
langlerho^2rangle=frac2n(n+1)(n+2)
$$
Using this, along with the fact that the distribution of $rho^2$ becomes increasingly sharp for large $n$, allows us to approximate the Gaussian as a constant.
$$
I_n(1)approxfrac1n!e^frac-n(n+1)(n+2), n>>1
$$
Numerically, it appears the relative error of this approximation tends toward zero as $n$ increases.
answered Aug 3 at 4:05
Kajelad
1,848619
edited Aug 3 at 6:00
answered Aug 3 at 4:05
Kajelad
1,848619
answered Aug 3 at 4:05
Kajelad
1,848619
answered Aug 3 at 4:05
Kajelad
1,848619
1,848619
Actually, your recurrence relation for I_n(r) looks like it might be solvable via Fourier analysis! You showed that I_n is the n-fold convolution of the function e^-x^2/2. Taking Fourier transforms, convolution becomes multiplication F[e^-x^2/2]^n. If I recall, the Gaussian is an eigenfunction of the Fourier transform to some constant, so we get C^n e^-nx^2/2. Then we have to inverse transform to recover I_n.
– Chris Jones
Aug 3 at 22:55
1
That was my first thought as well. Unfortunately, the bounds of the integral become somewhat problematic, since it isn't a "proper" convolution. If we let $I_0(r)=Theta(r)$, then it can be shown by induction that all subsequent $I_n(r)$ are zero for $r<0$. In that case, the integral becomes a convolution of $I_n$ with $e^-x^2/2Theta(x)$. It might be possible to attack this convolution analytically, but it gets a bit messy as far as I can tell.
– Kajelad
Aug 3 at 23:00
add a comment |Â
Actually, your recurrence relation for I_n(r) looks like it might be solvable via Fourier analysis! You showed that I_n is the n-fold convolution of the function e^-x^2/2. Taking Fourier transforms, convolution becomes multiplication F[e^-x^2/2]^n. If I recall, the Gaussian is an eigenfunction of the Fourier transform to some constant, so we get C^n e^-nx^2/2. Then we have to inverse transform to recover I_n.
– Chris Jones
Aug 3 at 22:55
1
That was my first thought as well. Unfortunately, the bounds of the integral become somewhat problematic, since it isn't a "proper" convolution. If we let $I_0(r)=Theta(r)$, then it can be shown by induction that all subsequent $I_n(r)$ are zero for $r<0$. In that case, the integral becomes a convolution of $I_n$ with $e^-x^2/2Theta(x)$. It might be possible to attack this convolution analytically, but it gets a bit messy as far as I can tell.
– Kajelad
Aug 3 at 23:00
Actually, your recurrence relation for I_n(r) looks like it might be solvable via Fourier analysis! You showed that I_n is the n-fold convolution of the function e^-x^2/2. Taking Fourier transforms, convolution becomes multiplication F[e^-x^2/2]^n. If I recall, the Gaussian is an eigenfunction of the Fourier transform to some constant, so we get C^n e^-nx^2/2. Then we have to inverse transform to recover I_n.
– Chris Jones
Aug 3 at 22:55
Actually, your recurrence relation for I_n(r) looks like it might be solvable via Fourier analysis! You showed that I_n is the n-fold convolution of the function e^-x^2/2. Taking Fourier transforms, convolution becomes multiplication F[e^-x^2/2]^n. If I recall, the Gaussian is an eigenfunction of the Fourier transform to some constant, so we get C^n e^-nx^2/2. Then we have to inverse transform to recover I_n.
– Chris Jones
Aug 3 at 22:55
1
1
That was my first thought as well. Unfortunately, the bounds of the integral become somewhat problematic, since it isn't a "proper" convolution. If we let $I_0(r)=Theta(r)$, then it can be shown by induction that all subsequent $I_n(r)$ are zero for $r<0$. In that case, the integral becomes a convolution of $I_n$ with $e^-x^2/2Theta(x)$. It might be possible to attack this convolution analytically, but it gets a bit messy as far as I can tell.
– Kajelad
Aug 3 at 23:00
That was my first thought as well. Unfortunately, the bounds of the integral become somewhat problematic, since it isn't a "proper" convolution. If we let $I_0(r)=Theta(r)$, then it can be shown by induction that all subsequent $I_n(r)$ are zero for $r<0$. In that case, the integral becomes a convolution of $I_n$ with $e^-x^2/2Theta(x)$. It might be possible to attack this convolution analytically, but it gets a bit messy as far as I can tell.
– Kajelad
Aug 3 at 23:00
add a comment |Â
Sign up or log in
StackExchange.ready(function ()
StackExchange.helpers.onClickDraftSave('#login-link');
);
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
StackExchange.ready(
function ()
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fmath.stackexchange.com%2fquestions%2f2863970%2fgaussian-measure-of-the-standard-simplex%23new-answer', 'question_page');
);
Post as a guest
Sign up or log in
StackExchange.ready(function ()
StackExchange.helpers.onClickDraftSave('#login-link');
);
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Sign up or log in
StackExchange.ready(function ()
StackExchange.helpers.onClickDraftSave('#login-link');
);
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Sign up or log in
StackExchange.ready(function ()
StackExchange.helpers.onClickDraftSave('#login-link');
);
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
1
$n=2$ is doable by rotating the square 45° and realizing the coordinates are independent Gaussians, which are in the square iff they are at most $sqrt2/2$ in magnitude. But there is no rotation of the octahedron which admits a similar simplification...
– Mike Earnest
Jul 27 at 1:47