Analysis of Keplers' platonic solid model

How wrong was he?

Copernican and Modern measurements come from: Orbital Mechanics: Theory and Applications By Tom Logsdon pg 6

Kepler's predictions come from: Wolfram Alpha: Kepler's Mysterium Cosmographicum



In [1]:

    
import pandas as pd



In [6]:

    
planets = ['Mercury', 'Venus', 'Earth', 'Mars', 'Jupiter', 'Saturn']

copernican_measurements = pd.Series(data = [.38, .72, 1, 1.52, 5.22, 9.07], index=planets)
platonic_predictions = pd.Series(data = [.408, .795, 1, 1.258, 3.775, 6.539], index=planets)
modern_measurements = pd.Series(data = [.39, .72, 1, 1.52, 5.2, 9.54], index=planets)

copernican_error = copernican_measurements - platonic_predictions
modern_error = modern_measurements - platonic_predictions



In [29]:

    
sum_linear_error = (platonic_predictions - copernican_measurements).sum()
print 'sum linear', sum_linear_error

mean_linear_error = (platonic_predictions - copernican_measurements).mean()
print 'mean linear', mean_linear_error

sum_sqared_error = ((platonic_predictions - copernican_measurements)**2).sum()
print 'sum squared', sum_sqared_error

mean_sqared_error = ((platonic_predictions - copernican_measurements)**2).mean()
print 'mean squared', mean_sqared_error

root_mean_sqared_error = sqrt(((platonic_predictions - copernican_measurements)**2).mean())
print 'root mean squared', root_mean_sqared_error









    



sum linear -4.135
mean linear -0.689166666667
sum squared 8.569039
mean squared 1.42817316667
root mean squared 1.19506199281



In [4]:

    
#copernicus vs modern
root_mean_sqared_error = sqrt(((copernican_measurements - modern_measurements)**2).mean())
print 'root mean squared', root_mean_sqared_error









    



root mean squared 0.192093727123



In [8]:

    
# compare kepler with copernicus, and modern measures

plt.figure(figsize=(10,10))

plt.subplot2grid(shape=(5,1), loc=(0,0), rowspan=4)
plt.box('off')
#plt.xticks([])
plt.ylabel('Mean Orbit Radius', fontsize=18)

x_cols = [-4.0, 0., 4.0]

#draw the boxes. Dont know why they wont show up.
for platonic, modern, error in zip(platonic_predictions, modern_measurements, modern_error):
    gca().add_patch(Rectangle((.5*x_cols[2]-.5*error,platonic), width=error, height=error, 
                            fill=True,  color='b', alpha=.5))

for platonic, modern, error in zip(platonic_predictions, copernican_measurements, copernican_error):
    gca().add_patch(Rectangle((.5*x_cols[0]-.5*error,platonic), width=error, height=error, 
                            fill=True,  color='r', alpha=.5))

arrowprops = dict(arrowstyle="-", linewidth = 0.75, color = 'black')
    
#mark the labels
for planet, radius in zip(copernican_measurements.index, copernican_measurements):
    plt.annotate(planet, xy=(.5*x_cols[0],radius), xytext=(x_cols[0],radius), 
                 fontsize=18, ha='center', va='center', 
                 arrowprops=arrowprops)

for planet, radius in zip(platonic_predictions.index, platonic_predictions):
    plt.annotate(planet, xy=(.5*x_cols[0],radius), xytext=(x_cols[1],radius), 
                 fontsize=18, ha='center', va='center',
                 arrowprops=arrowprops)

for planet, radius in zip(modern_measurements.index, modern_measurements):
    plt.annotate(planet, xy=(.5*x_cols[2],radius), xytext=(x_cols[2],radius), 
                 fontsize=18, ha='center', va='center',
                 arrowprops=arrowprops)


#gca().add_patch(Rectangle((1,1),1,1))    
    
plt.ylim(0,9)
plt.axis('equal')
plt.suptitle('Copernicus Measurement     Kepler Predictions     Modern Measurement', fontsize=18)

### summations
plt.subplot2grid(shape=(5,1), loc=(4,0), rowspan=1)
plt.box('off')
plt.xticks([])
plt.yticks([])
plt.plot([0,0],[1,1])
plt.ylabel('Sum', fontsize=18)









    Out[8]:





<matplotlib.text.Text at 0x1048335d0>



In [13]:

    
#Maybe now add the normal error measurements

