Probability of random sphere lying inside the unit ball

I’ll do this for $n=2$ first, where the calculation is more tangible, and then build on that for the general case.

Parameterize the three points $(x_i,y_i)$ by the centre $(x,y)$ of the circle passing through them, the radius $r$ of that circle and three angles $\phi_i$, so that $(x_i,y_i)=(x,y)+r(\cos\phi_i,\sin\phi_i)$. The Jacobian is

\begin{eqnarray} \frac{\partial(x_1,y_1,x_2,y_2,x_3,y_3)}{\partial(x,y,r,\phi_1,\phi_2,\phi_3)} &=& \left|\matrix{ 1&0&1&0&1&0\\ 0&1&0&1&0&1\\ \cos\phi_1&\sin\phi_1&\cos\phi_2&\sin\phi_2&\cos\phi_3&\sin\phi_3\\ -r\sin\phi_1&r\cos\phi_1&0&0&0&0\\ 0&0&-r\sin\phi_2&r\cos\phi_2&0&0\\ 0&0&0&0&-r\sin\phi_3&r\cos\phi_3\\ }\right|\;. \end{eqnarray}

The last three rows, which reflect how the coordinates change with the angular variables, are mutually orthogonal. Their orthogonal complement is spanned by the three vectors that reflect how the three points change with the radial variable. We can transform the matrix accordingly, multiplying by an orthogonal matrix with determinant $1$:

\begin{eqnarray} \pmatrix{ 1&0&1&0&1&0\\ 0&1&0&1&0&1\\ \cos\phi_1&\sin\phi_1&\cos\phi_2&\sin\phi_2&\cos\phi_3&\sin\phi_3\\ -r\sin\phi_1&r\cos\phi_1&0&0&0&0\\ 0&0&-r\sin\phi_2&r\cos\phi_2&0&0\\ 0&0&0&0&-r\sin\phi_3&r\cos\phi_3\\ } \\[10pt] \times\pmatrix{ \cos\phi_1&0&0&-\sin\phi_1&0&0\\ \sin\phi&0&0&\cos\phi_1&0&0\\ 0&\cos\phi_2&0&0&-\sin\phi_2&0\\ 0&\sin\phi_2&0&0&\cos\phi_2&0\\ 0&0&\cos\phi_3&0&0&-\sin\phi_3\\ 0&0&\sin\phi_3&0&0&\cos\phi_3 } \\[10pt] = \pmatrix{ \cos\phi_1&\cos\phi_2&\cos\phi_3&-\sin\phi_1&-\sin\phi_2&-\sin\phi_3\\ \sin\phi_1&\sin\phi_2&\sin\phi_3&\cos\phi_1&\cos\phi_2&\cos\phi_3\\ 1&1&1&0&0&0\\ 0&0&0&r&0&0\\ 0&0&0&0&r&0\\ 0&0&0&0&0&r }\;. \end{eqnarray}

Thus

$$ \frac{\partial(x_1,y_1,x_2,y_2,x_3,y_3)}{\partial(x,y,r,\phi_1,\phi_2,\phi_3)}=r^3\left| \matrix{ \cos\phi_1&\cos\phi_2&\cos\phi_3\\ \sin\phi_1&\sin\phi_2&\sin\phi_3\\ 1&1&1 } \right|\;. $$

This determinant is twice the area of the triangle formed by the three radial unit vectors $(\cos\phi_i,\sin\phi_i)$. We want

\begin{eqnarray} P(\text{circle in unit disk}) &=& \frac1{\pi^3}\iiint_\limits{\sqrt{x^2+y^2}+r\lt1}\mathrm dx\mathrm dy\mathrm dr\int_0^{2\pi}\!\int_0^{2\pi}\!\int_0^{2\pi}\mathrm d\phi_1\mathrm d\phi_2\mathrm d\phi_3\frac{\partial(x_1,y_1,x_2,y_2,x_3,y_3)}{\partial(x,y,r,\phi_1,\phi_2,\phi_3)} \\ &=& \frac2{\pi^2}\iint_\limits{\rho+r\lt1}\rho\mathrm d\rho\mathrm dr\int_0^{2\pi}\!\int_0^{2\pi}\!\int_0^{2\pi}\mathrm d\phi_1\mathrm d\phi_2\mathrm d\phi_3\frac{\partial(x_1,y_1,x_2,y_2,x_3,y_3)}{\partial(x,y,r,\phi_1,\phi_2,\phi_3)} \\ &=& \frac1{\pi^2}\int_0^1(1-r)^2\mathrm dr\int_0^{2\pi}\!\int_0^{2\pi}\!\int_0^{2\pi}\mathrm d\phi_1\mathrm d\phi_2\mathrm d\phi_3\frac{\partial(x_1,y_1,x_2,y_2,x_3,y_3)}{\partial(x,y,r,\phi_1,\phi_2,\phi_3)} \\\\ &=& \frac1{\pi^2}\int_0^1(1-r)^2r^3\mathrm dr\int_0^{2\pi}\!\int_0^{2\pi}\!\int_0^{2\pi}\mathrm d\phi_1\mathrm d\phi_2\mathrm d\phi_3\left| \matrix{ \cos\phi_1&\cos\phi_2&\cos\phi_3\\ \sin\phi_1&\sin\phi_2&\sin\phi_3\\ 1&1&1 } \right|\\ \\ &=& \frac1{60\pi^2}\int_0^{2\pi}\!\int_0^{2\pi}\!\int_0^{2\pi}\mathrm d\phi_1\mathrm d\phi_2\mathrm d\phi_3\left| \matrix{ \cos\phi_1&\cos\phi_2&\cos\phi_3\\ \sin\phi_1&\sin\phi_2&\sin\phi_3\\ 1&1&1 } \right|\;. \end{eqnarray}

