W2 calculates wrong errors

[sschmidt]

W2 computes the wrong standard deviations:

Instead of (examples for unweighted values)

\sqrt{ \frac{1}{n-1} \sum ( x_i - <x> )^2 }

it computes

\sqrt{ \frac{1}{n-1} \sum ( x_i )^2 }

See analysis.f90, function observable_get_error (obs)

   err = sqrt (obs%sum_squared_values / obs%sum_squared_weights) &
          / sqrt (real (obs%count - 1, default))

[sschmidt]

example file with two histograms, these histograms are mirrors of each other, they should hav the same error.

[sschmidt]

With r2042 I committed a tentative fix, the error is now not calculated using

\sqrt{ \frac{1}{n(n-1)} \sum_i ( x_i - <x> )^2 }

which would have the need to store all values of x_i for a correct calculation, but

\sqrt{ \frac{1}{n(n-1)} \sum_i ( x_i - <x>_i )^2 } with <x>_i the average of the x from x_1 to x_i

This should be alright as long as

• n is not small,
  
• the observable is not correlated with the event number.
  

Both should be sufficiently fulfilled.
I'll leave this open in case someone wants to implement a correct solution.

[ohl]

Replying to sschmidt:

\sqrt{ \frac{1}{n(n-1)} \sum_i ( x_i - <x>_i )^2 } with <x>_i the average of the x from x_1 to x_i

It's wasteful, because you have to calculate the average every time observable_record_value_* is called. But probably this doesn't matter much in the grand scheme of things.

More importantly, I'm not sure that this is an unbiased estimator of the error. Why not

\sqrt{ \frac{1}{n(n-1)} \sum_i x_i^2 - \frac{1}{n^2(n-1)} (\sum x_i )^2  }

which is both faster and unbiased (please check that I didn't screw up the powers of n ...)

[sschmidt]

Replying to ohl:

More importantly, I'm not sure that this is an unbiased estimator of the error. Why not

\sqrt{ \frac{1}{n(n-1)} \sum_i x_i^2 - \frac{1}{n^2(n-1)} (\sum x_i )^2  }

which is both faster and unbiased (please check that I didn't screw up the powers of n ...)

Yes that's correct and better (didn't think of this formula at first.) Even the powers of n are correct.

[kilian]

I implemented the formula in r2057. Please double-check, then close the ticket.

[ohl]

Looks good to me ...

[Juergen Reuter]

Milestone v2.0-rc4 deleted