1
$\begingroup$

I'm working on an algorithm (written in Python/Cython, but it reads like pseudo-code) that estimates the gradient of each point in noisy data, using a variable window size. It's working very well, but it seems that the algorithm is limited by the regression part. Here is what I use:

cdef double regression(np.ndarray[DTYPE_t, ndim=1] data, np.ndarray[DTYPE_t, ndim=1] time, unsigned int leftlim2, unsigned int rightlim2):

    # declaring some variables
    cdef unsigned int length, j
    cdef double x, y, sumx, sumy, xy, xx, result, a, b, invlen

    # resetting values
    length = 0
    sumx = 0
    sumy = 0
    xy = 0
    xx = 0

    # doing a loop from the left limit to the right limit
    for j from leftlim2 <= j < rightlim2: 
        x = time[j]
        y = data[j]
        sumx += x   # the sum of x
        sumy += y   # the sum of y
        xy += x*y   # the sum of x*y
        xx += x*x   # the sum of x^2

    # estimating the best-fit slope
    length = rightlim2 - leftlim2
    invlen = 1.0/length
    a = xy-(sumx*sumy)*invlen
    b = xx-(sumx*sumx)*invlen
    result = a/b
    return result

Inputs:

  • vectors/arrays of the data and time that was measured during an experiment. The data array contains noisy data of, for example, applied force, the time array contains equally spaced time recordings (0.1s, 0.2s, 0.3s, etc.)
  • the left and right limits of how much data has to be included for the regression, provided as indices (i.e. the data used for regression is given by data[leftlim2:rightlim2])

Output: the slope of a straight line (the a in y = a*x + b) approximating the dataset.

I'm only interested in the slope, not in the intercept, hence the use of a loop rather than regression using matrix-vector multiplications, which is numerically inconvenient (i.e. slow). I was wondering if anyone knows a way to increase the efficiency of the regression, without sacrificing accuracy. Perhaps there's a way to exploit the equal spacing of the time array? For example, the sum of x can be written as 0.5*dx*(l^2 - l) + x[0]*l where dx is the time step and l = rightlim2 - leftlim2. Is there a way to represent the sum of x*y or x*x in a similar fashion? (I have tried and failed :( )

$\endgroup$
1
  • $\begingroup$ If you are still interested in an answer, I highly recommend migrating the question to Computational Science. You can flag the post to ask a moderator for migration. $\endgroup$ – Post No Bulls Dec 30 '13 at 6:07
1
$\begingroup$

You want to fit a model Y = a + b X using linear regression and you are only concerned by the value of the slope "b".

From the normal equations, you have b = (S1 * S2 - N * S3) / (S2 * S2 - N * S4) where N is the number of data points, S1 the sum of the Y(i), S2 the sum of the X(i), S3 the sum of the X(i) * Y(i) and S4 the sum of the X(i) * X(i).

Now, you want to take advantage that the X(i) are equally spaced; this means that
X(i) = X(1) + (i - 1) dX
where X(1) is the first value and dX the step size. This allows you to directly compute the value of sums S2 and S4 using the fact that the X(i) are in arithmetic progression. Using this fact leads to

S2 = N * X(1) + N * (N + 1) * dX / 2

S4 = N * X(1) * X(1) + N * (N + 1) * X(1) * dX + N * (N + 1) * (2 * N + 1) * dX * dX / 6

If you do not know the number of data points but only the first value X_first, the last value X_last and the step size dX, you have

N = 1 + (X_last - X_first) / dX

and X(1) = X_first.

Similarly, it could be possible to expand S3 and write

S3 = X(1) * S1 + dX * Sum[(i - 1) Y(i),{i,1,N}]

but I do not think this would induce any savings in calculation time.

By the way, you can also get the intercept for free since a = (S1 - b * S2) / N

$\endgroup$
4
  • $\begingroup$ Thanks for your answer. I already found the same result for S2, but your analysis of S4 is new to me. $\endgroup$ – MPA Dec 30 '13 at 16:29
  • $\begingroup$ @MPA. The problem you posted is still more interesting if you consider a polynomial model since you will find the sum of the X(i), the sum of X(i)^2, ..., the sum of X(i)^n in the normal equations. Since the X(i) are in arithmetic progression, all these sums have simple closed forms. $\endgroup$ – Claude Leibovici Dec 31 '13 at 9:04
  • $\begingroup$ One can set X(1) to 0 or some other magic values to get simpler answer for the slope b. $\endgroup$ – fchen Apr 1 '15 at 12:58
  • $\begingroup$ Since the points are equally spaced it is easy to choose an offset on X such that $\sum X = 0$; then two of the terms disappear completely and you get $b=\frac{\sum xy}{\sum{x^2}}$. That is computationally much easier / faster (especially since, for known window, the denominator will be constant. $\endgroup$ – Floris Sep 19 '17 at 22:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.