МНК, статистика, дискретность

[B@R5uk]

МНК — вещь очень полезная и удобная. Более того, если подгоняемые величины имеют нормальное распределение, то МНК кроме всего прочего естественным образом следует из принципа максимального правдоподобия, то есть имеет статистическое обоснование.

Однако, подгоняемые величины далеко не всегда подчиняются распределению, близкому к нормальному. Более того, распределение может даже не быть непрерывным. Например, любые показания, снимаемые с прибора со шкалой (линейка, вольтметр и так далее) имеют дискретность. Хорошо, если в показаниях естественным (или принудительным) образом присутствует шум (величиной, большей цены деления), и их, набрав подходящее количество, можно усреднить. Полученные средние будут иметь распределение более-менее похожее на нормальное (центральная предельная теорема), и их можно будет без зазрения совести подставлять в МНК.

Но и так бывает тоже далеко не всегда. Вопрос: что в этом случае делать?

Вот, например, предельно утрированный случай. Имеется изображение линии и требуется оценить коэффициент наклона этой линии к горизонтали (с погрешностью):



Можно, конечно, перебрать все пиксели линии, составить пары координат икс-игрек и подставить их в МНК. Но что-то подсказывает мне, что такой подход будет неправильным. В частности, интуитивно понятно, что ключевую информацию в изображении несут те участки, где линия переходит с одной строчки на другую. Остальные же пиксели просто подтверждают нам, что, да, мы имеем дело именно с прямой; но информационной нагрузки в смысле вклада в вычисление коэффициента наклона прямой (и его погрешности) эти пиксели не имеют. Но одно дело интуиция, а другое — как это формализовать?

[wrest]

Я не понял насчёт линии. У линии просто меряете $\Delta x / \Delta y$

[ilghiz]

B@R5uk в сообщении #1567757 писал(а):

Вот, например, предельно утрированный случай. Имеется изображение линии и требуется оценить коэффициент наклона этой линии к горизонтали (с погрешностью):

Можно, конечно, перебрать все пиксели линии, составить пары координат икс-игрек и подставить их в МНК.

А чем вам этот подход не нравится? Зафитите каждую точку полиномом первой степени, получите два уравнения и два неизвестных с симметричной матрицей, и будет вам счастье. Можно коенечно total least squares применить, но конкретно в этом примере разницы в точности поиска наклона не будет.

[B@R5uk]

ilghiz, вопрос не в том, что нравится, а в том, как правильно.

wrest в сообщении #1567760 писал(а):

У линии просто меряете $\Delta x / \Delta y$

А можно формулу конкретную? Вот, например, хочу, чтобы точность была максимальна, поэтому беру две крайние точки линии (дельты будут максимальны, точность выше): по иксу дельта равна 129, по игреку — 8, получается наклон 0,0620 относительная погрешность пусть будет 1/129, абсолютная тогда получится 0,0005, то есть результат

$$k=0,0620\pm 0,0005$$

Правильно? Потому что, если я возьму дельты по ступенькам линии, то у меня они получатся равными 109 и 7 соответственно, и результат будет

$$k=0,0642\pm 0,0006$$

Какой из них будет точнее? Погрешности я практически из воздуха взял; даже результаты в пределах этих погрешностей не совпадают. Будет ли результат точнее, если я учту, что ступеньки имеют длину 8—15—16—16—15—16—16—15—13 (именно в таком порядке)? Если да, то как?

Из каких вообще соображений подходить к этим вопросам?

[wrest]

B@R5uk в сообщении #1567762 писал(а):

Будет ли результат точнее, если я учту, что ступеньки имеют длину 8—15—16—16—15—16—16—15—13 (именно в таком порядке)?

Да, ну там есть априорное знание что тангенс угла наклона будет $\dfrac{3}{15+16+16}\approx 0,064$ т.к. видно что паттерн повторяется.

[ilghiz]

Не, смотреть на крайние и на не крайние - не правильно, но у вас есть реально произвол, так как вам дана только часть этого отрезка. То есть по идее вы можете провести эту линию и через $8/129$ и через $3/47$ и через $7/109$. Все будет подходить, но МНК как раз даст наиболее точный результат.

Смотрите, классика

$$\min_{a,b} E_{a,b}; ~~ E_{a,b} = \sum_i ||a i + b - y_i||_2^2$$

$$\frac12 \frac{dE_{a,b}}{da} = a \sum_i i^2 + b \sum_i i = \sum_i i y_i$$

$$\frac12 \frac{dE_{a,b}}{db} = a \sum_i i + N b = \sum_i y_i,$$