The integral is $(2\pi)^3$ times twice the average area of a triangle formed by three points independently uniformly distributed on the unit circle, which is $\frac3{2\pi}$ (see e.g. MathWorld). Thus

$$ P(\text{circle in unit disk})=\frac1{60\pi^2}\cdot(2\pi)^3\cdot2\cdot\frac3{2\pi}=\frac25\;, $$

in agreement with your result.

Now let’s generalize this for arbitrary $n$. We have $n+1$ points in $n$ dimensions, and thus $n(n+1)$ degrees of freedom. We parametrize them by the centre ($n$ coordinates), the radius ($1$ coordinate) and $n-1$ angular variables for each point ($(n+1)(n-1)=n^2-1$ coordinates). Differentiating with respect to the $n^2-1$ angular variables yields $n^2-1$ mutually orthogonal direction vectors with magnitude $r$ (corresponding to the last $2^2-1=3$ rows above). Transforming to separate the space spanned by these vectors and its orthogonal complement yields the Jacobian

$$ \left|\matrix{ \vec u_1&\cdots&\vec u_{n+1}\\ 1&\cdots&1 }\right|\;, $$

where $\vec u_k$ is the radial unit vector for the $k$-th point. This is $n!$ times the volume of the simplex formed by these $n+1$ unit vectors. Thus, if we denote the surface area of the $n$-sphere by $S_n$, the volume of the $n$-ball by $V_n$ and $n!$ times the average volume spanned by $n+1$ points independently uniformly distributed on the $(n-1)$-sphere by $T_n$, we have

\begin{eqnarray} P(\text{$n$-sphere in unit $n$-ball}) &=& \frac{S_{n-1}}{V_n^{n+1}}\iint_\limits{\rho+r\lt1}\mathrm d\rho\mathrm dr\rho^{n-1}r^{n^2-1}\int\mathrm d\Omega^{n+1}\left|\matrix{ \vec u_1&\cdots&\vec u_{n+1}\\ 1&\cdots&1 }\right| \\ &=& \frac{S_{n-1}}{nV_n^{n+1}}\int_0^1\mathrm dr(1-r)^nr^{n^2-1}\int\mathrm d\Omega^{n+1}\left|\matrix{ \vec u_1&\cdots&\vec u_{n+1}\\ 1&\cdots&1 }\right| \\ &=& \frac{S_{n-1}^{n+2}}{nV_n^{n+1}}B\left(n+1,n^2\right)T_n\;, \end{eqnarray}

where

$$ B(x,y)=\frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)} $$

is the beta function.

The volume and surface area are given by (see e.g. Wikipedia)

$$ V_n=\frac{\pi^{\frac n2}}{\Gamma\left(\frac n2+1\right)} $$

and

$$ S_{n-1}=nV_n\;, $$

respectively. From What is the expected volume of the simplex formed by $n+1$ points independently uniformly distributed on $\mathbb S^{n-1}$? we have

$$ T_n=\Xi\left(\frac{n^2}2\right)\Xi\left(\frac n2\right)^{-n}\prod_{l=1}^{n-1}\Xi\left(\frac l2\right) $$

with

$$ \Xi(n):=\frac{\Gamma\left(n+\frac12\right)}{\Gamma(n)}\;. $$

Thus

\begin{eqnarray} P(\text{$n$-sphere in unit $n$-ball}) &=& \frac{S_{n-1}^{n+2}}{nV_n^{n+1}}B\left(n+1,n^2\right)T_n \\ &=& n^{n+1}\pi^\frac n2\frac{\Gamma\left(n^2\right)\Gamma(n+1)}{\Gamma\left(\frac n2+1\right)\Gamma\left(n^2+n+1\right)}\Xi\left(\frac{n^2}2\right)\Xi\left(\frac n2\right)^{-n}\prod_{l=1}^{n-1}\Xi\left(\frac l2\right)\;. \end{eqnarray}

For $n=2$ we can recover the above result,

$$ P(\text{circle in unit disk})=2^3\pi\frac{\Gamma(4)\Gamma(3)}{\Gamma(2)\Gamma(7)}\frac{\Xi(2)\Xi\left(\frac12\right)}{\Xi(1)\Xi(1)}=\frac25\;, $$

and for $n=3$ we obtain

$$ P(\text{sphere in unit ball})=3^4\pi^\frac32\frac{\Gamma(9)\Gamma(4)}{\Gamma\left(\frac52\right)\Gamma(13)}\frac{\Xi\left(\frac92\right)\Xi\left(\frac12\right)\Xi(1)}{\Xi\left(\frac32\right)\Xi\left(\frac32\right)\Xi\left(\frac32\right)}=\frac{24\pi^2}{1925}\approx0.123\;. $$