Model measurement and error in NumPy
I'd like to try the SciPy suite instead of Octave for doing the statistics in my lab experiments. Most of my questions were answered here, there is just another thing left:
I usually have an error attached to the measurements, in Octave I just did the following:
R.val = 10; R.err = 0.1; U.val = 4; U.err = 0.1;
And then I would calculate I with it like so:
I.val = U.val / R.val; I.err = sqrt( (1 / R.val * U.err)^2 + (U.val / R.val^2 * R.err)^2 );
When I had a bunch of measurements, I usually used a structure array, like this:
R(0).val = 1; R(0).err = 0.1; … R(15).val = 100; R(15).err = 9;
Then I could do R(0).val or directly access all of them using R.val and I had a column vector with all the values, for mean(R.val) for instance.
How could I represent this using SciPy/NumPy/Python?
The easiest is indeed to use NumPy structured arrays, that give you the possibility to define arrays of homogeneous elements (a record) composed of other homogeneous elements (fields).
For example, you could define
R = np.empty(15, dtype=[('val',float),('err',float)])
and then fill the corresponding columns:
R['val'] = ... R['err'] = ...
Alternatively, you could define the array at once if you have your val and err in two lists:
R = np.array(zip(val_list, err_list), dtype=[('val',float),('err',float)])
In both cases, you can access individual elements by indices, like R (which would give you a specific object, a np.void, that still gives you the possibility to access the fields separately), or by slices R[1:-1]...
With your example, you could do:
I = np.empty_like(R) I['val'] = U['val'] / R['val'] I['err'] = np.sqrt((1 / R['val'] * U['err'])**2 + (U['val'] / R['val']**2 * R['err'])**2)
You could also use record array, which are basic structured array with the __getattr__ and __setattr__ methods overloaded in such way that you can access the fields as attributes (like in R.val) as well as indices (like the standard R['val']). Of course, as these basic methods are overloaded, record arrays are not as efficient as structured arrays.
This kind of error propagation is exactly what the uncertainties Python package does. It does so transparently:
from uncertainties import ufloat R = ufloat(10, 0.1) U = ufloat(4, 0.1) I = U/R print I
(Disclaimer: I'm the author of this package.)