<jupyter_start><jupyter_text>Regression and the Age of the universe    Import scientific python packages<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
import pandas as pd<jupyter_output><empty_output><jupyter_text>Age of the universeThe inspiration for today's exercise fame from Lecture 16 of Prof. <NAME>'s [Python for Earth Science Students](https:github.com/ltauxe/Python-for-Earth-Science-Students) class and some of our materials is modified from that lecture. That whole course is open source and pretty great, if you want a deeper dive into python programming. Linear regression is a method for estimating the potential association between two variables. Today, we will use the retreat velocity of galaxies and supernova as a function of their distance as our example data set. Such data underlies what has come to be known as "Hubble's Law" (same Hubble as for the Hubble telescope). Hubble published these results in 1929 [Hubble, E. P. (1929) Proc. Natl. Acad. Sci., 15, 168–173.]  At the time, it was unclear whether the universe was static, expanding, or collapsing. Hubble hypothesized that if the universe were expanding, then everything in it would be moving away from us. The greater the distance between the Earth and the galaxy, the faster it must be moving.  So all that had to be done was to measure the distance and velocity of distant galaxies.  Easy-peasy - right?  To measure velocity, Hubble made use of the doppler shift. To understand how this works, recall that the pitch you hear as an ambulance approaches changes. During doppler shift, the ambulance's pitch changes from high (as it approaches) to low (as it recedes). The pitch changes because the relative frequency of the sound waves changes. The frequency increases as the ambulance approaches, leading to a higher pitch, and then decreases as it moves away, resulting in a lower pitch.  Just in case you haven't had this life experience, let's listen to such a [siren](https://www.youtube.com/watch?v=imoxDcn2Sgo) here.The same principle applies to light, but rather than hear a change in frequency, we observe a shift in the wavelength (the color) emitted by the galaxy. If a star or galaxy is moving away from us, its absorption bands are shifted towards longer wavelengths - the red end of the visible spectrum. The faster the star or galaxy travels away from the observer, the greater the shift will be to the red:So a star (or galaxy) moving away from us will have a red shift with the wavelength being spread out._[Figures from http://www.a-levelphysicstutor.com/wav-doppler.php](http://www.a-levelphysicstutor.com/wav-doppler.php)_Hubble measured the red shift of different galaxies and converted them to velocities. He then estimated the distance to these objects, which is harder to do (and he was pretty far off). Improving such data was a major motivation of the Hubble Space Telescope. Those data and continued improvement to approaches for estimating these distances and velocities and investigating additional types of celestial objects is a major focus of ongoing research.   Type 1a supernovae dataLet's import data from Freedman et al. (2000) of the distance and retreat velocity of type 1a supernovae. These supernovae are described as follows in a review paper that Freedman wrote in 2010 (https://doi.org/10.1146/annurev-astro-082708-101829):> One of the most accurate means of measuring cosmological distances out into the Hubble flowutilizes the peak brightness of SNe Ia. The potential of supernovae for measuring distances wasclear to early researchers (e.g., Baade, Minkowski, Zwicky), but it was the Hubble diagram ofKowal (1968) that set the modern course for this field, followed by decades of work by Sandage,Tammann, and collaborators (e.g., Sandage & Tammann 1982, 1990; see also the review byBranch 1998). Analysis by Pskovskii (1984), followed by Phillips (1993), established a correlationbetween the magnitude of a SN Ia at peak brightness and the rate at which it declines, thusallowing supernova luminosities to be “standardized.” This method currently probes farthestinto the unperturbed Hubble flow, and it possesses very low intrinsic scatter:*Freedman and Madore (2010) who then go onto describe how using Cepheid variable stars (a type of pulsating star) has allowed for the distances to be better calibrated.*> SNe Ia result from the thermonuclear runaway explosions of stars.From observations alone, the presence of SNe Ia in elliptical galaxies suggests that they do notcome from massive stars. Many details of the explosion are not yet well understood, but thegenerally accepted view is that of a carbon-oxygen, electron-degenerate, nearly-Chandrasekharmasswhite dwarf orbiting in a binary system with a close companion *Freedman and Madore (2010)*<jupyter_code>Supernova_data = pd.read_csv('Data/Freedman2000_Supernova1a.csv')
Supernova_data.tail()<jupyter_output><empty_output><jupyter_text>The ```VCMB``` column is velocity relative to the cosmic microwave background in km s$^{-1}$. The ```D(Mpc)``` column is the distance in Mpc which is the unit typically used for these measurements. 1 Mpc =  3.09 x 10$^{19}$ kmGo ahead and double-click on this cell to see how I am getting labels that have the proper superscripts.To create nice labels with superscripts, we can use latex formatting, which can also be done in a markdown cell.  For a superscript, first we need to encase the text in dollar signs and then use the ^ symbol to make the following text a superscript. If there is more than one number in the superscript, you must enclose what you want as the superscript in curly braces.For example, to print $10^3$, we use ```$10^3$``` and for 'per second' (s$^{-1}$): ```s$^{-1}$```<jupyter_code>plt.scatter(Supernova_data[...],Supernova_data[...],color='red',label='1A Supernovae data (Freedman et al. 2000)')
plt.ylabel('Velocity (km s$^{-1}$)')
plt.xlabel('Distance (Mpc)')
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Determining the slope of this line (the Hubble constant)We have distance on the x-axis in megaparsecs and velocity on the y-axis in km/s. The slope of this line is the Hubble constant: $v = H_o d$where $v$ is velocity, $d$ is distance, and $H_o$ is the Hubble constant. This looks a lot like the equation for a line through the data ($y=mx + b$) where $m$ is the slope and $b$ is the y-intercept.  In this case, the y-intercept should be 0 or nearly so, and $m$ is $H_o$.So how do we find the slope?Here is where we can use linear regression to find the "best fit" line through the data. The approach is to minimize the sum of the squares of the distances (residuals) between the points and a line through them. In this illustration below, the residuals are the vertical distance between each data point and the line:The approach in linear regression is to find the line that minimizes the squared value of these distances all added up. We determine the best-fit line through this least squares approach using the ```np.polyfit()``` function. A straight line is a first degree polynomial (*note that the function can fit higher order polynomials as well*).<jupyter_code>np.polyfit?<jupyter_output><empty_output><jupyter_text>Fitting a line with `np.polyfit()``np.polyfit()` can be used to calculate best fit lines (setting the degree (```deg```) to 1), or higher order curves (setting degree to 2 or higher) returning the slope and the intercept. Let's put it to use:<jupyter_code>np.polyfit(...)<jupyter_output><empty_output><jupyter_text>So $H_o$, the slope of the best-fit line, is 67.5 (in the odd units of kilometers per second per megaparsec).  Let's plot the best fit line on our graph.  We can assign the best fitting slope and y-intercept from **np.polyfit( )** to a variable (**m_b**).<jupyter_code>m_b= np.polyfit(Supernova_data['D(Mpc)'],Supernova_data['VCMB'],1)
print(m_b) #see if that worked<jupyter_output><empty_output><jupyter_text>**m_b**  seems to be an array of coefficients, where the first is the slope and the second is the y-intercept.  We can now use the function `np.polyval()` which will calculate new y values using the model of a linear fit. We can feed **m_b** into **np.polyval( )**, along with our x array to get a new set of y values which are the y values for the best-fit linear model. Then we can plot the model data as a black line along with the original data.<jupyter_code>model_y_values = np.polyval(m_b,Supernova_data['D(Mpc)'])
model_y_values
plt.scatter(Supernova_data['D(Mpc)'],Supernova_data['VCMB'],
            color='red',label='1A Supernovae data (Freedman et al. 2000)')
plt.scatter(Supernova_data['D(Mpc)'],model_y_values,
         color='black',marker='s',label='best fit to 1A Supernovae data')
plt.ylabel('Velocity (km s$^{-1}$)')
plt.xlabel('Distance (Mpc)')
plt.legend()
plt.show()
plt.scatter(Supernova_data['D(Mpc)'],Supernova_data['VCMB'],
            color='red',label='1A Supernovae data (Freedman et al. 2000)')
plt.plot(Supernova_data['D(Mpc)'],model_y_values,
         color='black',label='best fit to 1A Supernovae data')
plt.ylabel('Velocity (km s$^{-1}$)')
plt.xlabel('Distance (Mpc)')
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Using this linear model for predictionWhat would we predict that the velocity would be for a supernova that happened to be 350 Mpc?**_Code for you to write_****Use the ```np.polyval()``` function to come up with what the linear model predicts the velocity would be?**  Evaluating model fitWe'd also like to know who well this model fits our data (i.e. how correlated the data are). We'll use the $R^{2}$ correlation coefficient for this. $R^{2}$ is zero for uncorrelated data, and 1 for perfectly linear data (so no misfit between the model line and data). We'll use the scipy function `stats.linregress` to compute $R^{2}$.<jupyter_code>from scipy.stats import stats<jupyter_output><empty_output><jupyter_text>And use it, to get what is normally called the $R^2$ value, which when 1. represents perfect agreement.> Pearson correlation coefficient between several example X,Y sets. Source: https://en.wikipedia.org/wiki/Correlation_and_dependence<jupyter_code>stats.linregress(Supernova_data['D(Mpc)'],Supernova_data['VCMB'])<jupyter_output><empty_output><jupyter_text>Fitting a line with `stats.linregress`<jupyter_code>slope, intercept, rvalue, pvalue, stderr = stats.linregress(Supernova_data['D(Mpc)'],Supernova_data['VCMB'])
rvalue**2<jupyter_output><empty_output><jupyter_text>Not a bad fit!  We can have confidence that there is a strong correlation between distance and velocity. The universe is expanding.     Evaluting the fit through plotting residualsTo see how well the regression performs, the data scientist must measure how far off the estimates are from the actual values. These differences are called *residuals*.$$\mbox{residual} ~=~ \mbox{observed value} ~-~ \mbox{regression estimate}$$A residual is what's left over – the residue – after estimation. Residuals are the vertical distances of the points from the regression line. There is one residual for each point in the scatter plot. The residual is the difference between the observed value of $y$ and the fitted value of $y$, so for the point $(x, y)$,$$\mbox{residual} ~~ = ~~ y ~-~\mbox{fitted value of }y~~ = ~~ y ~-~\mbox{height of regression line at }x$$<jupyter_code>residual =
plt.scatter(Supernova_data['D(Mpc)'],residual,color='red')
plt.hlines(0,xmin=0,xmax=500)
plt.xlim(0,500)
plt.show()<jupyter_output><empty_output><jupyter_text>**The residual plot of a good regression shows no pattern. The residuals look about the same, above and below the horizontal line at 0, across the range of the predictor variable.**  Estimating the age of the universeTo calculate the age of the universe, we can use Hubble's law:  We had $v=H_o d$ as Hubble's law and we know that distance = velocity x time, or,  $d=vt$.  So, if we divide both sides by $v$ and  we get: 1=$H_o$t. Solving for $t$ (the age of the universe), we get $t=1/H_o$ [in some weird units.]<jupyter_code>t = 
print(t)<jupyter_output><empty_output><jupyter_text>But the units are weird (not years, Mpc s/km).  To fix this, we need to know how many kilometers are in a megaparsec.  As it happens,  there are 3.09 x 10$^{19}$ km/Mpc.So, we can calculate the age of the universe in seconds (**Age_sec**) by converting the megaparsecs to km:Age (s) = $t \frac{s \cdot Mpc}{km}$ x $3.09 x 10^{19} \frac {km}{Mpc}$<jupyter_code>Age_sec=
print(Age_sec)<jupyter_output><empty_output><jupyter_text>That's a lot of seconds!  We should convert seconds to years.  Here's another fun fact: there are approximately   $\pi$ x 10$^7$ seconds in a year. More precisely, there are 60 (s/min) x 60 (min/hr) x 24 (hr/day) x 365.25 (days/yr)<jupyter_code>s_yr=
print('%e'%(s_yr))<jupyter_output><empty_output><jupyter_text>Ok.  so not exactly $\pi \times 10^7$, but close....<jupyter_code>Age_yrs=
print(Age_yrs)<jupyter_output><empty_output><jupyter_text>And now in billions of years:<jupyter_code>print ('Age of the universe (in billions of years):')
print(Age_yrs*1e-9)<jupyter_output><empty_output><jupyter_text>**_Code for you to write_****Write a function that takes in a Hubble constant value and calculates the age of the Universe in billions of year**<jupyter_code>def age_of_universe(Hubble_constant):
    
    return age<jupyter_output><empty_output><jupyter_text>Using other data sets to estimate the Hubble constantDetermining the Hubble constant continues to be a major avenue of astrophysical research. In fact, Wendy Freedman's group just published a new study (https://arxiv.org/abs/1907.05922) that is summarized in this short video:https://www.youtube.com/watch?v=awcnVykOKZYFrom that paper here is a visualization of Hubble constant determinations over the past 18 years:Let's look at another data set from the 2000 study to see how different data sets can lead to different answers. Tully-Fisher Relation galaxy data> The total luminosity of a spiral galaxy (corrected to face-on inclination to account for extinction)is strongly correlated with the galaxy’s maximum (corrected to edge-on inclination) rotationvelocity. This relation, calibrated via the Leavitt Law or TRGB, becomes a powerful means of determiningextragalactic distances (Tully&Fisher 1977, Aaronson et al. 1986, Pierce&Tully 1988,Giovanelli et al. 1997). The TF relation at present is one of the most widely applied methods fordistance measurements *Freedman and Madore (2010)* **_Code for you to write_****Import the 'Data/Freedman2000_IBandTullyFisher.csv' file. Make a linear fit to determine the slope between `VCMB` and `D(Mpc)`. Calculate the implied age of the universe from these TF galaxy data alone.**  Going even further out into the universeLet's look at new data sets available for the classic Hubble problem.  I found one published by Betoule et al. in 2014 [http://dx.doi.org/10.1051/0004-6361/201423413](http://dx.doi.org/10.1051/0004-6361/201423413).   In this paper, data are plotted using the parameters $z$ and $\mu$ which are related to the red shift velocity and distance.  $z$ is the fractional shift in the spectral wavelength and $\mu$ is related to distance. Here is a plot from the Betoule et al. paper: _[Figure from Betoule et al., 2014.]  These data are type Ia supernova from different observation collaborations_  Notice that they plotted the data on a log scale. (This hides some surprising things.)It turns out that we have been looking at data that are low-z (that is relatively close and low red shift). We  need to convert $z$ and $\mu$ to distance and velocity to compare to the results we have considered thus far.  According to [http://hyperphysics.phy-astr.gsu.edu/hbase/Astro/hubble.html](http://hyperphysics.phy-astr.gsu.edu/hbase/Astro/hubble.html)velocity $v$ (as fraction of the speed of light, $c$) is given by${v\over c}= \bigl({{(z+1)^2-1}  \over {(z+1)^2+1}}\bigr)$where $c=3 \times 10^8$m s$^{-1}$.  And according to the Betoule et al. (2014) paper, $\mu$ relates to distance in parsecs $d$ like this:  $\mu=5\log(d/10)$.   Let's read in the data (available from this website:  http://cdsarc.u-strasbg.fr/viz-bin/qcat?J/A+A/568/A22sRM2.2), which are averages of the data shown in the figure above,and take a peek.<jupyter_code>Betoule_data = pd.read_csv('Data/mu_z.csv',header=1)
Betoule_data.head()<jupyter_output><empty_output><jupyter_text>Now we can plot it the same way as the cosmologists did in the paper, using $\mu$ and $\log z$:<jupyter_code>plt.scatter(Betoule_data[...],Betoule_data[...],color='blue')
plt.xlabel('z')
plt.ylabel('$\mu$')
plt.semilogx()
plt.show()<jupyter_output><empty_output><jupyter_text>To compare these new data with the previous considered data, we must do the following:  - Transform $z$  to velocity  - Transform  $\mu$ to distance using the equations provided. - Truncate the new dataset which goes to much farther distances than the 'old' data set<jupyter_code>c = 2.9979e8 / 1000 # speed of light in km/s
# the formula for v from z
Betoule_data['velocity'] = c * (((Betoule_data['z']+1.)**2-1.)/((Betoule_data['z']+1.)**2+1.)) 
Betoule_data['distance']=10000*(10.**((Betoule_data['mu'])/5.))*1e-9 # convert mu to Gpc
plt.figure(figsize=(8,6))
plt.scatter(Betoule_data['distance'],Betoule_data['velocity'],
         color='blue',label='1A Supernovae data (Betoule et al. 2014)')
plt.scatter(Supernova_data['D(Mpc)'],Supernova_data['VCMB'],
            color='red',label='1A Supernovae data (Freedman et al. 2000)')
plt.ylabel('Velocity (km s$^{-1}$)')
plt.xlabel('Distance (Mpc)')
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>These data sets are similar to one another for the "close" objects, but we can see that a linear model doesn't work well for objects that are at greater distances.To visualize this reality, let's plot the fit to the Freedman et al. 2000 data atop this plot (applying it to the Betoule distances using `np.polyval()`.<jupyter_code>model_y_values = np.polyval(m_b,Betoule_data['distance'])

plt.figure(figsize=(8,6))
plt.scatter(Betoule_data['distance'],Betoule_data['velocity'],
         color='blue',label='1A Supernovae data (Betoule et al. 2014)')
plt.scatter(Supernova_data['D(Mpc)'],Supernova_data['VCMB'],
            color='red',label='1A Supernovae data (Freedman et al. 2000)')
plt.plot(Betoule_data['distance'],model_y_values,
         color='black',label='1A Supernovae fit to Freedman data')
plt.ylabel('Velocity (km s$^{-1}$)')
plt.xlabel('Distance (Mpc)')
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Clearly this fit is quite poor.Let's make a first-order polynomial fit to all the Betoule data and then plot the residual:<jupyter_code>fit = np.polyfit(Betoule_data['distance'],Betoule_data['velocity'],1) 
y_values = np.polyval(fit,Betoule_data['distance'])

plt.subplot(2,1,1)
plt.scatter(Betoule_data['distance'],Betoule_data['velocity'])
plt.plot(Betoule_data['distance'],y_values,color='orange',)
plt.title('data and a polynomial degree 1 fit')
plt.ylabel('Velocity (km s$^{-1}$)')
plt.xlabel('Distance (Mpc)')

plt.subplot(2,1,2)
plt.scatter(Betoule_data['distance'],Betoule_data['velocity']-y_values)
plt.title('residuals of a polynomial degree 1 fit')
plt.ylabel('Residual velocity (km s$^{-1}$)')
plt.xlabel('Distance (Mpc)')

plt.tight_layout()
plt.show()<jupyter_output><empty_output><jupyter_text>There is a lot of structure to the residual of this degree 1 fit. Let's try a degree 2 polynomial fit (known as quadratic):$f(x)=ax^2+bx+c$<jupyter_code>fit = np.polyfit(Betoule_data['distance'],Betoule_data['velocity'],2) 
y_values = np.polyval(fit,Betoule_data['distance'])

plt.subplot(2,1,1)
plt.scatter(Betoule_data['distance'],Betoule_data['velocity'])
plt.plot(Betoule_data['distance'],y_values,color='orange',)
plt.title('data and a polynomial degree 2 fit')
plt.ylabel('Velocity (km s$^{-1}$)')
plt.xlabel('Distance (Mpc)')

plt.subplot(2,1,2)
plt.scatter(Betoule_data['distance'],Betoule_data['velocity']-y_values)
plt.title('residuals of a polynomial degree 2 fit')
plt.ylabel('Residual velocity (km s$^{-1}$)')
plt.xlabel('Distance (Mpc)')

plt.tight_layout()
plt.show()<jupyter_output><empty_output><jupyter_text>There is a lot of structure to the residuals of this degree 2 fit (and the residuals are still high). Let's try a degree 3 polynomial fit (known as cubic):$f(x)=ax^3+bx^2+cx+d$<jupyter_code>fit = np.polyfit(Betoule_data['distance'],Betoule_data['velocity'],3) 
y_values = np.polyval(fit,Betoule_data['distance'])

plt.subplot(2,1,1)
plt.scatter(Betoule_data['distance'],Betoule_data['velocity'])
plt.plot(Betoule_data['distance'],y_values,color='orange',)
plt.title('data and a polynomial degree 3 fit')
plt.ylabel('Velocity (km s$^{-1}$)')
plt.xlabel('Distance (Mpc)')

plt.subplot(2,1,2)
plt.scatter(Betoule_data['distance'],Betoule_data['velocity']-y_values)
plt.title('residuals of a polynomial degree 3 fit')
plt.ylabel('Residual velocity (km s$^{-1}$)')
plt.xlabel('Distance (Mpc)')

plt.tight_layout()
plt.show()<jupyter_output><empty_output><jupyter_text>Can a degree 4 polynomial fit do better?$f(x)=ax^4+bx^3+cx^2+dx+e$<jupyter_code>fit = np.polyfit(Betoule_data['distance'],Betoule_data['velocity'],4) 
y_values = np.polyval(fit,Betoule_data['distance'])

plt.subplot(2,1,1)
plt.scatter(Betoule_data['distance'],Betoule_data['velocity'])
plt.plot(Betoule_data['distance'],y_values,color='orange',)
plt.title('data and a polynomial degree 4 fit')
plt.ylabel('Velocity (km s$^{-1}$)')
plt.xlabel('Distance (Mpc)')

plt.subplot(2,1,2)
plt.scatter(Betoule_data['distance'],Betoule_data['velocity']-y_values)
plt.title('residuals of a polynomial degree 4 fit')
plt.ylabel('Residual velocity (km s$^{-1}$)')
plt.xlabel('Distance (Mpc)')

plt.tight_layout()
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10934065934065934, 'fraction_numerical': 0.018956043956043955, 'mean_word_length': 5.600483529767302, 'pattern_counts': {'":': 0, '<': 102, '<?xml version=': 0, '>': 106, 'https://': 5, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '25', 'dir': 'jupyter-structured-clean-dedup', 'id': '145630', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 7264, 'n_tokens_neox': 6652, 'n_words': 2602, 'path': 'week_09/Week09_Age_of_Universe_InClass.ipynb', 'repo_name': 'ds-connectors/EPS-88-FA20', 'size': '21844'}