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

  • Великая Октябрьская резня бензопилой

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

  • Очередная несуразность в единицах измерения

    Когда-то я написал программу PhysUnitCalc - калькулятор, умеющий работать с размерностями. Мне казалось, что я наступил уже на все грабли, которые…

  • Big Data, чтоб их... (3)

    "В предыдущих сериях": мой прибор выдаёт 6 значений: 3 координаты и 3 угла, т.е все 6 степеней свободы твёрдого тела. Причём ошибки измерения этих 6…

  • 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

  • Великая Октябрьская резня бензопилой

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

  • Очередная несуразность в единицах измерения

    Когда-то я написал программу PhysUnitCalc - калькулятор, умеющий работать с размерностями. Мне казалось, что я наступил уже на все грабли, которые…

  • Big Data, чтоб их... (3)

    "В предыдущих сериях": мой прибор выдаёт 6 значений: 3 координаты и 3 угла, т.е все 6 степеней свободы твёрдого тела. Причём ошибки измерения этих 6…