:

## Simulation of mean distance between random points in a hypercube

Maple
`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..7F:=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`
` `

﻿