#Introduction

I flew a lot of laps with my drone at home and timed them automatically: there were 8.5 s laps, a lot of 10 s laps, a fair number of 15 s laps, and very few 30 s and 60 s laps. Of course, I could just plot a histogram to visualize these results, but I wanted some kind of smooth curve, so that I wouldn't have to pick the bin width of the histogram, and so that this smooth curve would take into account how close the results are to each other.

I tried googling it and found nothing. Or I found something but didn't understand a thing (like this paper: it started with a bang well, but then wandered off somewhere and lost me). So I had to invent everything myself. I sat over a notebook for a bit, thought it through, then coded for a while and got a result (not right away, of course). In this post I'll explain how to approximate a one-dimensional distribution with a polynomial of degree N. The resulting algorithm has complexity $O(m)$, where $m$ is the number of points.

#Derivation

We have a lot of data points: $D_i$, $M$ of them in total. We want to approximate their distribution with a polynomial $P(x) = \sum_{i = 0}^n c_i x^i$.

What does it mean to approximate a distribution? It means that on any interval $[a_k, b_k]$ the integral of this function approaches the density of the data ($count_{a_k}^{b_k} D$ is the number of points of $D$ that fall into the interval $[a_k, b_k]$):

$$\int_{a_k}^{b_k} P(x) dx \approx \frac{count_{a_k}^{b_k} D}{M}$$ $$\int_{a_k}^{b_k} \left[\sum_{i = 0}^n c_i x^i \right] dx - \frac{count_{a_k}^{b_k} D}{M} \approx 0$$ $$\sum_{i = 0}^n c_i \cdot \frac{a_k^{i+1} - b_k^{i+1}}{i+1} - \frac{count_{a_k}^{b_k} D}{M} \approx 0$$

And if we take all the intervals, and there are K of them:

$$\sum_{k = 0}^K \left[ \sum_{i = 0}^n c_i \cdot \frac{a_k^{i+1} - b_k^{i+1}}{i+1} - \frac{count_{a_k}^{b_k} D}{M} \right]^2 \to \min$$

Here I squared the difference so that it's always greater than 0, and so that summing over different intervals gives a positive number. And so that the final sum over all intervals doesn't equal zero but tends to some minimum. In this equation we want to find the unknowns $c_i$.

If we expand the brackets and collect like terms, we get an equation of the following form:

$$a_{0,0} + a_{0,1} c_0 + \dots + a_{1, 1} c_0^2 + a_{1, 2} c_0 c_1 + \dots a_{n+1,n+1} c_n^2 \to min$$

You can see that this is a quadratic equation with many unknowns. And the minimum of such a function can be found directly: take the derivative and set it to 0:

$$\left\{ \begin{aligned} &-a_{0, 1} = a_{1, 1} c_0 + a_{1, 2} c_1 + \ldots + a_{1, n+1} c_n \\ &-a_{0, 2} = a_{2, 1} c_0 + a_{2, 2} c_1 + \ldots + a_{2, n+1} c_n \\ &\ldots \\ &-a_{0, n} = a_{n, 1} c_0 + a_{n, 2} c_1 + \ldots + a_{n+1, n+1} c_n \\ \end{aligned} \right.$$

So all we have to do is compute some coefficients $a_{i, j}$ and then solve the system of linear equations (SLE), and we get the coefficients of a polynomial that honestly approximates the given distribution.

#Algorithm

Input:

• $D_i$: an array of M values whose distribution we want to approximate.
• $N$: the degree of the resulting polynomial.
• $K$: the number of intervals.

Obviously, the degree of the polynomial should be small (5-10), so as not to "overfit" the data and end up with a plot that equals the histogram or something like that.

The number of intervals shouldn't be very large either. I figured that $5\cdot M$ would be plenty. In practice it turned out that $N \cdot M$ is optimal.

So the whole algorithm consists of two stages:

• Iterating over $K$ intervals to build the matrix and the vector.
• Solving a system of linear equations of size $N$, which is not very large.

There are two ways to iterate over the intervals:

• Randomly
• Systematically

In the random case, the most expensive operation is counting the number of points that fall into the interval, which takes $O(\ln n)$ if the data is sorted.

If instead we iterate over the intervals systematically, from the minimum to the maximum, we can use the two-pointer technique and know the number of points in the current interval in $O(1)$. And since the number of intervals is $O(n)$, this is the optimal option.

Altogether, the whole algorithm looks like this in pseudocode:

vec = /* 1D vector of size n+1 */
mat = /* 2D matrix of size n+1 */
minv = D[0];
maxv = D[-1];
seg_count = sqrt(K);
step = (maxv - minv) / seg_count;

apos = 0;
for ai in 0..seg_count+1 {
	a = minv + ai * step;
	while apos < M && a > D[apos] { apos += 1; }

	bpos = apos;
	for bi in ai+1..seg_count+1 {
		b = minv + bi * step;
		while bpos < M && b > D[bpos] { bpos += 1; }

		count = (bpos - apos) / M;

		for i in 0..n+1 {
			vec[i] += 2 * count * (b^(i+1) - a^(i+1))/(i+1);
		}

		for i in 0..n+1 {
			for j in 0..n+1 {
				to_add = (b^(i+1) - a^(i+1))/(i+1) 
				       * (b^(j+1) - a^(j+1))/(j+1);

				mat[i, j] += to_add;
				mat[j, i] += to_add;
			}
		}
	}
}

x = solve_SLE(vec, mat)

#Web interface

Links:

• Online interface in a separate tab
• Source code

I put in some effort and coded a web interface where you can play with this algorithm. And anyway, what's the point of talking about some algorithm you came up with if there's no code? If that paper had at least had code, never mind a web interface, maybe I wouldn't have had to invent anything myself...

My drone flight data is already loaded in there. Yes, I flew 300 laps in a day, and what are you going to do about it 😁. You can see how much my skill on this track leveled up the next day: faster times became much more frequent.