In a reply to Markiyan Hirnyk's recent question on this topic, Alec Mihailovs posted solutions

for n=1..7 for the mean distance between random points in an n-cube.

I also posted some code to simulate the value for arbitrary n.

In this post I'm going to show the results of simulations for n=1..10, along with histograms

of the results (with superimposed normal pdf curve) and normal probability plots, since the

simulations rely on the central limit theorem. % errors, for n=1..7 where the result is known,

and 95% confidence intervals for those where the result is not known (by me), are also shown.

ncubem:=proc(n::posint,N::posint,m)

# procedure simulates the mean value between two random points A and B

# in a n-cube of face length m, N times.

local A,B,s,i,RV,dim:

RV := Statistics:-RandomVariable(Uniform(0,m)):

s:=0:

for dim from 1 to n do

A[dim]:=Statistics:-Sample(RV,N);

B[dim]:=Statistics:-Sample(RV,N);

end do:

for i from 1 to N do

s:=s+sqrt(sum((A[dimm][i]-B[dimm][i])^2,dimm=1..n))

end do:

s/N;

end proc:

unitncubesim:=proc(n::posint,N::posint,monte::posint,exact)

# procedure samples the mean value of a unit n-cube N times

# exact is an expression for the known exact value. If this

# is provided, the statistical error is computet. If it is not

# then a 95% confidence interval for the mean is computed.

local sample,i,mu,sigma,qtl,ciL,ciU,p1,p2,p3,XX,A,pr1;

sample:=seq(ncubem(n,monte,1),i=0..N):

mu:=Statistics:-Mean([sample]);

sigma:=Statistics:-StandardDeviation([sample]):

printf("simulated value for %d-cube: %0.6f",n,mu);

if type(exact,algebraic) then

printf("; known value: %0.6f, %0.4f%% error",evalf(exact),100*evalf(abs(mu-exact)/exact,6));

else

qtl:=Statistics:-Quantile(Normal(0,1),evalf(1-((1-(95/100))/2))):

ciL:=evalf(mu-(qtl*sigma/sqrt(N)),5):ciU:=evalf(mu+(qtl*sigma/sqrt(N)),5):

printf(" with 95%% confidence interval [%0.6f,%0.6f]",ciL,ciU);

end if;

p1:=Statistics[Histogram]([sample],bincount=15):

XX:=Statistics:-RandomVariable(Normal(mu,sigma)):

p2:=plot(Statistics:-PDF(XX,x),x=mu-(4*sigma)..mu+(4*sigma)):

p1:=plots:-display(p1,p2);

p3:=Statistics:-ProbabilityPlot([sample],Normal(mu,sigma),thickness=2,color=BLACK):

print(plots:-display(Array([p1,p3])));

end proc:

f:=n->evalf(2^n*Int(sqrt(add((x||i)^2,i=1..n))*mul(1-x||i,i=1..n),

[seq(x||i=0..1,i=1..n)])):

#HT Alec Mihailovs for this solution for n=1..7

F:=seq(f(n),n=1..7):

Digits:=15:

t:=time():

unitncubesim(1,2000,500,F[1]);

simulated value for 1-cube: 0.333516; known value: 0.333333, 0.0549% error

unitncubesim(2,1500,500,F[2]);

simulated value for 2-cube: 0.521228; known value: 0.521405, 0.0339% error

unitncubesim(3,1000,500,F[3]);

simulated value for 3-cube: 0.661723; known value: 0.661707, 0.0024% error

unitncubesim(4,900,400,F[4]);

simulated value for 4-cube: 0.777830; known value: 0.777666, 0.0211% error

unitncubesim(5,800,300,F[5]);

simulated value for 5-cube: 0.878373; known value: 0.878531, 0.0180% error

unitncubesim(6,700,300,F[6]);

simulated value for 6-cube: 0.969059; known value: 0.968942, 0.0121% error

unitncubesim(7,650,300,F[7]);

simulated value for 7-cube: 1.052150; known value: 1.051584, 0.0542% error

unitncubesim(8,600,300,"");

simulated value for 8-cube: 1.128114 with 95% confidence interval [1.12700,1.12920]

unitncubesim(9,550,300,"");

simulated value for 9-cube: 1.200337 with 95% confidence interval [1.199100,1.201500]

unitncubesim(10,500,300,"");

simulated value for 10-cube: 1.267209 with 95% confidence interval [1.266000,1.268400]

time()-t:

%/60;

8.0751