A curve-fitting function in Eidos would be very useful for plotting

[bhaller]

It's common to have (x,y) scatterplot data and want to fit a nice best-fit curve through that data for visualization. It'd be great if Eidos had a way to do this. There are many ways of curve-fitting, of course. Some guarantee that the fitted curve actually passes through all the points in question; that is not what I have in mind (although doubtless useful in some situations). What I have in mind is fitting a curve, preferably without any predetermined functional form, that is in some sense a "best fit" to the data without responding too sensitively to any one data point. AFAICT this is quite complicated and requires iterative numerical approximation and so forth. Not necessarily out of the question, but I'd prefer not to have to figure it all out for myself, so I'm hoping for a simple C++ implementation of something open-source that I can just throw into Eidos.

I found https://github.com/ttk592/spline but I think (?) it only fits curves that exactly pass through all points; not what I want.

Then I found https://software.nasa.gov/software/MSC-26693-1 which is on GitHub at https://github.com/nasa/polyfit. It looks pretty reasonable? It only fits polynomials, of a degree that the user specifies, so that's not ideal; I'd rather have something that fits a simple function to data that follow a simple pattern, and fits a complex function to data that follow a complex pattern, automatically. I suppose I could add a front end to this library that tries polynomials up to some limit order (square root of number of data points or some such, with a max of ~20?), and chooses the lowest-order polynomial with a "good" R2 according to some heuristic. But ugh; I was hoping to avoid that kind of heuristics. It's also not ideal that its license is unclear and it says "The program has been reviewed and approved by export control for public release. However some export restriction exists. Please respect applicable laws." That is really ambiguous and I don't know whether I can legally incorporate it in SLiM. :-<

I've found other libraries that can do more than polyfit, but there is a huge step upward in complexity, such that they are not usable in Eidos; e.g., http://ceres-solver.org.

See #527 for the context in which this originally arose.

Ideas, @petrelharp or @vsudbrack or anybody else?

[vsudbrack]

The most common one I know is the LOESS
https://www.wikiwand.com/en/articles/Local_regression
and since Eidos now has inverse matrix, it can calculate it (although idk how many points it will be able to do at once because inverting large matrices is costly). This feature would be very exciting :)
It is implement on R with the function loess() from the package stats.

[bhaller]

OK, it took some sleuthing, but I think I found a usable C implementation of LOESS! The original paper is here:

https://sites.stat.washington.edu/courses/stat527/s13/readings/Cleveland_JASA_1979.pdf

and the C code is here:

https://www.netlib.org/a/dloess

I'll see if I can do something with this. :-> Thanks for the suggestion, it looks like exactly what I want. :-> It's surprising that there are (as far as I can find) no more modern open-source LOESS packages, on GitHub or elsewhere, in C or C++. I found ones in Python and even JavaScript, but no C/C++ except this 1979 (!) paper. Sheesh.

[gregorgorjanc]

R's stats:loess calls stats:::simpleLoess, which calls C_loess_raw and C_lowesp (this last one is Fortran), but there is quite a bit of R code to prepare inputs. The ancestor of LOESS is LOWESS, which in R is in stats:lowess, which calls C_lowess (with very minimal input prep).

[bhaller]

Yeah, I could try to extract C loess code from R or Python implementations if necessary, but usually that kind of work is extremely gross. :-> Usually the core function calls out all over the place to other C code, and you have to track down all those connections and surgically cut out what you need. I'm hoping for a self-contained C or C++ solution that I can just drop in.

[bhaller]

Well, the Cleveland 1979 loess code is a dead end, I think. It turns out to be a mix of C and Fortran (pages and pages of dense code with single-letter variable names and no comments), and has external dependencies on LINPACK (an ancient Fortran math library) and BLAS (which possibly I could connect up to the BLAS functions in the GSL, but the rest is a showstopper). I think the loess angle is looking difficult; it seems like it is quite a mathematically complex method, and thus tends to rely on pretty extensive external math library support.

