Error Analysis Notes

From Physics 111-Lab Wiki

Jump to: navigation, search

I. Error Analysis Exercise and Notes II. MatLab Scripts for curve fitting. If inside the Physics111 Library site GoTo Advanced Lab Reprints then click on MatLab folder.

Matlab is the accepted program to use for Physics 111-Lab. Use it in the 111-Lab or see Don for remote access from off campus.


Contents

Introduction

View the Video a must see video introduction to error analysis.

First run through one of the several good texts on error analysis, L. Lyons, "A Practical Guide to Data Analysis for Physical Science Students", Cambridge University Press (1994). There is no substitute for this reading. To give all details here would require us to write a whole book. We are only going to skim the subject of statistical errors.

Error of a Single Parameter

Measurement statistics

As we read about laboratory results, we often see the value of a measured parameter as written something like 2.10 ± 0.05. The value 2.10 in this example is usually an average of several measurements, because successive measurements of any particular quantity are never exactly the same. Where does ± 0.05 come from, and what does it mean?

When we make a measurement and get a value for a single quantity or parameter, we want to know how likely it is to be correct. So we measure it again and again, say N times in all, and look at the scatter in the measurements. For example, we might have y = 2.10 as the average of the 5 measurements yi = 2.18, 2.03, 1.95, 2.24, 2.10,

 \bar {y} = \frac {\Sigma y_i}{N} = \left [2.18 + 2.03 + 1.95 + 2.24 + 2.10 \right ]/5 = 2.10

We then calculate the error of the mean, with the equation

\sigma_{\bar{y}} = \frac {Standard~Deviation, \sigma}{\sqrt{N}} = \frac {\sqrt { \left [ \Sigma \left ( \bar {y} - y_i \right )^2 \right ] }} { \sqrt { N \left ( N - 1 \right ) }} = 0.05


Loosely speaking, this means that if we make another set of N measurements, the average of this new set has a 67% likelihood of falling within a range from 2.05 to 2.15, or more compactly written as 2.10 ± 0.05

Counting statistics

In experiments like Muon Lifetime, we record the number of events or counts in a given time interval, where the events are random, as in radioactive decay. Suppose we get N counts in a fixed time interval. If we take other samples for the same interval, we get some close but different values of N. After many such measurements, we can construct the standard deviation, as described above. We don't need to do this, because statistical considerations specify that if the counts come at truly random times, then the standard deviation is \sqrt {N}. So, we can make a single record of counts and write N ± \sqrt {N} for the counts in a given interval of time. The larger the N, the larger the standard deviation, but the fractional error grows smaller as the number of counts grows larger. The fractional error is \frac {1}{\sqrt {N}}. For example, to get a fractional error of 1%, N must be 10 000!

Curve Fitting with Two Parameters (Linear Regression)

In many laboratory experiments we have several parameters, rather than just one, that we want to determine. One approach is to fit a smooth analytical curve, a function of the parameters, to the datum points and extract the most likely values of the parameters from the fit. In cases where the curve is a straight line, we want to fit a data set to determine the slope and intercept of the line, each of which has some physical significance.

We start with the equation y = mx + b = f(x). We assume that the independent variable is x and we record a datum point yi at each value xi. We then select a straight line and adjust its slope "m" and intercept "b" to get a "best fit" to the measured yi with an error σi. By definition we get a best fit when the sum of the squares of the deviations of the fitted curve from the datum points at each xi is a minimum. In mathematical terms, adjust m and b so as to get a minimum of the sum:

\chi^2 = \Sigma \frac{1}{\sigma_i^2} \left [ y_i - f \left ( x_i \right ) \right ]^2 = \Sigma \frac{1}{\sigma_i^2} \left [ y_i - \left ( mx_i + b \right ) \right ]^2

We do this by taking the derivatives of the sum, χ2, with respect to m and b, setting the equations equal to zero, and solving the simultaneous equations for the "best" values of m and b. This is called a "least-squares fit". These values will have errors that can be calculated with techniques described in the texts. Lyons gives a worked-out example.

References

  1. L. Lyons, "A Practical Guide to Data Analysis for Physical Science Students", Cambridge University Press (1994). Louis Lyons.
  2. Yardley Beers, "Introduction to the Theory of Error", Addison Wesley Publishers, (1957), Yardley Beers This is a short book which is very good.
  3. P. Bevington, '"Data Reduction and Error Analysis for the Physical Sciences", McGraw-Hill. [An old standard that is pretty dry but straightforward. Chapter 5 is particularly important.] Computer Programs available:McGrawHill
  4. A. C. Melissinos and J. Napolitano, "Experiments in Modern Physics, 2nd Edition", Academic Press (2003).
  5. W. H. Press, et al., "Numerical Recipes in C: The Art of Scientific Computing, 2nd Edition", Cambridge University Press (1992); refer to Ch. 14—"Modeling of Data". [The Numerical Recipes in Pascal or FORTRAN books contain identical information. This book is the standard reference for doing scientific work on computers. Chapter 14 has a good introduction to the method of maximum likelihood, chi–square fitting, modeling data in general, error estimates of fit parameters, and, important for later experiments, the Monte Carlo method of simulation.]

We have available in the Physics 111-Lab Books 1, 2, 3, and 4; Number 1 and 2 are on-line above.

Problem Set

Problem 1

We want to measure the specific activity (number of decays per second) of a radioactive source so that we can use it to calibrate the equipment of the gamma-ray experiment. We use an electronic counter and a timer to measure the number of decays in a given time interval. In round numbers we measure 1000 decays in 5 minutes of observation. For how long should we measure in order to obtain a specific activity with a precision of 1%? Explain.

