Fitting Models to Data (Part 2)

Fitting Models to Data (Part 2)#

Prof. Ryan Trainor, Franklin & Marshall College

Objectives#

To introduce the role of measurement uncertainties in model evaluation
To compare logarithmic and linear scaling for visualizing data
To estimate the uncertainty on inferred model parameters

Learning Outcomes#

By the end of this workbook, students should be able to

Include measurement uncertainties in the presentation of data and model fitting
Compute the Chi-square value for a model and use it to evaluate the goodness of fit
Use logarithmic scaling to recognize and evaluate powerlaw relationships
Use comments to document their code and functions

Prerequisites#

This notebook assumes that the student has already completed Part 1 of this Model Fitting activity and is familiar with its prerequisites

Dependencies#

This notebook requires the standard libraries for scientific Python (numpy, SciPy) as well as the Astropy library. Astropy is included in the full installation of the Anaconda Distribution of Python, or it can be installed on its own by running conda install astropy or pip install astropy.

Note for instructors: this notebook has been written as an activity to be completed by a group of students in a lab or classroom setting. A solution with filled-in answers is available upon email request.

Review: Plotting the data#

Our first step is to create the same plot we made in the final step of Part 1 of the Model Fitting activity: i.e., to make a plot of luminosity vs. mass that includes the best-fit power-law function. You do not need to include the linear or quadratic curves we considered in that activity in this plot.

NOTES:

Feel free to copy and paste code from Part 1, but try to leave out any unnecessary commands: this is a good exercise to see what is the minimal amount of code you need to make this plot!
You will need the data file from Part 1 if you do not have it. You can download the data at this link or using the wget command in the cell below.
As you write or copy your code, practice using comments. Remember that comments are sections of code that are not read by python—they only exist to help you (or other people) understand what your code is doing. A comment starts with a # sign, and everything after the # on the same line is ignored by python. For each line or section of 2-3 lines that you copy from your Fitting Models to Data (Part 1) code, add a comment to explain what your code is doing in that section and why you need it.
While you’re looking at your previous lab, make sure you refresh your memory about the Residual Sum of Squares (RSS) and how it relates to the best-fit line. You’ll need this information later for Checkpoint 1.

# Download the file 'star-mass-lum-lumerr.dat' using wget (use this if you are running on Google Colab)
!wget --no-check-certificate 'https://drive.google.com/uc?export=download&id=1BWet8_t2gyNGtFNMP_1opqcYBHDinEfG' -O star-mass-lum-lumerr.dat

--2024-04-25 20:25:14--  https://drive.google.com/uc?export=download&id=1BWet8_t2gyNGtFNMP_1opqcYBHDinEfG
Resolving drive.google.com (drive.google.com)... 142.250.190.110, 2607:f8b0:4009:801::200e
Connecting to drive.google.com (drive.google.com)|142.250.190.110|:443... connected.
HTTP request sent, awaiting response... 

303 See Other
Location: https://drive.usercontent.google.com/download?id=1BWet8_t2gyNGtFNMP_1opqcYBHDinEfG&export=download [following]
--2024-04-25 20:25:14--  https://drive.usercontent.google.com/download?id=1BWet8_t2gyNGtFNMP_1opqcYBHDinEfG&export=download
Resolving drive.usercontent.google.com (drive.usercontent.google.com)... 

172.217.5.1, 2607:f8b0:4009:809::2001
Connecting to drive.usercontent.google.com (drive.usercontent.google.com)|172.217.5.1|:443... connected.

HTTP request sent, awaiting response...

200 OK
Length: 7790 (7.6K) [application/octet-stream]
Saving to: ‘star-mass-lum-lumerr.dat’


star-mass-lum-lumer   0%[                    ]       0  --.-KB/s               
star-mass-lum-lumer 100%[===================>]   7.61K  --.-KB/s    in 0s      

2024-04-25 20:25:14 (98.2 MB/s) - ‘star-mass-lum-lumerr.dat’ saved [7790/7790]

# Re-create the final plot from Part 1, omitting unecessary steps

Displaying error bars#

Once you have reproduced your plot from Part 1, look at the data compared to the best-fit line. You should see that the line does not perfectly go through all of the points. Does this mean our model is a bad model? Maybe, but maybe not—if there are measurement errors in our data, then it may be that our measurement errors are causing the data and model to disagree.

Thankfully, our data include an estimate of our uncertainty on each of the luminosity measurements (which we have ignored in Part 1). These uncertainties are given in the 4th column of the Astropy Table that includes our data. This column is named ‘Lum_err_Lsun’, so we can print the contents of this column by using the name of the Table indexed by the column name. Test this below:

# Print the column named 'Lum_err_Lsun' in your Astropy Table
print(?)

  Cell In[3], line 2
    print(?)
          ^
SyntaxError: invalid syntax

Given these uncertainties, we can now add error bars to the plot we created previously. You can do this by replacing the plt.plot() command that you used to plot the measured masses and luminosities with the plt.errorbar() command. This command has the following syntax:

plt.errorbar(x,y,yerr,marker=?,ls=?,color=?)

