nabbla (nabbla1) wrote,
nabbla
nabbla1

Categories:

Как не надо считать среднеквадратичное отклонение

Если у нас есть набор значений xn (n=1..N), и мы хотим посчитать для них среднеквадратичное отклонение, есть две классические формулы:

std1.png

std2.png

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

Кнут предостерегает - делать этого нельзя, можно обжечься, ведь при большом среднем (по сравнению со среднеквадратичным) у нас будут вычитаться два близких значения, можно ненароком и отрицательную дисперсию схлопотать! Это опасение я и студентам передал в прошлом году, после чего спокойненько продолжал пользоваться второй формулой и получал нормальные результаты, дескать, при 80-битных числах с плавающей точкой и маленьких выборках можно не заморачиваться. (Да, еще есть формула оценки ср. кв по выборке, отличающаяся от верхей формулы знаменателем N-1 вместо N, но речь сейчас не о том)

Ага, ЩАЗ. Вчера-таки оно выстрелило - извлечение корня из отрицательного числа, чего в принципе быть не должно. Ну, подумал - да, не надо было столько выборок, таки накопилась ошибка. Ладно, лень разбираться, пока поставлю 10 выборок, грубо, но жить можно. Не тут-то было, на 10 (!) выборках оно ломается. Вот уж чего не ожидал, все мои кватернионные алгоритмы на 32-битном DSP и то за 40 минут сильно ошибки не накапливают, норма остается в норме, да что же это такое-то!?

Оказалось, "убить" 2-ю формулу можно элементарно, повторяемость 50%. Вот как это делается.


Даем ей 10 одинаковых значений - у нас должна получиться нулевая дисперсия, но из-за ошибок в вычислениях она может получиться как немножко положительной, так и немножко отрицательной. Конечно, если мы какие-то простенькие значения возьмем, вроде 1 или 2, то при возведении в квадрат ошибки вообще не возникнет, все в порядке. А вот на 0.11 уже получим отрицательную дисперсию, поскольку это в двоичной форме бесконечная дробь, и такая, что дисперсия у нас вылезет отрицательная. Совершенно без разницы, какая у нас точность - 32 бита, 64, 80 или 1000 бит, такая ошибка все равно может проявиться.

Который раз так получается - чтобы свести с ума какой-нибудь метод, нужно дать ему абсурдно простое задание. Простой QuickSort (без случайного выбора) тратит N2 операций на сортировку уже отсортированного массива, мой хитрючий кватернионный алгоритм проглючил только когда ему сказали двигаться только по углу места от 0 до 110, и назад. Симуляторы электронных цепей неплохо рассчитывают здоровенные цифровые и аналоговые схемы, но на обычном диодном мостике довольно любопытно "залипают", об этом как-нибудь расскажу.

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

И ссылки на первоисточники,
раз: http://zach.in.tu-clausthal.de/teaching/info_literatur/Welford.pdf
два: https://yadi.sk/d/TEnFbsbYcrebS
(Дональд Кнут - Искусство программирования, том 2, нас интересуют страницы 258-259)
Tags: Монте-Карло для чайников, математика, моделирование
Subscribe

Recent Posts from This Journal

  • Нахождение двух самых отдалённых точек

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

  • Слишком общительный счётчик

    Вчера я чуть поторопился отсинтезировать проект,параметры не поменял: RomWidth = 8 вместо 7, RamWidth = 9 вместо 8, и ещё EnableByteAccess=1, чтобы…

  • Балансируем конвейер QuatCore

    В пятницу у нас всё замечательно сработало на симуляции, первые 16 миллисекунд полёт нормальный. А вот прошить весь проект на ПЛИС и попробовать "в…

  • Post a new comment

    Error

    Anonymous comments are disabled in this journal

    default userpic

    Your IP address will be recorded 

  • 12 comments

Recent Posts from This Journal

  • Нахождение двух самых отдалённых точек

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

  • Слишком общительный счётчик

    Вчера я чуть поторопился отсинтезировать проект,параметры не поменял: RomWidth = 8 вместо 7, RamWidth = 9 вместо 8, и ещё EnableByteAccess=1, чтобы…

  • Балансируем конвейер QuatCore

    В пятницу у нас всё замечательно сработало на симуляции, первые 16 миллисекунд полёт нормальный. А вот прошить весь проект на ПЛИС и попробовать "в…