I think this is looking like a substantial project. Would perhaps be a good project for a grad student with good math and CS chops. I'll tag it as such and move on, for now. If anybody reading this has ideas for a clear path forward, please speak up. :->

[petrelharp]

I suspect a usable implementation of loess would be, indeed, a Project. The easiest ways I can think of to do this would be

1. local averaging: loess is a local averaging method, but kinda fancy; but there are many ways to produce a "curve which is averages in windows" that would get something just as good for visualization; either (a) at a grid of x values, average the nearest k y values; or (b) at a grid of x values, average all y values within some distance
2. linear model: the formula for least-squares fitting of a linear model is straightforward (with matrix algebra); you could fit a k-th order polynomial to the data (better would be splines, but then you'd have to implement those...)

[bhaller]

@petrelharp Are there functions in R that you'd recommend that I look at, that I would want to pattern Eidos after in this area? I mean, obviously there are a MILLION ways to fit a model to data in R. :-> I'm asking if there is something that you would particularly recommend as being simple, with a good API and an approach such as the ones you mention that would be relatively tractable. :-> Or if you ever feel inspired and with time on your hands, maybe you could explore doing this. Probably either of the approaches you mention would now be do-able in pure Eidos code as a user-defined function, given the new matrix algebra stuff I added to Eidos in SLiM 5.1...? So it mostly just requires some math chops to make it happen, I think. It'd be a nice addition.

[petrelharp]

Hm. I can't find anything that's quite so dead-simple in R, because there's so many not-quite-as-simple-but-better smoother (eg lowess or loess or smooth.spline) The simplest I can find is supsmu, based on this FORTRAN code. From the description:

     ‘supsmu’ is a running lines smoother which chooses between three
     spans for the lines. The running lines smoothers are symmetric,
     with ‘k/2’ data points each side of the predicted point, and
     values of ‘k’ as ‘0.5 * n’, ‘0.2 * n’ and ‘0.05 * n’, where ‘n’ is
     the number of data points.  If ‘span’ is specified, a single
     smoother with span ‘span * n’ is used.

     The best of the three smoothers is chosen by cross-validation for
     each prediction. The best spans are then smoothed by a running
     lines smoother and the final prediction chosen by linear
     interpolation.

However, that's not great because (a) it's O(n) when really you want something O(number of points you need to draw a good line); and (b) it gets fancy with cross-validation.

Here's a very simple "give me the mean of the nearest p% of points" smoother:

x <- runif(1000)
y <- 2 * x + rnorm(length(x), sd=0.2)

smooth1 <- function (x, y, p) {
    # note this could be done easily by sorting x so then
    # xk would be a vector like (1, 1, ..., 1, 2, 2, ..., 2, ...)
    k <- round(length(x) * p)
    breaks <- quantile(x, seq(0, 1, length.out=k+2)[-c(1, k+2)])
    xk <- findInterval(x, c(-Inf, breaks, Inf))
    stopifnot(all(range(xk) == c(1, k+1)))

    xx <- tapply(x, xk, mean)
    yy <- tapply(y, xk, mean)
    return(list(x=xx, y=yy))
}


xy1 <- smooth1(x, y, p=0.01)
xy2 <- smooth1(x, y, p=0.1)

plot(x, y, pch=20)
lines(xy1, lwd=2, col='red')
lines(xy2, lwd=2, col='green')

[bhaller]

@petrelharp I see. The code you supply looks very useful, thanks! Since I just shipped 5.1 it'll (hopefully) be a while before the next release, so this issue might stagnate for a while, but I'll get around to it. :-> Maybe I'll wrap both your smoother and the PolyFit polynomial smoother into one function, with a "method" parameter or some such, allowing for future expansion if I ever do manage to get fancier things like loess or splines in there. :-> I think what you've provided should already be quite useful!

[bhaller]

I've just discovered that the GSL has support for B-spline smoothing, which would fit the bill here quite well I think:

https://www.gnu.org/software/gsl/doc/html/bspline.html#overview

I don't time right now to look into it further, but I suspect this is a very good path forward.