Chapter 1, Section 1.5

Enter the orbital periods of the five major uranian satellites (Miranda, Ariel, Umbriel, Titania and Oberon).

Uperiod = {1.413, 2.520, 4.144, 8.706, 13.463};

Take the logarithm (base 10) of each period.

LUperiod = Log[10,Uperiod]

Carry out a linear least squares fit to the data using Mathematica's Fit command.  Here i=1 denotes Miranda, i=2 is Ariel, i=3 is Umbriel, i=4 is Titania and i=5 is Oberon.

fn = Fit[LUperiod, {1,i}, {i}]

Now calculate values of this function for each value of i.

calc = Table[fn, {i,1,5}]

Compare these with the observed values to calculate the root mean square deviation for the actual uranian system.

χUranus = Sqrt[Apply[Plus,(LUperiod-calc)^2]/5]

Now define a function CalculateChi to find the equivalent value of χ for any list of five periods given in period.

CalculateChi[period_]:=
Module[{fun, calc, lper=Log[10,period]},
fun = Fit[lper, {1,i},{i}];
calc = Table[fun, {i,5}];
Sqrt[Apply[Plus, (lper-calc)^2]/5]
]

Set the first orbital period to be that of Miranda.

first = Uperiod[[1]]

Set a seed for use in the random number generator.  This allows the same values to be reproduced.

SeedRandom[257];

Define a function PeriodGen that will generate a set of five orbital periods using fixed lower and upper limits for successive ratios of periods.

PeriodGen[lower_, upper_]:=
Module[{range = upper-lower},
NestList[(lower+Random[]*range)*#&, first, 4]
]

Define a function GenerateChi that will use CalculateChi to calculate the value of χ for the random orbital periods produced by PeriodGen.

GenerateChi[n_, lower_, upper_]:=
Map[CalculateChi[#]&, Table[PeriodGen[lower,upper],{n}]]

For example, calculate five orbital periods, starting with Miranda's period and with successive ratios lying between 1.546 and 2.101 ...

PeriodGen[1.546, 2.101]

... and calculate χ for this set of periods.

CalculateChi[%]

As a further example, generate 10 sets of orbital periods with lower=1.0 and upper=1.2, calculating χ for each:

GenerateChi[10, 1.0, 1.2]

Now calculate the number of sets of five orbital periods out of (10^5)with lower=1.546 and upper=2.101 (the observed values for the uranian system) that have χ ≤ χUranus.  NB: This can take a while!

χsmaller=
Length[Select[
GenerateChi[100000, 1.546, 2.101], (#<=χUranus)&]]

Hence, an estimate of the probability of generating as good a fit as the uranian system by chance is:

probability=N[χsmaller/100000]

This provides the "+" point in Fig.1.4.