Re:Re:Hi, I have answered below: On 29 November 2010 13:35, Dr. Klaus-Peter Neidig <peter.neidig@xxxxxxxxxxxxxxxxx> wrote:
Hi Edward, after reading the mails again and checking the code again I think it would be good to do another cross-check, perhaps best with an artificial example: suppose we have mixing times 10, 10, 10 (3 replicates) 20 30, 30 (2 replicates) 40 50 and spectra with 3 peaks. A: I would take peak1, then check for the largest difference in intensity in the first 3 replicates and in the 2 replicates. The largest found difference is then assigned as a systematic error to peak1 at all mixing times.
Is an absolute value is used for finding the maximum difference? So say in the mixing time of 10, the intensities are 1, 4, 2 so the max difference is 3. And in the 30 has intensities of 0.9 and 3. So the error of peak 1 would be 3. Is this correct?
The procedure is repeated for peak2 and peak3. Most likely peak1, peak2 and peak3 then have different systematic errors. That would be like a worst case scenario per peak and is what currently was active in the PDC.
Oh, this is different to what I though was happening. I assumed that this value was averaged, or the maximum of peaks 1, 2, and 3 was taken as the error for all peaks. This would probably result in very noisy error estimates as the amount of error sampling is only 2 (in this case). This will also overestimate the error. I can set up a test to see what happens. I am using a Gaussian centered at 20 with a standard deviation of 1. I take 2 points (peak intensities from duplicated spectra) randomly from the distribution and calculate the difference. This is then repeated N times (number of duplicated spectra) and the maximum difference is taken. I repeat this M times for determining some stats, where M = 1e6. The average error estimate is: N = 2, ave(error) = 1.594 N = 3, ave(error) = 1.875 N = 4, ave(error) = 2.070 N = 5, ave(error) = 2.221 N = 6, ave(error) = 2.339 N = 100, ave(error) = 3.884 N = 2, sd(error) = 0.852 N = 3, sd(error) = 0.829 N = 4, sd(error) = 0.806 N = 5, sd(error) = 0.786 N = 6, sd(error) = 0.770 N = 100, sd(error) = 0.566 These are all log-normal distributions. As you can see, the error estimate is always on average overestimated, and the more duplicated spectra N, the worse this becomes. The spread of values is also a worry. With 3 duplicated spectra, the resultant errors are 1.875 +/- 0.829 which means that the error estimates are all over the place. The perfect estimate would be 1.000 +/- 0.000, as the sd is exactly 1. The Python script, an xmgrace, and a PS plot for N=3 can be found at https://gna.org/task/index.php#comment1. Does this test match the technique that the PDC uses? Of course, there could be a bug in my script (https://gna.org/support/download.php?file_id=11436)!
B: The absolute worst scenario would be to check for the largest ever ocuring difference and assign that to all 3 peaks at all mixing times. This we used in earlier versions of the PDC.
In this case there is no sd(error). So using 100 peaks and repeating the calculation a number of times: [edau@localhost peak_intensity_errors]$ python sampling_test2.py N = 2, error = 5.005 N = 3, error = 4.485 N = 4, error = 3.863 N = 5, error = 4.268 N = 6, error = 5.392 N = 100, error = 5.353 [edau@localhost peak_intensity_errors]$ python sampling_test2.py N = 2, error = 4.161 N = 3, error = 4.056 N = 4, error = 5.135 N = 5, error = 3.890 N = 6, error = 4.412 N = 100, error = 5.469 [edau@localhost peak_intensity_errors]$ python sampling_test2.py N = 2, error = 3.550 N = 3, error = 4.184 N = 4, error = 4.574 N = 5, error = 4.182 N = 6, error = 5.080 N = 100, error = 5.973 [edau@localhost peak_intensity_errors]$ python sampling_test2.py N = 2, error = 4.549 N = 3, error = 3.824 N = 4, error = 5.337 N = 5, error = 4.507 N = 6, error = 4.113 N = 100, error = 6.410 [edau@localhost peak_intensity_errors]$ python sampling_test2.py N = 2, error = 3.966 N = 3, error = 4.667 N = 4, error = 4.770 N = 5, error = 4.708 N = 6, error = 5.058 N = 100, error = 5.607 The Python script for this is attached at https://gna.org/task/index.php?7180#comment2. This technique always overestimates the value of 1 and is hyper-variable, but increasing with increasing number of duplicated time points.
C: Your variance averaging would work as follows: loop over the first 3 replicates, calculate variance for each peak loop over the 2 replicates, calculate variance for each peak finally take the average over 6 variances and assign this as a systematic error to all peaks at all mixing times Is this correct ?
This is correct. Using 100 peaks, duplicates for 1 time point (ignoring the multiple relaxation periods for now), and M = 1e6 for some stats on the error estimate, the results are: ave(error) = 0.998 sd(error) = 0.071 This is for Gaussians of mu = 20 and sigma = 1. As you can see, the standard deviation of 1 is nicely estimated (plus and minus 0.071 for different attempts). This is the best error estimate. The Python script is at https://gna.org/task/index.php#comment3. The PS plot is at https://gna.org/task/download.php?file_id=11442. If we now extend this to multiple relaxation periods, that would be equivalent to M ~= 10 but with variance averaging rather than sd averaging. The error estimate would not be as good as 0.998 but the more replicated spectra, the closer the estimate will be to the real value.
I think A: and C: should be offered. After implementing this I see my data that the C: produces much larger errors than A and the errors of the fitted parameters get much larger, sometimes larger than the values itself. This seems not acceptable.
From my tests, the errors should be bigger with the A procedure. This
is strange! It might be worth seeing if the C procedure implementation matches that in https://gna.org/support/download.php?file_id=11440. I may have also have made a mistake within the logic of these testing scripts!
In the mean I also spoke to Wolfgang Bermel. My statement that peaks with longer T2 seem to very more in replicates is not obvious to him. When looking to the data it turned out that I had looked at peaks that do not only have longer T2 but also have larger resonance offsets. Wolfgang suspects that their might be an offset effect. The statement that spectra taken at longer mixing times show a lower base plane noise seems to be not true in general. The trend I have seen in one of the spectra is not fully convincing so we checked a couple of other spectra all taken with the latest release pulse programs. With those data we would not really confirm the trend. We concluded to have an eye on that from now on.
Ah, I remember also talking to him at ICMRBS 2008 about this and he said the same. The single scan interleaving might have fixed this. Maybe it was an old sample heating problem in the R2 that I remembered. Or maybe the problem was that I measured on a Varian machine ;) But thanks for asking. I am continually asking if people see this. Some do, and some don't, but it's hard to tell what is happening. Cheers, Edward
Attachment:
sampling_test.ps
Description: PostScript document