Fall Semester 2006
This problem is an introduction to linear fitting.
File ~p6720/data/f16.dat gives airspeed "v" versus the force of air resistance "f" for a single-engine jet flying at 40,000 ft above sea level. (The data are very crudely derived from the performance specifications of an F-16.) The airspeed is known precisely while the fractional errors in the force are at the 2.5% level (sigma_f/f=0.025).
A model for the air resistance force as a function of airspeed is F = C v β
Use the Numerical Recipes routine fit.cpp to perform a linear least squares fit of this model to the flight data. Hints for doing this are given in the October 25th lab exercise . Obtain from your fit the following information:
# best-fit beta: XXX +/- XXX # best-fit C: XXX +/- XXX # chi-square: XXX # Are the data and model consistent? Y/NAdditional hints: You first need to linearize the model. You then may translate the fractional errors in the force to absolute errors in the linear fit. (You are just "propagating errors" here. Think of sigma_x of some variate x as being a small finite change in x. With this in mind, a helpful relation might then be that delta_x/x is approximately dx/x=dln(x). This holds very well when delta_x is small compared to x.)
Finally, as mentioned in lecture, you can get a good sense of whether a model fits some data by looking at chi^2
chi^2 = Sum_{i=1 to N) ((y_i - model(x_i))/sigma_i)^2
where N is the number of data points (x_i,y_i), and model(x) is the fitting function (in this case a line). If chi^2 is about equal to number of "degrees of freedom" nu=N-M (where M is the number of model parameters -- it is 2 here) then the fit is good. If chi^2 is much larger, worry about a bad fit. If chi^2 is much smaller, worry that you have overestimated your uncertainties.
Submit your C++ code f16.cpp. It should compile with
CC -o f16 f16.cppYour code should assume that all necessary additional NR files are in the working directory. (Do not link to the p6720/nr library.)
Write the following to standard output upon execution:
# best-fit beta: XXX +/- XXX # best-fit C: XXX +/- XXX # chi-square: XXX # Are the data and model consistent? Y/NYou need not generate any plots for this exercise (although it is often a good idea to check your fit graphically as well as numerically).
Here we deal with generalized linear least squares fitting.
A file called ~p6720/data/grtest.dat has fictitious data from a pair of black holes in orbit around each other. Einstein's theory of General Relativity (GR) predicts a fit to this data according to the model
y = a_0 + a_1 * x + a_2 * x^2 + a_3 * x^3
where x and y are given in the datafile (along with errors in y)
and the coefficients a_i depend on Newton's gravitational constant,
black hole masses, and orbital parameters (among other things, possibly).
Note: this a linear model -- it is linear in the coefficients, a_i. Each a_i is multiplied by a "basis function", fn(x). Here, fn_i(x) [multiplying the coefficient a_i] is just x^i. In principle, the fn_i could be any function of x.
There are other competing descriptions of gravity such as the Brans-Dicke theory (see Misner, Thorne & Wheeler, Gravitation). Often these are extensions of GR which incorporate one or more extra parameters; in the limit of these parameters "going to zero", such theories become identical to GR.
Suppose an alternative (non-GR) model "GR++" predicts that the orbital data for the black-hole pair are dscribed by
y = a_0 + a_1 * x + a_2 * x^2 + a_3 * x^3 + a_4 * x^4.
The problem here is to find the
best-fit parameters a_i for both (GR and non-GR) models
and provide the following:
best-fit value of a_0 from GR: XXX+/-XXX best-fit value of a_0 from GR++: XXX+/-XXX chi2 from GR: XXX chi2 from GR++: XXX Goodness-of-fit q from GR: XXX Goodness-of-fit q from GR++: XXX best-pick model: GR/GR++Thus, you are to give the best-fit value of the constant term in both models, chi2 values, Goodness-of-fit values (see NR's Sections 14.3 (follow discussion on the "chi-square probability function") and 15.1 (subheading: "Chi-Square Fitting"). On the basis of the Goodness of Fit choose one model or the other as your "best-pick", the one that you think best describes the data. Give a few sentences justifying your answer.
In coming up with the best-pick model, apply a sharp "Occam's razor": Cut out all but the simplest of models, unless otherwise warranted by the data.
The fit you are being asked to perform here is a generalized linear least squares (GLLS) procedure ("linear" because the models are linear in the fitting parameters, and "general" becuase there can be an arbitrary number of parameters).
To perform the fits with the GR and GR++ models, you may take advantage of available high-quality codes written in Fortran, C, or C++. Here, use the Numerical Recipes code ~p6720/nr/recipes_cpp/recipes/lfit.cpp. An example driver routine is given in p6720/nr/recipes_cpp/demo/src/xlfit.cpp. The algorithm behind lfit.cpp is described fully in Section 15.4 of Numerical Recipes.
You should write a code grtest.cc which calls the lfit routine, assuming that all necessary additional NR files are in the working directory. (Do not link to the p6720/nr library.) In order to determine the goodness of fit, you will need to explicitly call the function NR::gammq (fit.cpp does this call for us, lfit.cpp does not). Remember that the incomplete gamma function takes two arguments, the first is 0.5*(number of degrees of freedom), the second is 0.5*chi**2.
Have your code read in the data grtest.dat (its also fine if you just hardwire the numbers in), and write the following information to standard output:
best-fit value of a_0 from GR: XXX+/-XXX best-fit value of a_0 from GR++: XXX+/-XXX chi2 from GR: XXX chi2 from GR++: XXX Goodness-of-fit q from GR: XXX Goodness-of-fit q from GR++: XXX best-pick model: GR/GR++
Remember that you can fit anything you want with this type of code so long as the fitting function is linear in the coefficients a[i]. Otherwise, the basis functions fx[i] can be anything you want!
Submit your code grtest.cc.