Here, x is your array of masses, y is your array of luminosities, and yerr is your array of luminosity uncertainties. marker allows you to specify a symbol to mark the points, ls allows you to set a linestyle that connects the points, and color allows you to specify the color. Let’s use orange squares, and let’s not use a line to connect them. Remember that you can find information on format options by googling “matplotlib plot formats” or going to the following link: https://matplotlib.org/3.1.1/api/_as_gen/matplotlib.pyplot.plot.html. You may also find ls='none' useful).

Also, continue adding comments to explain what you’re doing and changing with each step

# Create your plot with errorbars here:

Error bars are useful for determining whether the model is statistically consistent with the data. Note that we cannot know the actual error on each measurement: if we knew the the error on a given measurement, we would simply subtract the error to get the true value!

Instead, our uncertainty (as shown by our error bars) represents the magnitude of the typical error we expect for a single measurement: i.e., it represents the the expected magnitude of the residual (the difference between the data and the model, introduced in Part 1) for a true model. For a good model, we expect that 68% of the points will have error bars that include the model line, and 95% of the points will be less than 2x the error bar away from the model line.

CHECKPOINT 1#

If you are completing this activity in class, stop here and check your plot with your instructor or TA before you continue. You may want to think about the following questions in advance:

What does the “residual sum of squares” (RSS) mean, and how does it relate to your best-fit line?
Which stars in your plot have the the most precisely-determined luminosities? Which have the most uncertain luminosities?
Does your best-fit line seem consistent with your data, given the uncertainties in your fit?

Using a logarithmic plot#

One thing you might notice about your plot above is that most of the points are all clustered together at very small values in the bottom-left corner. This is a common issue in astronomy, where we often deal with numbers that cover many orders of magnitude in scale—if you make your plot big enough to accommodate the big numbers, then you can barely see the differences among the small numbers!

One way to deal with this is to use a plot with logarithmic scaling: this will change our x and y axes to show powers of 10 instead of linear values. Python has an easy way to do this by adding a single command to the script for our plot.

Sometimes we have data that is evenly distributed on either the x or y axes but heavily skewed on the other axis. In these situations, you can use plt.semilogx() to make the x axis logarithmic while leaving the y axis on a linear scale, while plt.semilogy() does the opposite. If both axes should be shown logarithmically, we can use plt.loglog(). Re-create your plot from above (with the errorbars) and try adding one of these commands. Test all three commands and decide as a group which one makes the data easiest to see.

# Create your plot with errorbars and logarithmic scaling here:

If you have not already included your best-fit powerlaw relationship, make sure to display it as well. What shape does it make? Hopefully you find that it looks like a straight line (i.e., the power law looks linear when plotted logarithmically).

One other thing you may find is that our power-law model that looked so good in linear space suddenly looks like a poor fit, especially for the low masses: it systematically misses all of the low-mass points! What’s going on here?

Remember that we made our fit to minimize the RSS: we wanted a line that was “as close as possible” to all of the points, and each point was treated equally. This means that being off by 1 L\(_\odot\) from our largest luminosity was just as bad as being off by 1 L\(_\odot\) from our smallest luminosity. However, that doesn’t really make sense—1 L\(_\odot\) is a tiny fraction of the total luminosity of our bright stars, but it’s a huge fraction for our faint stars. Also, and more importantly, remember what we saw for our error bars: the largest luminosities are the most uncertain, whereas the small values are more precisely known in our data. With this in mind, we should focus on trying to get a better fit to the points we know precisely, even if it means we get a slightly worse fit to the other points. That is, we need to give the points with small uncertainties additional weight, so that they count more in the fit.

In order to create a fit that weights points according to their uncertainties, we will need to discuss the \(\chi^2\) value (pronounced “\(\chi\)-squared”, where \(\chi\) is a Greek letter that physicists pronounce like “KAI”).

\(\chi^2\) values#

The \(\chi^2\) value of a fit is very similar to the RSS. The only difference is that each residual is weighted by the inverse of its uncertainty, so that points with small uncertainties are given larger weight. It is defined as follows:

\[ \chi^2=\sum_{i=1}^N \left(\frac{\Delta y_i}{\sigma_i}\right)^2=\sum_{i=1}^N \left(\frac{y_{i,\mathrm{obs}}-y_{i,\mathrm{pred}}}{\sigma_i}\right)^2 \]

As in Part 1, \(\Delta y_i\) is the residual of point i, which is the difference between the observed value at point i (\(y_{i,\mathrm{obs}}\)) and the value at point i predicted by the model (\(y_{i,\mathrm{pred}}=f(x_i)\), where \(f(x)\) represents our model). \(\sigma_i\) is the uncertainty in the y value for that particular data point. Note that if don’t divide by \(\sigma_i\) (or if \(\sigma=1\) for every point) then the \(\chi^2\) value is equivalent to the RSS.

In the box below, do the following three things:

Create a function that will calculate the \(\chi^2\) value for arbitrary arrays yobs, ypred, and sigma.
Also create a second function that computes our old friend the RSS based on yobs and ypred.
Finally, use your function to find the \(\chi^2\) value for the best-fit power-law model you found above.

# Fill in the functions below

