39.1: LSF Example - Tracking the Planets

The following table lists the average distance from the sun to each of the first seven planets, using Earth’s distance as a unit of measure (AUs).

Mercury	Venus	Earth	Mars	Jupiter	Saturn	Uranus
0.39	0.72	1.00	1.52	5.20	9.54	19.2

The following is a plot of the data:

# Here are some libraries you may need to use

%matplotlib inline
import matplotlib.pylab as plt
import numpy as np
import sympy as sym
import math
sym.init_printing()

distances = [0.39, 0.72, 1.00, 1.52, 5.20, 9.54, 19.2]
planets = ['Mercury', 'Venus', 'Earth', 'Mars', 'Jupiter', 'Satern','Uranus']
ind = [1.0,2.0,3.0,4.0,5.0,6.0,7.0]

plt.scatter(ind, distances);
plt.xticks(ind,planets)
plt.ylabel('Distance (AU)')

Note

That the above plot does not look like a line, and so finding the line of best fit is not fruitful. It does, however look like an exponential curve (maybe a polynomial?). The following step transforms the distances using the numpy log function and generates a plot that looks much more linear.

log_distances = np.log(distances)

plt.scatter(ind,log_distances)
plt.xticks(ind,planets)
plt.ylabel('Distance (log(AU))')

For this question we are going to find the coefficients (\(c\)) for the best fit line of the form \(c_1 + c_2 i = \log d\), where \(i\) is the index of the planet and \(d\) is the distance.

The following code constructs this problem in the form \(Ax=b\) and define the \(A\) matrix and the \(b\) matrix as numpy matrices

A = np.matrix(np.vstack((np.ones(len(ind)),ind))).T
b = np.matrix(log_distances).T
sym.Matrix(A)

sym.Matrix(b)

Do This

Solve for the best fit of \(Ax=b\) and define a new variable \(c\) which consists of the of the two coefficients used to define the line (\(\log d = c_1 + c_2 i\))

Do This

Modify the following code (as needed) to plot your best estimates of \(c_1\) and \(c_2\) against the provided data.

## Modify the following code

est_log_distances = (c[0] + c[1]*np.matrix(ind)).tolist()[0]
plt.plot(ind,est_log_distances)

plt.scatter(ind,log_distances)
plt.xticks(ind,planets)
plt.ylabel('Distance (log(AU))')

We can determine the quality of this line fit by calculating the root mean squared error between the estimate and the actual data:

rmse = np.sqrt(((np.array(log_distances)  - np.array(est_log_distances)) ** 2).mean())
rmse

Finally, we can also make the plot on the original axis using the inverse of the log (i.e. the exp function):

est_distances = np.exp(est_log_distances)
plt.scatter(ind,distances)
plt.plot(ind,est_distances)
plt.xticks(ind,planets)
plt.ylabel('Distance (AU)')

The asteroid belt between Mars and Jupiter is what is left of a planet that broke apart. Let’s the above calculation again but renumber so that the index of Jupyter is 6, Saturn is 7 and Uranus is 8 as follows:

distances = [0.39, 0.72, 1.00, 1.52, 5.20, 9.54, 19.2]
planets = ['Mercury', 'Venus', 'Earth', 'Mars', 'Jupiter', 'Satern','Uranus']
ind = [1,2,3,4,6,7,8]

log_distances = np.log(distances)

Do This

Repeat the calculations from above with the updated model. Plot the results and compare the RMSE.

## Copy and Paste code from above

## Copy and Paste code from above
est_log_distances = (c[0] + c[1]*np.matrix(ind)).tolist()[0]

est_distances = np.exp(est_log_distances)
plt.scatter(ind,distances)
plt.plot(ind,est_distances)
plt.xticks(ind,planets)
plt.ylabel('Distance (AU)')

rmse = np.sqrt(((np.array(log_distances)  - np.array(est_log_distances)) ** 2).mean())
rmse
## Copy and Paste code from above

This model of planet location was used to help discover Neptune and prompted people to look for the “missing planet” in position 5 which resulted in the discovery of the asteroid belt. Based on the above model, what is the estimated distance of the asteroid belt and Neptune (index 9) from the sun in AUs?

Hint

You can check your answer by searching for the answer on-line.

#Put your prediction calcluation here