Problem 2

You are given two measurements of distance and their associated uncertainties: A \pm \sigma_A and B \pm \sigma_B . Calculate the propagated uncertainty in the (a) total distance A + B, (b) difference A - B, (c) the perimeter 2A + 2B and (d) the area  A \times B .

Problem 3

In this problem we will be generating and analyzing lists of normally distributed random numbers. The distribution we are sampling has true mean 0 and standard deviation 1.

  1. If we sample this distribution N=5 times, what do we expect the mean to be? How about the standard deviation? Whats the error on the mean?
  2. Using Matlab, generate a list of N=5 normally distributed random numbers (the command randn(N,M) will generate M lists of length N). Calculate the mean, standard deviation and the error on the mean. Is this what you expected?
  3. Now find the mean and standard deviation for M=1000 experiments of N=5 measurements each. Plot a histogram of the distribution of the means. About how many experiments are within 1 sigma of the true mean of 0. About how many are within 2 sigma? Is this what you expected?
  4. Now repeat questions 1-3 for N=10,50,100,1000.

Problem 4

In this problem we will repeat the above process, but now using lists of exponentially distributed random numbers. The probability of selecting a random number between x and x+dx is \propto e^{-x}dx.

  1. What do you expect to be the mean of the distribution? What do you expect to be the standard deviation? (Note: The standard deviation is defined exactly as it is for a normal distribution, but the 1 sigma = 68% rule no longer applies to a exponential distribution.) What do you expect to be the error on the mean for N=100? Given M=1000 lists of N=100 random numbers, what do you expect the distribution of means to look like? What is the error on the mean?
  2. Make a list of N=100 exponentially distributed random numbers, (Hint: this can be easily done starting with a uniform distribution of random numbers.) Calculate the mean and standard deviation.
  3. Make M=1000 lists of N=100 exponentially distributed random numbers. Make a histogram of the means. Does the distribution of means look as you thought? What is the standard deviation of the means. Does this agree with what you thought?
  4. Repeat the previous steps for N=1000 & 10000. Does the error on the mean scale as you thought?

This is a demonstration of the Central Limit Theorem.

Problem 5

You are given a dataset (Media:Peak.zip) from a gamma-ray experiment consisting of ~1000 hits. For each hit, the energy of the gamma-ray is recorded. We will assume that the energies are randomly distributed about a common mean, and that each hit is uncorrelated to others. Read the dataset from the enclosed file and:

  1. Produce a histogram of the distribution of energies. Choose the number of bins wisely, i.e. so that the width of each bin is smaller than the width of the peak, and at the same time so that the number of entries in the most populated bin is relatively large. Since this plot represents randomly-collected data, plotting error bars would be appropriate (hint: use errorbar function in Matlab).
  2. Compute the mean and standard deviation of the distribution of energies and their statistical uncertainties.
  3. Fit the distribution to a Gaussian function, and compare the parameters of the fitted Gaussian with the mean and standard deviation computed above.
  4. How consistent is the distribution with a Gaussian? In other words, compare the histogram from (1) to the fitted curve, and compute a goodness-of-fit value, such as χ2 / df.

Problem 6

In the optical pumping experiment we measure the resonant frequency as a function of the applied current (local magnetic field). Consider a mock data set:

Current I (Amps) 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2
Frequency f (MHz) 0.14 0.60 1.21 1.94 2.47 3.07 3.83 4.16 4.68 5.60 6.31 6.78
  1. Plot a graph of the pairs of values. Assuming a linear relationship between I and f, determine the slope and the intercept of the best-fit line using the least-squares method with equal weights, and draw the best-fit line through the data points in the graph.
  2. From what s/he knows about the equipment used to measure the resonant frequency, your lab partner hastily estimates the uncertainty in the measurement of f to be σf = 0.01 MHz. Estimate the probability that the straight line you found is an adequate description of the observed data if it is distributed with the uncertainty guessed by your lab partner. (Use Matlab to calculate it or look it up in a table. For example, see table C-4 in Bevington). What can you conclude from these results? Repeat the analysis assuming your partner estimated the uncertainty to be σf = 1 MHz. What can you conclude from these results?
  3. Assume that the best-fit line found in the previous exercise is a good fit to the data. Estimate the uncertainty in measurement of y from the scatter of the observed data about this line. Again, assume that all the data points have equal weight. Use this to estimate the uncertainty in both the slope and the intercept of the best-fit line. This is the technique you will use in the Optical Pumping lab to determine the uncertainties in the fit parameters.
  4. Now assume that the uncertainty in each value of f grows with f: σf = 0.03 + 0.03 * f (MHz). Determine the slope and the intercept of the best-fit line using the least-squares method with unequal weights (weighted least-squares fit).

Problem 7

  1. In the muon lifetime experiment we obtain a histogram for the decay rate as a function of the time after the muon enters the detector and announces its presence. We expect the distribution (the histogram) to be described by an exponential function. Rather than fitting with an exponential function, it is more convenient to plot the logarithm of the decay rate as a function of time and then fit a straight line to it. Since each data point (xi,yi) has a statistical error, σi, associated with it, qualitatively, what happens to these errors when the semi-log histogram (xi,logyi) is plotted? Explain and illustrate. what happens, quantitatively? Assume yi is reasonably large.
  2. In a separate experiment, you calculate that  \log E_0 = 2.1 \pm .5 . What is the E0 and its experimental bounds? (Note that .5 is not small compared to 2.1).
Personal tools