#Введение

Я налётывал дома на дроне очень много кругов, и автоматически замерял их время, были круги по 8.5с, много кругов по 10с, и какое-то количество по 15с, и очень мало по 30с, 60с. Конечно, можно построить гистограммку, чтобы визуализировать эти результаты, но я хотел какую-то гладкую кривую, чтобы мне не приходилось подбирать длину ведра в гистограмме, и чтобы эта гладкая кривая учитывала насколько близко находятся результаты рядом друг с другом.

Пытался нагуглить - не нашёл. Или находилось что-то, но я ничего не понял (как, например, эта статья - начали за здравие хорошо, а потом ушли куда-то и я не понял). Поэтому пришлось изобретать всё самому. Посидел немного над тетрадкой, подумал, затем покодил и получил результат (не сразу, конечно же). В этой статье расскажу как аппроксимировать одномерное распределение полиномом N-ой степени. Итоговый алгоритм получился по сложности $O(m)$, где $m$ - количество точек.

#Вывод

У нас есть много данных: $D_i$, всего их $M$. Мы хотим аппроксимировать их распределение полиномом $P(x) = \sum_{i = 0}^n c_i x^i$.

Что значит аппроксимировать распределение? Это значит что на произвольном отрезке $[a_k, b_k]$ интеграл данной функции стремится к плотности данных: ($count_{a_k}^{b_k} D$ - количество точек $D$, попавших в отрезок $[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$$

И если взять все отрезки, количество которых K:

$$\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$$

Здесь я взял разницу в квадрат, чтобы она всегда была больше 0, и чтобы в сумме на разных отрезках получалось положительное число. И чтобы итоговая сумма по всем отрезкам не равнялась нулю, а стремилась к некоторому минимуму. В этом уравнении мы хотим найти неизвестные $c_i$.

Если раскрыть скобки, и сложить слагаемые, то получится уравнение следующего вида:

$$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$$

Видно что это уравнение второй степени с множеством неизвестных. А найти минимум такой функции можно напрямую: взяв производную и найдя её 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.$$

Итого, надо просто посчитать какие-то коэффициенты $a_{i, j}$ и затем решить СЛАУ, и у нас получятся коэффициенты полинома, который честным образом аппроксимирует данное распределение.

#Алгоритм

Входные данные:

• $D_i$ - массив значений размера M, распределение которых надо аппроксимировать.
• $N$ - степень итогового полинома.
• $K$ - количество отрезков.

Степень полинома, очевидно, надо брать небольшую (5-10), чтобы не "переобучаться" на данных, и не делать график, который равен гистограмме или ещё чему-то такому.

Количество отрезков нужно брать тоже не очень большое. Я прикинул, что $5\cdot M$ вполне достаточно. В итоге на практике оказалось что оптимальным получается $N \cdot M$.

Значит весь алгоритм состоит из двух этапов:

• Перебор $K$ отрезков и получение матрицы и вектора.
• Решение СЛАУ размером $N$, то есть не очень большой.

Перебирать отрезки можно двумя способами:

• Рандомно
• Организованно

В рандомном случае самой ресурсоёмкой операцией получается подсчёт числа точек, попадающих на этот отрезок, что занимает $O(\ln n)$, если данные отсортированы.

Если же пытаться перебирать отрезки организованно, от минимума к максимуму, то можно использовать метод двух указателей, и за $O(1)$ знать количество точек, находящихся на данном отрезке. А так как количество отрезков равно $O(n)$, то это - оптимальный вариант.

Итого весь алгоритм выглядит так в псевдокоде:

vec = /* одномерный вектор размера n+1 */
mat = /* двумерная матрица размера 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)

#Веб-интерфейс

Ссылки:

• Онлайн-интерфейс в отдельной вкладке
• Исходный код

Я постарался и накодил веб-интерфейс где с этим алгоритмом можно поиграться. И вообще какой смысл выпускать говорить о том что ты придумал какой-то алгоритм, без кода? Вот если бы в той статье был хотя бы код, я не говорю про веб-интерфейс, может я бы сам ничего не придумывал бы...

Там уже вставлены данные моих полётов на дроне. Да, я отлетал 300 кругов за день, и что вы мне сделаете 😁. Видно как сильно прокачался навык на данной трассе на следующий день, что более быстрые времена стали намного чаще.