$$y_0,...,y_{129} = 0,...0:8; 1,...,1:15; 2,...,2:16; 3,...,3:16; 4,...,4:15; $$

$$ 5,...,5:16; 6,...,6:16; 7,...,7:15; 8,...,8:13$$

$$ 723905 a + 8385 b = 46566 $$

$$ 8385 a + 130 b = 540 ~~ | ~~ -64.5 $$

$$a = 23472/366145 \simeq 0.0641$$

Хотя более правильно - выбросить первый и последний сегменты, так как у нас нет информации где они закончились, тогда

$$y_0,...,y_{108} = 1,...,1:15; 2,...,2:16; 3,...,3:16; 4,...,4:15; 5,...,5:16; 6,...,6:16; 7,...,7:15$$

и

$$a = 3371/53955 \simeq 0.0625$$

Теперь о точности. Надо взять максимально и минимально возможное значение от угла этой возможной линии. Тут только через total least squares. Глаза слипаются, не сейчас.

wrest в сообщении #1567770 писал(а):

Да, ну там есть априорное знание что тангенс угла наклона будет $\dfrac{3}{15+16+16}\approx 0,064$ т.к. видно что паттерн повторяется.

Не, боюсь, что тут не все так однозначно, так как паттерн может быть дробной последовательностью, в том числе 15-16-16-15-16-16-15-16-16-16-15-16-16 или 15-16-16-15-16-16-15-16-15-16-16

[wrest]

ilghiz в сообщении #1567771 писал(а):

Не, боюсь, что тут не все так однозначно, так как паттерн может быть дробной последовательностью,

Может, поэтому я взял два знака, а не три.

[B@R5uk]

ilghiz в сообщении #1567771 писал(а):

но МНК как раз даст наиболее точный результат

На основании чего?

ilghiz в сообщении #1567771 писал(а):

более правильно - выбросить первый и последний сегменты

Опять же, на основании чего? Как строго показать, что вот такое вот подмножество имеющихся точек не содержит дополнительной информации в дополнению к информации, содержащейся в остальных точках?

[Евгений Машеров]

Правильный подход - заметить, что отклонения измеренных значений от линии имеют равномерное распределение и считать методом максимального правдоподобия для этого распределения. Разумеется, это предполагает, что иной ошибки, кроме ошибки квантования, нет.
Практичный подход - плюнуть на "неоптимальность" и считать МНК, полагая, что это когда есть тяжёлые хвосты, МНК опасен, а если вообще хвостов нет, то вполне работает.
Ну и совсем рабоче-крестьянский - находим точки разрыва, полагаем, что истинное значение в них посредине, проводим линию через наиболее отдалённые точки.
Задача на студенческую курсовую - генерировать данные и сравнить точности этих подходов.

[wrest]

B@R5uk в сообщении #1567762 писал(а):

Погрешности я практически из воздуха взял;

Ну погрешности исходных данных там -- полквадратика чисто из-за их квантования. Вот из этого и исходить.

[B@R5uk]

Евгений Машеров в сообщении #1567777 писал(а):

Правильный подход - заметить, что отклонения измеренных значений от линии имеют равномерное распределение и считать методом максимального правдоподобия для этого распределения. Разумеется, это предполагает, что иной ошибки, кроме ошибки квантования, нет.

Пока рассматриваю утрированную задачу выше. Случайного шума перед квантованием нет. ММП будет означать, что для каждой экспериментальной точки значение модели лежит на некотором отрезке, то есть у нас на каждую экспериментальную точку будет два неравенства (равномерное распределение). Потому что в противном случае (модельное значение выходит за пределы отрезка) функция правдоподобия обратится в ноль.

Другими словами. Вот у меня модель прямой:

$$y=\lfloor ax+b\rfloor,\quad|a|\le\frac{1}{2}$$

для каждой экспериментальной точки $\left\{x_k,\;y_k\right\}$ по ММП я получаю пару неравенств:

$$0\le ax_k+b-y_k<1$$

В результате, искомые значения моих коэффициентов прямой — это множество, изображение которого на плоскости a—b является многоугольником, ограниченным парами прямых, соответствующих экспериментальным точкам.

Более, того. Вот у меня есть первая "ступенька" прямой на картинке. Координаты крайних пикселей этой ступеньки (нижний левый угол — $\left\{0,\;0\right\}$) будут $\left\{8,\;20\right\}$ и $\left\{22,\;20\right\}$, если я не обсчитался. Из этих экспериментальных точек получаются следующие неравенства на искомые коэффициенты прямой:

$$0\le 8a+b-20<1$$

$$0\le 22a+b-20<1$$

