Surface Area of Unit Sphere

Paradoxically, surface area first grows, and then starts to decrease as the number of dimensions goes over 7

In[8]:=

unitarea := (2Pi^(d/2))/Gamma[d/2]

In[105]:=

Plot[unitarea, {d, 0, 15}, AxesLabel→ {"dimension", "area"}] ;

[Graphics:HTMLFiles/index_3.gif]

Point of maximum surface area

In[24]:=

FindRoot[D[area, d] == 0, {d, 5}]

Out[24]=

{d→7.25695}

Volume of Unit Sphere

In[1]:=

<<Graphics`ContourPlot3D`

In[5]:=

ContourPlot3D[Cos[Sqrt[x^2 + y^2 + z^2]], {x, -2, 2}, {y, -2, 2}, {z, -2, 2}, PlotPoints→ {5, 5}]

[Graphics:HTMLFiles/index_8.gif]

Out[5]=

-Graphics3D -

Similar thing happens with volume of unit sphere

In[9]:=

volume = Integrate[unitarea * r^(d - 1), {r, 0, 1}, Assumptions→ {Re[d] >0}]

Out[9]=

(2 π^(d/2))/(d Gamma[d/2])

In[20]:=

Plot[volume, {d, 0, 15}, AxesLabel→ {"dimension", "volume"}] ;

[Graphics:HTMLFiles/index_13.gif]

In[15]:=

volume/.d→5//N

Out[15]=

5.26379

Fraction in corners

As dimension grows, increasingly large fraction of mass is concentrated in the corners of unit cube (volume outside of the volume contained inside of an inscribed sphere)

corners = 1 - Integrate[unitarea * r^(d - 1), {r, 0, 1/2}, Assumptions→ {Re[d] >0}]

Out[39]=

1 - (2^(-d) π^(d/2))/Gamma[1 + d/2]

In[43]:=

Plot[corners, {d, 1, 10}] ;

[Graphics:HTMLFiles/index_19.gif]

In[45]:=

corners/.d→7//N

Out[45]=

0.963088

Visualization of multidimensional cube

7 dimensional cube,most 96% of the mass is concentrated in"corners". For a d dimensional cube, square of the diagonal is d, so as the dimension increases, distance from the center of the cube, to the corners goes to infinity (whereas distance to any facet stays constant), so corners become long spikes

[Graphics:HTMLFiles/index_22.gif]

E[{|X|}^2] for standard Normal

For standard normal(identity covariance matrix, 0 mean), expected norm squared is just d (number of dimenions)

In[22]:=

2 Plot[d, {d, 1, 10}, AxesLabel→ {"dimensions", E[{|X |} }] ;

[Graphics:HTMLFiles/index_25.gif]

This can be derived by expanding the sum inside the norm^2 so that integrals decouple and using integration by parts to individual integrals

E[{|X|}] for standard Normal

In[25]:=

expdev = 2^(1/2) Gamma[d/2 + 1/2]/Gamma[d/2] ;

In[30]:=

Plot[expdev, {d, 1, 100}, AxesLabel→ {"dimensions", "E[{|X|}]"}] ;

[Graphics:HTMLFiles/index_29.gif]

This can be derived by converting to hyperspherical coordinates, and then using integration by parts/setting up a recurrence to integrate sin^ix and x^iexp(-x^2)

Using the duplication formula for Gamma function we can rewrite

In[78]:=

expdev/.{Gamma[x_ + 1/2] → (Pi^(1/2) d !)/(2^d (d/2) !), Gamma[x_] → (x - 1) !}

Out[78]=

(2^(1/2 - d) π^(1/2) d !)/((-1 + d/2) ! d/2 !)

Another simplification

In[79]:=

%/.(-1 + x_) ! →x !/x

Out[79]=

(2^(-1/2 - d) d π^(1/2) d !)/(d/2 !)^2

Now we can use Stirling's approximation

In[80]:=

%/.x_ ! →x^xExp[-x] (2π x)^(1/2)

Out[80]=

d^(1/2)

Actually (d - 1/2)^(1/2) makes an even better approximation, as can be seen from this graph (dashed is the approximation)

In[101]:=

Plot[{(d - 1/2)^(1/2), expdev}, {d, 1, 10}, PlotStyle→ {{Dashing[{0.1}], Thickness[0.01]}, {Thickness[0.005]}}, AxesLabel→ {"dimensions", "E[{|X|}]"}] ;

[Graphics:HTMLFiles/index_41.gif]

Distribution of mass as function of radius for standard Normal

In[107]:=

density := 1/(2Pi)^(d/2) Exp[-1/2r^2] ;

In[114]:=

area = unitarea * r^(d - 1) ;

In[115]:=

density * area

Out[115]=

(2^(1 - d/2) ^(-r^2/2) r^(-1 + d))/Gamma[d/2]

For one dimension, just right half of the standard normal

In[145]:=

Plot[density * area/.{d→1}, {r, 0, 6}] ;

[Graphics:HTMLFiles/index_47.gif]

For 4 and 16 dimensions, the mass starts to drift right

In[241]:=

Plot[density * area/.{d→4}, {r, 0, 6}] ;

[Graphics:HTMLFiles/index_49.gif]

In[144]:=

Plot[density * area/.{d→16}, {r, 0, 6}, PlotRange→ {0, 0.6}] ;

[Graphics:HTMLFiles/index_51.gif]

As was found in previous section, the mean of this distribution is about (d - 1/2)^(1/2)

We can use our expression for exact mean of the distribution of {|X|} to find the variance

In[167]:=

d - expdev^2

Out[167]=

d - (2 Gamma[1/2 + d/2]^2)/Gamma[d/2]^2

In[168]:=

Plot[%, {d, 1, 10}] ;

[Graphics:HTMLFiles/index_56.gif]

The variance of {|X|} converges to 1/2 . Through Chebyshev ' s inequality, this means that at least 75% of the mass falls in the band between (d - 1/2)^(1/2) - 1 and (d - 1/2)^(1/2) + 1

Here's a numerical simulation to confirm this

In[203]:=

band[dd_] := Module[{ps}, ps = Sequence @@ ({area/(2Pi)^(d/2) Exp[-r^2/2], {r, (d - 1/2)^(1/2) - 1, (d - 1/2)^(1/2) + 1}}/.d→dd) ; NIntegrate[Evaluate[ps]]]

In[214]:=

Plot[band[d], {d, 2, 100}, PlotRange→ {0, 1}] ;

[Graphics:HTMLFiles/index_60.gif]

In[247]:=

band[1000000]

Out[247]=

0.842701

It seems to converge to a value larger than predicted by Chebychev bound (.75) but much smaller than predicted using normality assumption (.95)

Created by Mathematica  (May 18, 2006) Valid XHTML 1.1!