plt.figure(figsize=(6,10))

ax1 = plt.subplot()
ax1.set_frame_on(False)
ax1.set_xticks([])

x_cols = [-6.5, 0., 5.0]

for platonic, modern, error in zip(platonic_predictions, copernican_measurements, copernican_error):
    ax1.add_patch(Rectangle((.5*x_cols[0]-.5*error,platonic), width=error, height=error, 
                            fill=True,  color='r', alpha=.5))

arrowprops = dict(arrowstyle="-", linewidth = 0.75, color = 'black')
    
#mark the labels
for planet, radius in zip(copernican_measurements.index, copernican_measurements):
    plt.annotate(planet+': '+str(radius), xy=(.5*x_cols[0],radius), xytext=(x_cols[0],radius), 
                 fontsize=18, ha='left', va='center', 
                 arrowprops=arrowprops)

for planet, radius in zip(platonic_predictions.index, platonic_predictions):
    plt.annotate(planet+': '+str(radius), xy=(.5*x_cols[0],radius), xytext=(x_cols[1],radius), 
                 fontsize=18, ha='right', va='center',
                 arrowprops=arrowprops)
    
for platonic, copernican in zip(platonic_predictions, copernican_measurements):
    plt.plot([.5*x_cols[0],.5*x_cols[0]], [platonic, copernican], 'k--' )


ax1.set_yticks([])
ax1.set_ylim(0,9)
ax1.axis('equal');



In [26]:

    
#Only compare kepler with copernicus

plt.figure(figsize=(6,10))
cla()

ax1 = plt.subplot()
ax1.set_frame_on(False)
#ax1.set_xticks([])
#ax1.set_yticks([])
sse = ax1.add_patch(Rectangle((0,1), width=2*root_mean_sqared_error, height=3*root_mean_sqared_error, 
                            fill=True,  color='r', alpha=.5))

mse = ax1.add_patch(Rectangle((0,1), width=root_mean_sqared_error, height=root_mean_sqared_error, 
                            fill=False,  color='k', alpha=1))

ax1.add_patch(Rectangle((0+root_mean_sqared_error,1), width=root_mean_sqared_error, height=root_mean_sqared_error, 
                            fill=False,  color='k', alpha=1))

ax1.add_patch(Rectangle((0,1+root_mean_sqared_error), 
                        width=root_mean_sqared_error, height=root_mean_sqared_error, 
                            fill=False,  color='k', alpha=1))

ax1.add_patch(Rectangle((0+root_mean_sqared_error,1+root_mean_sqared_error), 
                        width=root_mean_sqared_error, height=root_mean_sqared_error, 
                            fill=False,  color='k', alpha=1))


ax1.add_patch(Rectangle((0,1+2*root_mean_sqared_error), width=root_mean_sqared_error, height=root_mean_sqared_error, 
                            fill=False,  color='k', alpha=1))

ax1.add_patch(Rectangle((0+root_mean_sqared_error,1+2*root_mean_sqared_error), 
                        width=root_mean_sqared_error, height=root_mean_sqared_error, 
                            fill=False,  color='k', alpha=1))

rmse = plt.plot([0,root_mean_sqared_error],[.9,.9], '-b', marker='o')

plt.legend((sse, mse, rmse), 
           ('Sum of Squared Error', 'Mean Squared Error', 'Root Mean Squared Error'))



plt.axis('equal')
plt.ylim(0,9)









    



/Users/houghtons/Downloads/matplotlib/lib/matplotlib/legend.py:613: UserWarning: Legend does not support [<matplotlib.lines.Line2D object at 0x104f7ac10>]
Use proxy artist instead.

http://matplotlib.sourceforge.net/users/legend_guide.html#using-proxy-artist

  (str(orig_handle),))






    Out[26]:





(0, 9)

Theil Statistics:

note that the way scipy calculates the std dev is different from that used in stermans theil statistics paper. Its irrelevant with large datasets, but with small data, makes a difference.



In [33]:

    
#calculate bias component

bias = ((platonic_predictions.mean() - copernican_measurements.mean())**2)/mean_sqared_error
bias









    Out[33]:





0.33255819779401913



In [38]:

    
#calculate the variance component

# the standard deviation in scipy uses a different normalization than we want, so calculate our own.
std_platonic = sqrt(((platonic_predictions-platonic_predictions.mean())**2).mean())
std_copernican = sqrt(((copernican_measurements-copernican_measurements.mean())**2).mean())

