add a function like R's filter() function for computing running averages for plotting etc.

[bhaller]

@vsudbrack suggests that a running average function would be useful; he is doing:

function (float)running_mean(lif y, [integer$ window = 3])
{
    half_window = asInteger(floor(window / 2));
    y_mean = rep(NAN, length(y));
    for (i in half_window:(length(y) - half_window - 1))
        y_mean[i] = mean(y[(i - half_window):(i + half_window)]);
    return y_mean;
}

R's filter() function would be perhaps good functionality to aim for, at least in part. Food for thought.

[petrelharp]

A more efficient implementation would do something like (not tested, not returning the right length vector, but you get the idea):

n = length(y);
ys = cumSum(y);
y_mean = (ys[window:(n-1)] - ys[0:(n-window)])/window;

I suspect that the general filter method in R uses a fast Fourier transform under the hood, since otherwise this sort of operation can be expensive.

[bhaller]

Hmm. Well, I'm not going to add FFT code to SLiM, I don't think...? :-> My impression is that this functionality is being requested for custom plotting purposes, so probably we don't need to go over the top here; I think @vsudbrack is just looking to smooth out plotted data to decrease noise. Shall I just add a running average function? I could add the implementation that @vsudbrack supplied, perhaps renamed calcRunningMean(); or @petrelharp perhaps you could supply a complete working example that is more efficient, based upon the code you gave? (I guess the problem with it is how the ends of the vector are treated? I don't know what would be appropriate there; does the code revert to Vitor's algorithm for the first and last window-1 entries?) Efficiency would indeed be nice, since it is nice to be able to update custom plots in every tick and have them redisplay without slowing the model down appreciably.

[bhaller]

The R doc notes:

convolve(, type = "filter") uses the FFT for computations and so may be faster for long filters on univariate series, but it does not return a time series (and so the time alignment is unclear), nor does it handle missing values. filter is faster for a filter of length 100 on a series of length 1000, for example.

So it sounds like FFT is not needed, even if we do want to have more functionality than just a simple running mean. Do we? We could certainly replicate the R filter() function with a little elbow grease; I just suspect that is overkill. So, step one: I need clarification as to what is really needed here, without requesting more than will actually be useful for people, since I have limited time.

[bhaller]

OK, having pondered this, I think I understand it better and I propose to implement a modified version of R's filter() function:

(float)filter(numeric x, float filter, [lif$ outside = F])

x is the sequence (typically a time series) to be filtered.

filter is a vector of float filter coefficients in reverse time order, following R; for a simple central moving average of width W, just pass rep(1,W)/W. R suggests but does not require that filter be an odd length, but I will require it, so that the filter is centered correctly; I can always relax that requirement later.

outside is not in R. In R, values outside x are (AFAICT) always treated as NA, and so the result of the filter has NA values at one or both ends; it only computes the filter results for positions where it has complete data, and otherwise produces NA. In Eidos I plan to offer a couple of options, governed by this outside parameter. If outside is F, values outside x will be considered NAN (since Eidos doesn't have NA), and the result of calculations involving those outside values will therefore be NAN; this is basically the same as R's behavior, but with NAN. If outside is T, values outside x will be excluded from the calculation (the filter value covering that position will be considered to be 0), and the other values in the filter will be adjusted to sum to the same total as before to compensate, giving the positions inside x more weight. Finally, if outside is integer or float, that value will be used as the value of x for all positions of outside x; one might pass an expected value or mean value in this way, to be used for all outside positions.

R has further parameters, for which I'll assume a single value and will not include them in the Eidos API for the time being, but they could be added later:

method in R can be convolution or recursive; I will support only convolution.

sides in R can be 1 (past values only) or 2 (centered); R's default is 2, and I will support only that mode.

circular in R can allow x to be treated as a circular buffer; I will not support that.

init in R is used only for recursive filters, which I will not support, but it is kind of similar to my outside parameter.

So. My main question is whether outside is too mathematically objectionable, and I should just leave it out. It feels useful to me; it's annoying getting NANs at the start and end instead of the function just computing an estimate based upon fewer data points. For now I will leave it out, and just implement the default outside=F case that produces NANs.

My other question is whether any of the functionality that I'm leaving out is deemed important enough to include. I don't see much utility, in SLiM, for recursive filters, or for the side=1 "past values only" behavior, or for the circular=T circular buffer behavior, or for supporting a filter with an even number of elements. Any dissent?

[bhaller]

OK, this has been added. See above for comments. I did not do the outside parameter for now; I decided to leave it as a strict subset of R's filter() function until someone asks me to extend it further. I tried using it with a plot, to smooth a plotted curve, and it worked great. :->

[petrelharp]

A bit late, but I like outside.

My suggestion was not that you implement FFT, but rather that you only support the running mean case (so, like R's runmean) since that can be done efficiently without FFT. Since you've done the full filter I think you should add a warning to the documentation saying that the runtime of this operation is (length of vector) * (length of filter) so will be very slow for large filters.

[bhaller]

A bit late, but I like outside.

OK, I'll work on it. :->

My suggestion was not that you implement FFT, but rather that you only support the running mean case (so, like R's runmean) since that can be done efficiently without FFT. Since you've done the full filter I think you should add a warning to the documentation saying that the runtime of this operation is (length of vector) * (length of filter) so will be very slow for large filters.

Sure. A speed test:

x = runif(1000000);

f = rep(1/999, 999);
f = runif(999, 0.99/999, 1.01/999);

s = clock();
y = filter(x, f);
e = clock();

catn("elapsed: " + (e - s));

With the first filter, a simple running mean (which the function detects and optimizes), the runtime is 0.006121 seconds. With the second filter (a simple running mean with slight variation, which is then not optimized)), the runtime is 1.10188 seconds. This is for a filter length of 999, which is the longest filter allowed by the code at present. (There's no strong reason to impose a max, but it seemed like it might catch bugs, and a filter length of 999 already seems longer than anybody is likely to use in practice...? But I could remove the max if it is not liked. :->)

[bhaller]

OK, outside is added. outside=T works very nicely for plotting! :->