def chi2(?,?,?):
    return ?
    
def RSS(?,?):
    return ?

chi2_orig = ?
print(chi2_orig)

If your calculation is working, you find a \(\chi^2\) value of approximately 5000. This is a lot! Remember: if our uncertainties and our model are both accurate, then on average our residuals will be equal in magnitude to our uncertainties: \((\Delta y)^2=\sigma^2\), and thus for a true model:

\[ \chi^2=\sum_{i=1}^N \left(\frac{\Delta y_i}{\sigma_i}\right)^2\approx\sum_{i=1}^N \left(1\right)^2\approx N \]

Since our table has 235 stars in it (you can check this with print(len(data)) if you named your Table data), then we should have \(\chi^2\approx 235\). The fact that the model above has a much larger \(\chi^2\) value means that our residuals are much larger than our uncertainties on average, and thus that there is a significant disagreement between our data and the model that cannot be explained by measurement errors alone.

\(\chi^2\) minimization#

Now our goal is to create a new fit that takes our uncertainties into account. We can do this by passing an optional argument called sigma to curve_fit. In the cell below use curve_fit() as you have previously, but this time set sigma equal to the array of luminosity uncertainties:

popt,pcov=?
A_new,alpha_new=popt
print(A_new,alpha_new)

In the cell below, now create a plot showing the data (with errorbars) and two power-law curves: one with your original best-fit parameters, and one with your new best-fit parameters. You may want to look at your plot in the log-log view and the linear view as you consider the questions below.

# Plot the luminosity vs. mass for the real data with errorbars
# include the model lines from the original fit and the fit weighted by the uncertainties

CHECKPOINT 2#

If you are completing this activity in class, stop here and check your plot with your instructor or TA before you continue. You may want to think about the following questions in advance:

Which line seems to be consistent with more of your points (i.e., within the uncertainties)?
Which line is more precise at predicting the value of the star with \(M=40\) M\(_\odot\)? Is this the best fit overall?
Which line has a smaller \(\chi^2\) value? Which has a smaller sum of squared residuals? You can use the cell below to do some calculations to check. Is either \(\chi^2\) value consistent with a “good” fit?

# Check chi^2 and RSS of original model and new model

Uncertainties in the model#

Because our data has some uncertainty on it, our model parameters have some uncertainty as well. What does this mean? Let’s imagine that we made new measurements of all of the luminosities of these stars. We might expect to find slightly different luminosities the second time, with the differences between the first and second measurements (i.e., the difference in measurement error) being similar in size to our uncertainties.

If we then fit a best-fit line to our new measurements, we would expect that we might get a slightly different “best” set of parameters A and \(\alpha\) as well. How different might we expect them to be? This question is asking about the uncertainties on our best-fit parameter values. If we can calculate uncertainties on the parameters, then we can be reasonably sure that anyone else calculating a similar fit would get parameters that are consistent with our to within ~2x the uncertainty (technically, the values should be within 2x uncertainties 95% of the time).

The uncertainties on our best-fit parameters are given by the covariance matrix, which is the pcov variable we’ve been ignoring so far. A covariance matrix is an \(N\times N\) matrix, where \(N\) is the number of parameters in our model. In our case, our parameters are A and alpha, so \(N=2\) and the covariance matrix looks like this:

\[\begin{split} \begin{bmatrix} \sigma_{AA} & \sigma_{A\alpha} \\ \sigma_{\alpha A} & \sigma_{\alpha\alpha} \end{bmatrix} \end{split}\]

The diagonal elements of this matrix (\(\sigma_{AA}\) and \(\sigma_{\alpha\alpha}\)) represent the variance (i.e., the square of the uncertainty) on our estimated parameters, which is to say that the uncertainties on our parameters estimated by our \(\chi^2\) minimization are the square-roots of pcov[0,0] and pcov[1,1] for A and alpha, respectively.

Given the above, complete the code block below to print the uncertainties:

A_err = ?
alpha_err = ?

print("The best-fit parameters are A={:4.2f}+/-{:4.2f} and alpha={:4.2f}+/-{:4.2f}".format(A_new,
                                                                                            A_err,
                                                                                            alpha_new,
                                                                                            alpha_err))

One way to visualize the uncertainties on A and alpha is to plot many different variations of our power law model where we slightly change the values of A and alpha in a manner consistent with the uncertainties on those parameters. The function below generates 100 random pairs (A,alpha) that are consistent with the uncertainties on those parameters:

def rand_A_alpha(A,alpha,A_err,alpha_err):
    '''Given input values of A and alpha and uncertainties on each value, generate
    100 pairs (A,alpha) such that the generated parameters form a gaussian
    distribution about the input values with a 1-sigma width given by the 
    uncertainty.'''

    n_rand = 100
    rand_A_array = A + A_err * np.random.normal(size=n_rand)
    rand_alpha_array = alpha + alpha_err * np.random.normal(size=n_rand)

    return zip(rand_A_array,rand_alpha_array)

In the cell below, create a new version of your plot with the logarithmic scaling, and use the rand_A_alpha() function above to generate 100 model lines that each use randomly-generated values of A and alpha.