#variance component
variance = (std_platonic - std_copernican)**2/mean_sqared_error
variance









    Out[38]:





0.66341051413639607



In [42]:

    
# calculate the correlation component

corrcoef = np.corrcoef(platonic_predictions, copernican_measurements).min()

correlation = 2*(1-corrcoef)*std_platonic*std_copernican
correlation









    Out[42]:





0.0057573774480853511



In [43]:

    
#check that they sum (close) to 1
variance + bias + correlation









    Out[43]:





1.0017260893785005

plot bias error



In [46]:

    
#Maybe now add the normal error measurements

plt.figure(figsize=(6,10))

ax1 = plt.subplot()
ax1.set_frame_on(False)
ax1.set_xticks([])

x_cols = [-6.5, 0., 5.0]

arrowprops = dict(arrowstyle="-", linewidth = 0.75, color = 'black', alpha=0)
    
#mark the labels
for planet, radius in zip(copernican_measurements.index, copernican_measurements):
    plt.annotate(planet+': '+str(radius), xy=(.5*x_cols[0],radius), xytext=(x_cols[0],radius), 
                 fontsize=18, ha='left', va='center', alpha=.5, 
                 arrowprops=arrowprops)

for planet, radius in zip(platonic_predictions.index, platonic_predictions):
    plt.annotate(planet+': '+str(radius), xy=(.5*x_cols[0],radius), xytext=(x_cols[1],radius), 
                 fontsize=18, ha='right', va='center', alpha=.5,
                 arrowprops=arrowprops)
    
#leave this in a sa kludgey way to keep the box the right shape    
for platonic, copernican in zip(platonic_predictions, copernican_measurements):
    plt.plot([.5*x_cols[0],.5*x_cols[0]], [platonic, copernican], 'k--', alpha=0 )

arrowprops = dict(arrowstyle="-", linewidth = 0.75, color = 'black', alpha=1)

plt.annotate('mean', xy=(.5*x_cols[0],copernican_measurements.mean()), 
             xytext=(x_cols[0],copernican_measurements.mean()), 
                 fontsize=18, ha='left', va='center', alpha=1,
                 arrowprops=arrowprops)

plt.annotate('mean', xy=(.5*x_cols[0],platonic_predictions.mean()), 
             xytext=(x_cols[1],platonic_predictions.mean()), 
                 fontsize=18, ha='right', va='center', alpha=1,
                 arrowprops=arrowprops)
    
error = platonic_predictions.mean() - copernican_measurements.mean()
ax1.add_patch(Rectangle((.5*x_cols[0]-.5*error,copernican_measurements.mean()), width=error, height=error, 
                            fill=True,  color='b', alpha=.5))    
    
ax1.set_yticks([])
ax1.set_ylim(0,9)
ax1.axis('equal')









    Out[46]:





(-3.594583333333333, -2.9054166666666665, 0.38, 9.0700000000000003)

Plot the variance component



In [50]:

    
#Maybe now add the normal error measurements

plt.figure(figsize=(6,10))

ax1 = plt.subplot()
ax1.set_frame_on(False)
ax1.set_xticks([])

x_cols = [-6.5, 0., 5.0]


arrowprops = dict(arrowstyle="-", linewidth = 0.75, color = 'black', alpha=0)
    
#mark the labels
for planet, radius in zip(copernican_measurements.index, copernican_measurements):
    plt.annotate(planet+': '+str(radius), xy=(.5*x_cols[0],radius), xytext=(x_cols[0],radius), 
                 fontsize=18, ha='left', va='center', alpha=.25, 
                 arrowprops=arrowprops)

for planet, radius in zip(platonic_predictions.index, platonic_predictions):
    plt.annotate(planet+': '+str(radius), xy=(.5*x_cols[0],radius), xytext=(x_cols[1],radius), 
                 fontsize=18, ha='right', va='center', alpha=.25,
                 arrowprops=arrowprops)
    
#leave this in as a kludgey way to keep the box the right shape    
for platonic, copernican in zip(platonic_predictions, copernican_measurements):
    plt.plot([.5*x_cols[0],.5*x_cols[0]], [platonic, copernican], 'k--', alpha=0 )

arrowprops = dict(arrowstyle="-", linewidth = 0.75, color = 'black', alpha=0.25)