Помножим первую систему неравенств на $(1-\beta)$, а вторую — на $\beta$. До тех пор, пока $0\le\beta\le 1$ оба множителя неотрицательны, и знаки в неравенствах сохраняются. Поэтому результаты умножения можно сложить:

$$0\le (8+16\beta)a+b-20<1$$

Получается система неравенств, в которой численный коэффициент перед параметром a является произвольным числом из отрезка $[8,\;22]$. Этот набор неравенств покрывает в том числе и те, которые соответствуют "экспериментальным точкам"-пикселям из середины рассматриваемой "ступеньки" прямой. Другими словами, в силу линейности модели пиксели внутри ступеньки являются прямым следствием пикселей на концах ступеньки и дополнительной информации о коэффициентах модели не несут.

Евгений Машеров, спасибо большое! Наступило некоторое прояснение. Во всяком случае, ответ на вопрос про информационное наполнение нашёлся.

Евгений Машеров в сообщении #1567777 писал(а):

Задача на студенческую курсовую - генерировать данные и сравнить точности этих подходов.

Кстати, вот интересная же задача. Проблема только алгоритмизовать нахождение коэффициентов по ММП. Пока не представляю, как систематически обрабатывать кучу неравенств. Одно только представление данных (в виде многоугольника) уже что-то нетривиальное по сравнению с МНК.

[Евгений Машеров]

Проблема с равномерным распределением, что появляется неоднозначность.
Сдвигаем параллельно или с наклоном, но чтобы не выскакивало "за флажки" - не меняем функцию правдоподобия.

[B@R5uk]

Ну, на мой взгляд, неоднозначность — ничего нового. С зашумлёнными данными мы же тоже не знаем, какую величину приняла случайная составляющая в каждом конкретном сэмпле. Если трактовать функцию максимального правдоподобия (после соответствующей нормировки) как плотность вероятности того, что параметры модели имеют данное конкретное значение, то для каждого параметра можно получить его собственную функцию распределения, наивероятнейшее значение, среднее значение и дисперсию. Здесь то же самое, только здесь наивероятнейших значений — целая область, поэтому смысловую нагрузку несёт среднее значение. Ну и самое главное для меня: можно посчитать погрешность параметров!

-- 26.10.2022, 14:45 --

Поскольку я не умею решать системы линейных неравенств (в частности, не знаю алгоритма исключения лишних неравенств), я подошёл к проблеме методом научного тыка на графике:

Используется синтаксис Matlab M

%   URL: dxdy.ru/topic151133.html
%   File: line_fit.m
%   Date: 2022.10.26

clc
clearvars
format compact

wdata = [
       0    19
       7    19
       8    20
      22    20
      23    21
      28    21
      39    22
      54    22
      55    23
      69    23
      70    24
      85    24
      86    25
     101    25
     102    26
     116    26
     117    27
     130    27
];
dnum = size (wdata, 1);

flags = [
    2
    5
    14
    17
];

aa = (0.061 : 0.00001 : 0.066)';

color = abs (randn (2 * dnum, 3));
color = 0.8 * color ./ repmat (sqrt (sum (color .^2, 2)), 1, 3);
%color = zeros (dnum, 3);

plot (aa (1), wdata (1, 2) - wdata (1, 1) * aa (1), 'w')
hold on
legstr = cell (2 + numel (flags) + 1, 1);
legstr {1} = '';
count = 1;
for k = 1 : dnum
    if ~ismember (k, flags)
        continue
    end
    x = wdata (k, 1);
    y = wdata (k, 2);
    bb0 = y - x * aa;
    bb1 = y - x * aa + 1;
    %plot (aa, bb, 'Color', color (k, :))
    plot (aa, bb0, 'Color', color (2 * k - 1, :), 'LineWidth', 2)
    plot (aa, bb1, 'Color', color (2 * k - 0, :), 'LineWidth', 2)
    count = count + 1;
    legstr {count} = ['\{ ', num2str(x), ', ', num2str(y), ', 0 \}'];
    count = count + 1;
    legstr {count} = ['\{ ', num2str(x), ', ', num2str(y), ', 1 \}'];
end
hold off
xlim ([0.0636 0.0642])
ylim ([19.52 19.56])
legend (legstr)
grid on
 

В результате областью допустимых параметров модели прямой оказался четырёхугольник:

$$7a+b-19<1$$

$$23a+b-21\ge 0$$

$$101a+b-25<1$$

$$117a+b-27\ge 0$$

То есть информация о параметрах модели зашита всего в четырёх пикселях картинки (что, на мой взгляд, очень забавно; думал, будет больше). Ниже я пометил их красным:



Осталось посчитать распределение для параметра a, его среднее и дисперсию. Наиболее вероятное значение тоже, кстати, будет вполне конкретным числом, поскольку ФМП зависит от двух параметров, по одному из которых производится интегрирование, чтобы получить распределение другого.

[Schrodinger's cat]

Отбрасываем конечные отрезки не вошедшие полностью и считаем среднее арифметическое длин всех оставшихся горизонтальных отрезков.
Тангенс угла будет равен отношению одного пикселя (все орезки по вертикали различаются на 1 пиксель) к средней длине отрезка (количество пикселей в горизонтальном отрезке).
Правда с шумом уже не покатит.

Евгений Машеров в сообщении #1567777 писал(а):

Ну и совсем рабоче-крестьянский - находим точки разрыва, полагаем, что истинное значение в них посредине, проводим линию через наиболее отдалённые точки

Мы не знаем из какой точки грани пикселя вышла линия и в какую точку грани другого пикселя она вошла.
На представленном примере погрешность будет гулять по x.
А на линиях с большим наклоном, приближенных к вертикальным, наоборот погрешность будет гулять по y.

А так, согласен с ilghiz, МНК по всем пикселям даст наиболее точный результат. Считая, что координаты пикселя в его центре.
Уже как то заморачивался с подобными задачами.

[B@R5uk]

Досчитал я свою модельную задачку. Получается следующее:



Значение наклона прямой (с удвоенным стандартным отклонением в качестве погрешности для надёжности):

$$a=0,06386\pm 0,00019$$

Используется синтаксис Matlab M

%   URL: dxdy.ru/topic151133.html
%   File: line_fit_2.m
%   Date: 2022.10.26

clc
clearvars
format compact

al = 0.0636;
ar = 0.0643;

wdata = [
       7    19    1
      23    21    0
     101    25    1
     117    27    0
];

xd = wdata (:, 1);
yd = wdata (:, 2) + wdata (:, 3);

indx = [
    1   4
    1   3
    2   3
    2   4
];

vrtx = zeros (4, 2);
for k = 1 : 4
    vrtx (k, :) = [xd(indx (k, :)), ones(2, 1)] \ yd (indx (k, :));
end

aa = sort ([vrtx(1 : 3, 1); al; ar]);
pp = [0; 0; 2 / (aa (4) - aa (2)); 0; 0];
amean = sum (aa (2 : 4)) / 3;
adev = sqrt ((sum (aa (2 : 4) .^ 2) - aa (2) * aa (3) - aa (2) * aa (4) - aa (3) * aa (4)) / 18);

color = [
    0 0 0
    0.6 0.3 0
    0.7 0 0.6
    0 0 1
];

subplot (3, 1, 1 : 2)
fill (vrtx (:, 1), vrtx (:, 2), [0, 1, 1])
hold on
for k = 1 : 4
    a = [al ar];
    b = yd (k) - a * xd (k);
    plot (a, b, 'Color', color (k, :), 'LineWidth', 3)
end
hold off
xlim ([al ar])
ylim ([19.52 19.56])
grid on
legend ('Model parameters allowed values', 'Pixel \{7, 19\}', 'Pixel \{23, 21\}', 'Pixel \{101, 25\}', 'Pixel \{117, 27\}')
ylabel ('Parameter {\it b}')
title ('Parameters space')

subplot (3, 1, 3)
plot (aa, pp, 'Color', color (4, :), 'LineWidth', 3)
hold on
plot ([1 1] * amean, ylim, 'Color', color (3, :), 'LineWidth', 2)
plot ([1 1] * (amean - adev), ylim, 'Color', color (1, :), 'LineWidth', 2)
plot ([1 1] * (amean + adev), ylim, 'Color', color (1, :), 'LineWidth', 2)
hold off
xlim ([al ar])
grid on
xlabel ('Parameter {\it a}')
title ('Distribution function of parameter {\it a}')
legend ('Distribution', 'Mean {\it a}', 'Deviation')
 

Schrodinger's cat в сообщении #1567846 писал(а):

Отбрасываем конечные отрезки не вошедшие полностью и считаем среднее арифметическое длин всех оставшихся горизонтальных отрезков.

Я такое пробовал выше. Результат получается завышенным и мимо кассы совсем (хотя надуманная погрешность реально возможные значения покрывает):

$$a=0,0642\pm 0,0006$$

С помощью МНК коллега выше тоже посчитал. Результат получился заниженным и тоже мимо кассы. С погрешностью этого расчёта вообще непонятно что делать.

ilghiz в сообщении #1567771 писал(а):

$$a = 3371/53955 \simeq 0.0625$$