## #125 Verification of DescStats.createMeanObservationForRange()

closed-fixed
9
2010-05-20
2010-02-02
David Benn
No

I would appreciate review of the code that yields mean curve datapoints, namely: DescStats.createMeanObservationForRange()

The AAVSO web-based light curve generator (LCG) generates significantly different (1 sigma) error bars that VStar. Why? Take, for example the last 200 days (from today) of Eps Aur in LCG and VStar.

## Discussion

• David Benn - 2010-02-02

Some questions:

o Does LCG use Standard Error of the Average calcs as per http://www.aavso.org/education/vsa/Chapter10.pdf?
o Does LCG use both V and Visual bands to calculate the mean curve datapoints and so std error of the average values?

• Aaron Price - 2010-02-03

The lcg is supposed to show stdev based on the data in the bin. The mean curve can be either V or visual, the user gets to select which one in the interface. It's entirely possible my code is wrong too.

• David Benn - 2010-02-03

It would be great to look at this in depth when time permits. Thanks Aaron.

• David Benn - 2010-02-20

Apart from asking the question: is the algorithm/approach correct, it would also be worth coming up with some simple test data that we can use to compare VStar against LCG.

• David Benn - 2010-03-24

Matt provided this via email:

Sara communicated with Matt who provided this information:
Here is what I would do, written in natural language:
1) Choose a bin size (e.g. 7 days)
2) Place the center of the first bin at a time (binsize/2) from the first data point (i.e. if the first data point is 2455200.0, put the bin center at 2455203.5).
3) Sum all of the magnitudes in the bin to find the average magnitude, AND sum all of the jds to find the average time.
4) Compute the variance of the magnitudes in each bin using
VAR(i) = Sum (mag(j)-avemag(i))^2
5) Divide the variance of each bin by the number of data points in that bin, minus 1 (N-1) (i.e. Var(i)/(N(i)-1) )
6) Take the square root of that value, to get the std. deviation:
sigma(i) = sqrt (Var(i)/(N(i)-1))
7) Repeat, shifting the bin center one bin width up in time each time (e.g. 2455203.5, 2455210.5, etc.)
The averaged light curve should then be JD(ave),mag(ave), and sigma for each bin.

• David Benn - 2010-03-24

I need to fix the 2nd part of step 3. Currently, oddly, I take the first and last JD in the bin and take the average of these two.

An additional point of departure is that I go on to calculate the "standard error of the average" as described in HoA ch 10 and use *that* as my +/- error value for the bin.

See discussion of this here also:

https://sourceforge.net/apps/mediawiki/vstar/index.php?title=Talk:Computing_Means

Now, that would result in mean curves having:

a. non 1-sigma error bars
b. much smaller error bars

So... Are we saying that I should stop at SD and just use that as my +/- error bar value? If so, it appears that I have abused/misunderstood the role of the "standard error of the average" in this context.

• David Benn - 2010-03-24
• priority: 7 --> 9

• David Benn - 2010-03-26

The latest from Matt is that standard error of the average is suitable here. Further input (e.g. from Aaron, Mike, others) would be beneficial.

• Aaron Price - 2010-03-27

The debate between stdev and SE is philosophical with no right answer. My understanding is that when you have a large sample size stdev is preferred and when you have smaller samples, SE is preferred. I think the only important thing is that the display clearly describe whether it is stdev or SE.

I originally meant for the LCG to show stdev, but maybe I messed up and coded SE. Don't use it as a test. (I'll look at the code for it when we roll out the new web site.)

Regarding this question. I'd always go with Matt's suggestion over anyone else's. As he likes to say: "Math is hard!"

• Aaron Price - 2010-03-27

side note: stdev in the LCG is computed using perl math libraries, so I bet it is accurate. The difference is likely from how the bin range is determined.

• David Benn - 2010-03-27

Okay, thanks for the input Aaron.

You guys are the domain experts and obviously that's why I'm asking for your input. My concerns re: std dev vs std err of the average are:

1. Error bars will be much smaller with the latter. For example, let's say you have a std dev of 0.12 for a set of values whose mean is calculated for a bin. If the bin size is 10 (number of data-points in the bin), the std error of the average is 0.12/sqrt(10) or 0.038. So +/-0.12 vs +/- 0.038. Pretty different error bars visually.

2. In commentary on Citizen Sky we see frequent references to 1-sigma error bars, in particular in looking at LCG plots. This relates to the probability of an actual value falling within the +/- std dev error bar value range assuming a bell curve, right? Now, as you say in your CS post on light curves in which you mention 1-sigma: "A good rule of thumb is to see if you can draw a horizontal and straight line between the error bars. If you can, it means there is no real variation in the data. If you can't, then the variation you see is more likely than not to be real.". As per 1 above, if the error bars are small, visually determining this is not easy.

Perhaps if you could compare a mean curve in VStar and its error bars against an LCG plot, or what you generally would expect, this might provoke more discussion.

There is some interesting commentary about the interpretation of error bars and the importance of argument here:

http://kicp.uchicago.edu/~kksm/lecture_6.pdf

Re: the side note, when you say Perl math libraries, do you mean something other than the intrinsic Perl math operations (IEEE 754 floating point), such as a CPAN arbitrary precision library? Java also uses IEEE 754 floating point libs. There are certainly all kinds of subtleties associated with FP math as I'm sure you know. Having said that, I managed to get the same results with my recent Java port of AAVSO DC DFT Fortran code, as the original Fortran code. But that probably has everything to do with the underlying, now-ubiquitous IEEE 754 nature of the operations.

Re: bin size, in VStar, the user can specify the bin size of course. :)

• David Benn - 2010-03-27

By the way, I know we're saying we should not hold LCG up as the golden reference, but I mention this because we make frequent reference to it on CS.

• Aaron Price - 2010-03-27

If it was just me, I'd go with stdev because that is what I've seen most often in the physical sciences. However, I want to reiterate that Matt is the closest we have to an expert on this issue and feel that we should defer to his decision.

The perl library I use is Statistics::Descriptive. LCG also allows the user to set a bin size. But the start and end points of where the bin counting begins may not be correct. GNUPlot has a weird issue where it sometimes likes to plot the first point of a plot DIRECTLY on the first y-axis and the last point DIRECTLY on the x-axis. This makes them nearly impossible to see on the plot and generated lots of e-mails to me along the lines of "where is my data point?". So to fix it I have had to add kludges, which may be affecting how bins are calculated. Plus, I'm a terrible programmer. So never use any of my software for any diagnostic. I'm hoping with my pending PhD that I'll be out of the programming business forever. :)

• David Benn - 2010-03-28

Thanks for the further advice Aaron. :)

• David Benn - 2010-03-28

Fixed mean JD/phase calculation and added a unit test. Assuming we really should be using std error of the average instead of std deviation, then this tracker can be marked as fixed. Marking it as such.

http://vstar.svn.sourceforge.net/viewvc/vstar?view=rev&revision=459

Sara, close this if happy.

• David Benn - 2010-03-28
• status: open --> open-fixed

• Aaron Price - 2010-03-31

FYI: I just consulted with Grant Foster and he agrees that we should use standard error instead of stdev. So I plan to adjust the LCG accordingly.

• David Benn - 2010-05-15

Wow, this was a great discussion reading back over it. Thanks.

Are we happy to close this tracker? Matt's info recorded below under the 2010-03-25 00:40:18 CST entry is what VStar does.

• David Benn - 2010-05-15
• assigned_to: david_benn --> aaronprice808

• Sara Beck - 2010-05-20

Thanks for staying with this one, David!

• Sara Beck - 2010-05-20
• status: open-fixed --> closed-fixed