plt.annotate('mean', xy=(.7*x_cols[0],copernican_measurements.mean()), 
             xytext=(x_cols[0],copernican_measurements.mean()), 
                 fontsize=18, ha='left', va='center', alpha=1,
                 arrowprops=arrowprops)

plt.annotate('mean', xy=(.3*x_cols[0],platonic_predictions.mean()), 
             xytext=(x_cols[1],platonic_predictions.mean()), 
                 fontsize=18, ha='right', va='center', alpha=1,
                 arrowprops=arrowprops)


plt.annotate('+1 std', xy=(.3*x_cols[0],platonic_predictions.mean()+std_platonic), 
             xytext=(x_cols[1],platonic_predictions.mean()+std_platonic), 
                 fontsize=18, ha='right', va='center', alpha=1,
                 arrowprops=arrowprops)

plt.annotate('+1 std', xy=(.7*x_cols[0],copernican_measurements.mean()+std_copernican), 
             xytext=(x_cols[0],copernican_measurements.mean()+std_copernican), 
                 fontsize=18, ha='left', va='center', alpha=1,
                 arrowprops=arrowprops)


#copernican std dev
plt.plot([.7*x_cols[0],.7*x_cols[0]], [copernican_measurements.mean()-std_copernican,
                                       copernican_measurements.mean()+std_copernican], 
         'k-', alpha=.25, linewidth=4)

plt.plot([.3*x_cols[0],.3*x_cols[0]], [platonic_predictions.mean()-std_platonic,
                                       platonic_predictions.mean()+std_platonic], 
         'k-', alpha=.25, linewidth=4)

#black copernican std dev
plt.plot([.65*x_cols[0],.65*x_cols[0]], [5, 5+std_copernican], 
         'k-', alpha=1, linewidth=4)

# black platonic std dev
plt.plot([.35*x_cols[0],.35*x_cols[0]], [5, 5+std_platonic], 
         'k-', alpha=1, linewidth=4)

# square error
error = std_platonic-std_copernican
ax1.add_patch(Rectangle((.5*x_cols[0]-.5*error,5+std_copernican), width=error, height=error, 
                            fill=True,  color='g', alpha=.5))  

# guide lines
plt.plot([.7*x_cols[0],.65*x_cols[0]], [copernican_measurements.mean(), 5], 
         'k-', alpha=.25, linewidth=2)

plt.plot([.7*x_cols[0],.65*x_cols[0]], [copernican_measurements.mean()+std_copernican, 5+std_copernican], 
         'k-', alpha=.25, linewidth=2)

plt.plot([.3*x_cols[0],.35*x_cols[0]], [platonic_predictions.mean(), 5], 
         'k-', alpha=.25, linewidth=2)

plt.plot([.3*x_cols[0],.35*x_cols[0]], [platonic_predictions.mean()+std_platonic, 5+std_platonic], 
         'k-', alpha=.25, linewidth=2)

plt.plot([.65*x_cols[0],.5*x_cols[0]], 
         [5+std_copernican, 5+std_copernican], 
         'k-', alpha=.25, linewidth=1)

plt.plot([.35*x_cols[0],.5*x_cols[0]], 
         [5+std_platonic, 5+std_platonic], 
         'k-', alpha=.25, linewidth=1)

ax1.set_yticks([])
ax1.set_ylim(0,9)
ax1.axis('equal');

Show how each error relates to the mse



In [52]:

    
#Draw the boxes for the side

plt.figure(figsize=(6,10))
cla()

ax1 = plt.subplot()
ax1.set_frame_on(False)

ax1.add_patch(Rectangle((0,1), width=root_mean_sqared_error, height=root_mean_sqared_error, 
                            fill=False,  color='k', alpha=1))


bias_error = (platonic_predictions.mean() - copernican_measurements.mean())**2

ax1.add_patch(Rectangle((0,1), width=bias_error/root_mean_sqared_error, height=root_mean_sqared_error, 
                            fill=True,  color='b', alpha=.5))

variance_error = (std_platonic-std_copernican)**2

# comment this out to just get the bias error in the box
ax1.add_patch(Rectangle((0+bias_error/root_mean_sqared_error,1), 
                        width=variance_error/root_mean_sqared_error, height=root_mean_sqared_error, 
                            fill=True,  color='g', alpha=.5))

plt.axis('equal')
plt.ylim(0,9);



In [ ]: