# Sampling from intersection of $n$-sphere and hyperplane

Edit

How do I uniformly sample from intersection of $$n$$-sphere centered at 0 and an arbitrary hyperplane for $$n\ge3$$?

Below is an approach using FindInstance to sample $$x$$ with $$\langle x, a\rangle=b$$ and $$\|x\|=r$$ but it's slow and not uniform, any tips?

Clear[x];
d = 3;
r = 2;
a = Boole@Array[# == 1 &, d]; (*1,0,0*)
b = 0.5;
c = Normalize[ConstantArray[1, d]]; (*1,1,1*);
x = Array[x0, d];
eqPlane = x . a == b;
eqSphere = Total[x^2] == r^2;
vals = x /. FindInstance[{eqPlane, eqSphere}, x, Reals, 200];
Graphics3D@Point[vals]


Original

Suppose $$a,b,x$$ are unit vectors in $$\mathbb{R}^d$$ with $$a$$, $$b$$, $$d$$, $$\langle a, x\rangle$$ known. How do I plot the distribution of $$\langle b, x\rangle$$ assuming that all feasible $$x$$ are "equally likely"?

Below is a naive attempt using geometric regions for $$d=3$$. I'm looking for something faster that works for $$d \in [3,1000]$$.

d = 3;
a = Boole@Array[# == 1 &, d]; (* 1, 0, 0 *)
b = Normalize[ConstantArray[1, d]]; (* 1, 1, 1 *)
ax = 0.5;

eps = 0.1;
shell = SphericalShell[ConstantArray[0, d], {1 - eps, 1 + eps}];
plane = Hyperplane[a, ax];
intersect =
DiscretizeRegion[RegionIntersection[shell, plane],
AccuracyGoal -> 2, MaxCellMeasure -> {"Length" -> 1}];
xs = RandomPoint[intersect, 1000];
Histogram[# . b & /@ xs]


• Apologies if this isn't helpful, but the intersection of a plane with a spherical shell is always an annulus, right? Since b is always the same in all components {1/Sqrt[d], ...,1/Sqrt[d]}, dotting it with x just gives the (scaled by 1/Sqrt[d]) Total of x's components. So the probability x is a specific value $p(x = x_0)$ would just be the integral over the annular region $R$ for all coordinates that sum to $x_0$ right? So for the 3D case: $$p(x = x_0) \propto \iiint_\limits{R} \left(\mathrm{x_1+x_2 + x_3 == x_0}\right) ~\mathrm{d}x_1 \mathrm{d}x_2 \mathrm{d}x_3$$
– ydd
Commented Oct 15, 2023 at 14:03
• ...I am not sure how to set up this integral yet, but maybe it would be faster than random sampling and histogramming? And it seems like it could generalize to higher dimensions. Also I think you can always reduce dimension by 1, because the region where the hyper-plane and shell intersect will always lie together in that hyper-plane right?
– ydd
Commented Oct 15, 2023 at 14:03
• So this answer gives R code to visualize histogram for 3 dimensions and some theory, I was looking for Mathematica code to redo the diagram in higher d stats.stackexchange.com/a/628799/511 Commented Oct 15, 2023 at 15:28
• Is $<b, x>$ the dot product of $b$ and $x$? I ask because $<b, x>$ is not standard statistical notation (at least not to me). If so, after the update in your question, is the question still the distribution of the dot product of $b$ and $x$ given the value of the dot product of $a$ and $x$ where $a$ is an arbitrary unit vector?
– JimB
Commented Oct 17, 2023 at 20:23
• This whole question boils down to the distribution of one Cartesian component for random variates drawn from a $(d-2)$-sphere. Sketch of proof: use coordinate freedom to set $\vec{a} = a \hat{e}_2$ and $\vec{b} = b_1 \hat{e}_1 + b_2 \hat{e}_2$. The intersection of your hyperplace $\vec{a} \cdot \vec{x} = c$ and a $(d-1)$-sphere will be a $(d-2)$-sphere. And for points on this intersection, we have $\vec{b} \cdot \vec{x} = b_1 x_1 + b_2 x_2 = b_1 x_1 + b_2 c/a$. So if you find the distribution of $x_1$ on a $(d-2)$-sphere, the desired distribution is just a linear transform of that. Commented Oct 17, 2023 at 20:58

Let $$a$$ and $$b$$ be known unit vectors of dimension $$d$$. $$X$$ is a random unit vector on a $$d$$-dimensional sphere. We want to know the distribution of $$B=b.X$$ given that $$a.X=A$$.

Without loss of generality we only have to consider unit vectors of $$a$$ and $$b$$ in the following form:

$$a=(1,0,\ldots,0)$$

$$b=(\rho,\sqrt{1-\rho^2},0,\ldots,0)$$

where we have $$a.b=\rho$$. We write for the first two dimensions of the random unit vector:

$$(U,V)=(X_1,X_2)/\sqrt{X_1^2+X_2^2+\cdots +X_d^2}$$

with $$X_i \sim N(0,1)$$ (and all random variables $$X_i$$ independent of each other). This has the same distribution as

$$(X_1,X_2)/\sqrt{X_1^2+X_2^2+\chi^2}$$

where $$\chi^2 \sim \chi^2_{d-2}$$ and is independent of $$X_1$$ and $$X_2$$.

Then we find the joint distribution of $$(a.X,b.X)=(A,B)$$ where

$$(A,B)=(U,\rho U+\sqrt{1-\rho^2}V)$$

The marginal distribution of $$A=a.X$$ is implemented as a PERTDistribution in Mathematica. The joint distribution of $$A$$ and $$B$$ is divided by the marginal distribution of $$A$$ to obtain the conditional distribution of $$B$$ given $$A$$.

Below is the Mathematica code to obtain the joint, marginal, and conditional distributions.

(* Joint distribution of first two normalized values and the norm *)
distuvz = TransformedDistribution[{x1/Sqrt[x1^2 + x2^2 + z],
x2/Sqrt[x1^2 + x2^2 + z],
Sqrt[x1^2 + x2^2 + z]}, {x1 \[Distributed] NormalDistribution[],
x2 \[Distributed] NormalDistribution[],
z \[Distributed] ChiSquareDistribution[d - 2]}];
jointPDFuvz = FullSimplify[PDF[distuvz, {u, v, w}],
Assumptions -> {w > 0, u^2 + v^2 < 1}]
(* (2^(1-d/2) \[ExponentialE]^(-(w^2/2)) (-((-1+u^2+v^2) \
w^2))^(d/2))/(π (-1+u^2+v^2)^2 w Gamma[-1+d/2]) *)

(* Integrate out z *)
pdfuv = Integrate[jointPDFuvz, {w, 0, ∞}, Assumptions -> d > 0]
(* ((-2+d) (1-u^2-v^2)^(1/2 (-4+d)))/(2 π) *)

(* Construct joint distribution of first two normalized values *)
dist2 = ProbabilityDistribution[pdfuv, {u, -1, 1}, {v, -Sqrt[1 - u^2], Sqrt[1 - u^2]}];


Note that I couldn't get Mathematica to obtain pdfuv directly but had to include Sqrt[x1^2+x2^2+z] and then integrate it out.

(* Now construct distribution of A=a.X and B=b.X *)
dist3 = TransformedDistribution[{u, ρ u + Sqrt[1 - ρ^2] v}, {u, v} \[Distributed] dist2];
pdfAB = FullSimplify[PDF[dist3, {A, B}] // PiecewiseExpand, Assumptions -> {-1 < A < 1,
A ρ - Sqrt[(1 - A^2) (1 - ρ^2)] < B < A ρ + Sqrt[(1 - A^2) (1 - ρ^2)], d > 0}]

(* Example of joint distribution of a.X and b.X *)
parms = {d -> 7, ρ -> 0};
Plot3D[pdfAB /. parms, {A, -1, 1}, {B, -1, 1},
PlotRange -> {{-1, 1}, {-1, 1}, {0, Automatic}},
PlotLabel -> "Joint distribution of A and B given that d=7 and ρ=0"]


(* Marginal distribution of A=a.X *)
distAz = TransformedDistribution[{x1/Sqrt[x1^2 + z],
Sqrt[x1^2 + z]}, {x1 \[Distributed] NormalDistribution[0, 1],
z \[Distributed] ChiSquareDistribution[d - 1]}];
pdfA = Integrate[PDF[distAz, {A, z}], {z, 0, ∞}, Assumptions -> -1 < A < 1 && d > 0];
(* ((1-A^2)^(-(3/2)+d/2) Gamma[d/2])/(Sqrt[π] Gamma[1/2 (-1+d)]) *)

(* Example of marginal distribution of A *)
parms = {d -> 7};
Plot[pdfA /. parms, {A, -1, 1},
PlotLabel -> "Marginal distribution A given that d=7"]


(* Finally the conditional distribution of B=b.X given A=a.X *)
pdfBGivenA = FullSimplify[pdfAB/pdfA // PiecewiseExpand, Assumptions -> {-1 < A < 1,
A ρ - Sqrt[(1 - A^2) (1 - ρ^2)] < B < A ρ + Sqrt[(1 - A^2) (1 - ρ^2)]}]
pdfBGivenA = FullSimplify[pdfBGivenA, Assumptions -> -1 + A^2 + B^2 - 2 A B ρ + ρ^2 < 0]

(* ((1 - A^2 - B^2 + 2 A B ρ - ρ^2)^(1/2 (-4 + d)) ((-1 + A^2) (-1 + ρ^2))^((3 - d)/2) Gamma[1/2 (-1 + d)])/(Sqrt[π] Gamma[-1 + d/2]) *)


This is the same as

FullSimplify[PDF[PERTDistribution[{A ρ - Sqrt[(1 - A^2) (1 - ρ^2)],
A ρ + Sqrt[(1 - A^2) (1 - ρ^2)]}, A ρ, d - 4], B],
Assumptions -> -1 + A^2 + B^2 - 2 A B ρ + ρ^2 < 0]


An example of the conditional distribution:

parms = {d -> 15, ρ -> 3/4, A -> 0.5};
Plot[pdfBGivenA /. parms, {B,
A ρ - Sqrt[(1 - A^2) (1 - ρ^2)] /. parms,
A ρ + Sqrt[(1 - A^2) (1 - ρ^2)] /. parms},
PlotRange -> {{-1, 1}, {0, Automatic}},
PlotLabel -> "Conditional distribution of B given A=0.5, d=15, and ρ=3/4",
PlotRangeClipping -> False]


Here is how you randomly sample circle x^2+y^2==1 and sphere x^2+y^2+z^2==1. It is same for any n-sphere in n-variables x^2+y^2+...+q^2==1.

Normalize/@RandomVariate[NormalDistribution[], {1000, 2}];
Graphics[Point[%]]

Normalize/@RandomVariate[NormalDistribution[], {1000, 3}];
Graphics3D[Point[%]]


Normalize /@ RandomVariate[NormalDistribution[], {300, 2}];
Show[ContourPlot3D[
Evaluate[{eqPlane, eqSphere}], {x0[1], -2, 2}, {x0[2], -2,
2}, {x0[3], -2, 2}, ContourStyle -> Opacity[0.5], Mesh -> None],
Graphics3D[{Point[(Prepend[Sqrt[15]/2 #, 0.5] & /@ %)]}]]


Update:

To show that my method works for general case I chose random numbers for a and b.

a = Normalize[{-5, 3, -8}]
b = Normalize[{-9, 10, 6}]
pl = {x, y, z} . a == 0
sp = x^2 + y^2 + z^2 == 1

Normalize /@ RandomVariate[NormalDistribution[], {100, 2}];
Show[ContourPlot3D[
Evaluate[{pl, sp}], {x, -2, 2}, {y, -2, 2}, {z, -2, 2},
Mesh -> None, ContourStyle -> Opacity[0.5]],
Graphics3D[{Point[{-((8 #[[1]])/Sqrt[89]) - (15 #[[2]])/(
7 Sqrt[178]), -(1/7) Sqrt[89/2] #[[2]], (
35 Sqrt[178] #[[1]] - 24 Sqrt[89] #[[2]])/(
623 Sqrt[2])} & /@ %]}]]

Normalize /@ RandomVariate[NormalDistribution[], {100000, 2}];
Histogram[# .
b & /@ ({-((8 #[[1]])/Sqrt[89]) - (15 #[[2]])/(
7 Sqrt[178]), -(1/7) Sqrt[89/2] #[[2]], (
35 Sqrt[178] #[[1]] - 24 Sqrt[89] #[[2]])/(
623 Sqrt[2])} & /@ %)]


{-(5/(7 Sqrt[2])), 3/(7 Sqrt[2]), -((4 Sqrt[2])/7)}

{-(9/Sqrt[217]), 10/Sqrt[217], 6/Sqrt[217]}

-((5 x)/(7 Sqrt[2])) + (3 y)/(7 Sqrt[2]) - (4 Sqrt[2] z)/7 == 0

x^2 + y^2 + z^2 == 1


• Simpler: Normalize /@ RandomVariate[NormalDistribution[], {1000, 3}] Commented Oct 17, 2023 at 20:05
• This seems hardwired for a specific hyperplane, rather than intersection with a generic hyperplane Commented Oct 17, 2023 at 21:55
• @Yaroslav Bulatov: It is you who chose the example. It is exactly your example in your Edit section of your question where you had problems with uniform sampling. Commented Oct 17, 2023 at 23:20
• @Yaroslav Bulatov: "This seems hardwired for a specific hyperplane, rather than intersection with a generic hyperplane" - isn't it straightforward that if you can uniformly sample any hypersphere you can sample uniformly any intersection of hyperplane with hypersphere? Commented Oct 17, 2023 at 23:31

This answer is not in competition for the bounty in that I don't have time right now to provide the necessary theoretical justifications for the conditional distribution of the original question.

The marginal distributions of $$a.X$$ and $$b.X$$ are identical with that distribution being implemented in Mathematica as

PERTDistribution[{-1, 1}, 0, d-3]


For all possibilities where $$a.b=\rho$$, the distribution of $$b.X$$ given that $$a.X=u$$ (i.e., the original question), one can generate a large random sample from the bivariate distribution of $$a.X$$ and $$b.X$$ using

$$a=(1,0,\ldots,0)$$

$$b=(\rho,\sqrt{1-\rho^2},0,\ldots,0$$

and then examine a slice of that bivariate distribution where $$u-\delta< a.X for a small value of $$\delta$$.

n = 1000000;
d = 7;
ρ = 0.6;
SeedRandom[12345];
x1 = RandomVariate[NormalDistribution[], n];
x2 = RandomVariate[NormalDistribution[], n];
chi = RandomVariate[ChiSquareDistribution[d - 2], n];
(* Normalize *)
x = {#[[1]], #[[2]]}/Sqrt[#[[1]]^2 + #[[2]]^2 + #[[3]]] & /@ Transpose[{x1, x2, chi}];
ax = x[[All, 1]]; (* a = {1, 0, 0, ..., 0} *);
bx = ({ρ, Sqrt[1 - ρ^2]} . #) & /@ x;
data = Transpose[{ax, bx}];
Histogram3D[data, {0.02}, "PDF", SphericalRegion -> True, RotationAction -> "Clip",
ImageSize -> Large, PlotRange -> {{-1, 1}, {-1, 1}, {0, All}}]


Now take a slice of that distribution where $$a.X=0.35$$.

u = 0.35;
δ = 0.01
data2 = Select[data, u - δ < #[[1]] < u + δ &];
Show[Histogram[data2[[All, 2]], "FreedmanDiaconis", "PDF", PlotRange -> {{-1, 1}, {0, All}}],
Plot[{PDF[PERTDistribution[{u ρ - Sqrt[(1 - u^2) (1 - ρ^2)], u ρ + Sqrt[(1 - u^2) (1 - ρ^2)]}, u ρ, d - 4], v]},
{v, -1, 1}, PlotRange -> {{-1, 1}, {0, Automatic}}]]


Without the necessary justification the following provides a good "visual fit": The conditional distribution for $$b.X$$ given that $$a.X=u$$ and $$a.b=\rho$$ is approximated by

PERTDistribution[{u ρ - Sqrt[(1 - u^2) (1 - ρ^2)], u ρ + Sqrt[(1 - u^2) (1 - ρ^2)]}, u ρ, d - 4]


for $$d>3$$. For $$d=3$$

ArcSinDistribution[{u ρ - Sqrt[(1 - u^2) (1 - ρ^2)], u ρ + Sqrt[(1 - u^2) (1 - ρ^2)}]


would be appropriate.

For any value of $$d$$, the means of $$U=a.X$$ and $$V=b.X$$ are zero with both variances $$1/d$$ as long as $$a$$ and $$b$$ are known unit vectors. The correlation of $$a.X$$ and $$b.X$$ is $$a.b$$ where $$.$$ means the dot product. When $$d$$ is large, then $$a.X$$ and $$b.X$$ has an approximate bivariate normal distribution.

Therefore (again with large values of $$d$$), the distribution of $$V$$ given that $$U=u$$ is approximately normal with mean $$u(a.b)$$ and variance $$(1-(a.b)^2)/d^2$$. (This follows from $$U$$ and $$V$$ having approximately a bivariate normal distribution and is an answer to your original question about the conditional distribution.)

But how do I justify the first paragraph? I don't have a formal proof at this point but consider the following simulations (while stealing from @azerbajdzan 's answer):

n = 100000;
d = 100;
SeedRandom[12345];
x = Normalize[#] & /@ RandomVariate[NormalDistribution[], {n, d}];
a = Normalize@RandomVariate[NormalDistribution[], d];
b = Normalize@RandomVariate[NormalDistribution[], d];
ax = (a . #) & /@ x;
bx = (b . #) & /@ x;

TableForm[{{Mean[ax], Variance[ax], Skewness[ax], Kurtosis[ax]},
{Mean[bx], Variance[bx], Skewness[bx], Kurtosis[bx]},
{0, 1./d, 0, 3}},
TableHeadings -> {{"a.X", "b.X", "Normal[0,1/d]"}, {"Mean", "Variance", "Skewness", "Kurtosis"}}]


Correlation[ax, bx]
(* 0.113887 *)
a . b
(* 0.11328 *)


By setting to smaller values of $$d$$ everything matches up except the skewness and kurtosis coefficients.

Special case: $$a.b=0$$

@azerbajdzan 's answer deals with the "updated" question as to how to sample. Here (and above) is an attempt to get an exact answer to the "original" question.

The distribution of $$a.X$$ is found in @whuber 's answer https://stats.stackexchange.com/questions/169812/inference-on-p-left-left-sum-i-1nx-i-right-sum-i-1nx-i2/169818#169818. That result says that the distribution of $$a.X$$ is the same for any unit vector $$a$$. Using Mathematica that result can be (essentially) duplicated.

TableForm[
Table[{d, TransformedDistribution[w[1]/Sqrt[Sum[w[i]^2, {i, d}]],
Table[w[i] \[Distributed] NormalDistribution[0, 1], {i, d}]]}, {d, 3, 10}],


Text removed: I had stated that if $$a.b=0$$ (i.e., $$a.X$$ and $$b.X$$ would be uncorrelated), then the distributions of $$a.X$$ and $$b.X$$ would be independent. That is not so. This is a case where zero correlation does NOT imply independence. This should have been obvious to me as not all combinations of $$a.X$$ and $$b.X$$ are possible.

What is needed to find the general result when $$a.b \neq 0$$ is the bivariate distribution of $$a.X$$ and $$b.X$$. But as noted above, for large values of $$d$$, there is a good approximation for that distribution.

• It is not "(while stealing from @Roman 's comment)" - it is stealing from my answer. @Roman's comment was only that multiplying by 1/Sqrt[x^2 + y^2] can be simplified with Normalize. Commented Oct 18, 2023 at 7:42
• @azerbajdzan. Oops! Sorry about that. I’ll fix that statement when I get back to a computer.
– JimB
Commented Oct 18, 2023 at 8:04
• @azerbajdzan Done.
– JimB
Commented Oct 18, 2023 at 15:35
• @JimB thanks for the mathematica code, as for shape of distribution of the original question, it's supposed to be shifted beta stats.stackexchange.com/a/628799/511 Commented Oct 19, 2023 at 1:02
• And to be clear, the PERTDistribution is a reparametrized beta distribution which includes the shifted beta distribution.
– JimB
Commented Oct 19, 2023 at 3:16