LINUX.ORG.RU
ФорумTalks

[арифметика] машинные вычисления. криворукость или погрешность?


0

0

возник вопрос по подсчету СКО случайной величины. существуют 2 формулы, которые аналитически равносильны.

S=sqrt(M(x^2)-M(x)^2)

и

S=M(x^2-M(x)^2)

M-матожидание.

интересует вопрос, если считать все это на x86 в double то в каком случае погрешность будет больше? скорость не важна. вопрос возник в связи с тем что считая по первой сормуле СКО иногда вылазит за размах выборки в обе стороны. насколько такое поведение для СВ может быть нормальным?


у дабла ноль вроде где-то около 10е-16. т.е. альфа * (1 + ноль) = альфа.

вполне возможно, что где-то, пока идёт счёт М - достигается такое условие.

сколько у вас точек у выборки? какой разброс величин?

gunja
()
Ответ на: комментарий от gunja

в общем, там считается статистика по выборке значений в количестве от ~500 до ~80000 в зависимости от условий.

разброс - от ~200 до ~500000 - максимальный найденный. знаю что данные странные. при попытке считать то же самое в целых дисперсия оказывалась иногда отрицательной...

iero
() автор топика
Ответ на: комментарий от gunja

еще подумалось, спасет ли в данном случае какой-либо вариант округления последних цифр (все равно это время в микросекундах и на выходе такая точность не требуется)?

iero
() автор топика
Ответ на: комментарий от iero

округление не спасёт, пока всё не разделить на одну и ту же величину.
вообще, лучше смотреть на ваш код, на то, как именно вы это делаете. 
требования к времени я не знаю ваши, но самое простое решение, которое приходит в голову:

1) сортируем массив по возрастанию.
2) вычисляем число элементов в массиве
3) на каждом шаге при расчёте матожидания (средняя на сколько я помню) - к результату прибавляем величина / кол-во элементов.

да, операций станет в разы больше. точность должна улететь ещё сильнее.
но где-то же ошибка. код покажите этого участка.

gunja
()
Ответ на: комментарий от iero

>>еще подумалось, спасет ли в данном случае какой-либо вариант округления последних цифр?

От чего спасет? Хотите увеличить точность округлением?

gkrellm
()
Ответ на: комментарий от gunja

        QVector<double> Averages(numbers);
        QVector<double> Deviats(numbers);
        QVector<int> Counts(numbers);
        Deviats.fill(0);
        Averages.fill(0);
        Maxs.fill(0);
        Counts.fill(0);
        Mins.fill(10000000);
        int total=qwmodel.rowCount();
        for (int i=0;i<total;i++)
        {
            currentindex=qwmodel.record(i).value("number").toInt(NULL);
            currentval=qwmodel.record(i).value("rtt").toInt(NULL);
            Averages[currentindex]+=currentval;
            Counts[currentindex]++;
            Deviats[currentindex]+=currentval*currentval;

        }
        for(int j=0;j<numbers;j++)
        {
            Averages[j]/=Counts[j];
            Deviats[j]=sqrt(Deviats[j]/Counts[j]-Averages[j]*Averages[j]);
        }
знаю что получается костыльная группировка по полю "number", другого способа считать СКО не нашел.

iero
() автор топика
Ответ на: комментарий от iero

то ли глаза у меня уставшие. то ли всё верно. что рекомендовал бы сделать - скобками обозначить явный приоритет операций (лично я его просто не помню, возможно там загвоздка). проверить, что numbers всегда меньше возвращаемых currentindex. и ещё, явные приведения типа сделайте. не поможет, но кто знает.

во втором цикле делите всегда на Counts[j], которые из кода не всегда ненулевые.

больше идеев немае

gunja
()

Откройте для себя формулу Kahan

summation of ten million fractions 

    int main(void){
        int i;
        float sum;
        float c,y,t;
  
        sum = 0.0;
        for(i=1; i<=10000000; i++){
            sum = sum + 1.0/((float)i);
        }
        printf("decreasing order: %f\n",sum);
  
        sum = 0.0;
        for(i=10000000; i>0; i--){
            sum = sum + 1.0/((float)i);
        }
        printf("increasing order: %f\n",sum);
  
        /* sum formula suggested by W. Kahan */
        sum = 1.0/((float)1);
        c = 0.0;
        for(i=2; i<=10000000; i++){
            y = 1.0/((float)i) - c;
            t = sum + y;
            c = (t - sum) - y;
            sum = t;
        }
        printf("kahan summation1: %f\n",sum);
  
        sum = 1.0/((float)10000000);
        c = 0.0;
        for(i=9999999; i>0; i--){
            y = 1.0/((float)i) - c;
            t = sum + y;
            c = (t - sum) - y;
            sum = t;
        }
        printf("kahan summation2: %f\n",sum);
    }

    /* output:

    float:

    decreasing order: 15.403683
    increasing order: 16.686031
    kahan summation1: 16.695311
    kahan summation2: 16.695311

    double:

    decreasing order: 16.695311
    increasing order: 16.695311
    kahan summation1: 16.695311
    kahan summation2: 16.695311

    */

sign
()
Ответ на: комментарий от iero

приведение попробовать 

- Averages[currentindex] +=currentval;
+ (Averages[currentindex]) += static_cast<double> (currentval); 

- Averages[j]/=Counts[j];
- Deviats[j]=sqrt(Deviats[j]/Counts[j]-Averages[j]*Averages[j]); 

+ (Averages[j]) /= static_cast<double>(Counts[j]);
+ Deviats[j]=sqrt(Deviats[j]/static_cast<double>(Counts[j]) - Averages[j]*Averages[j]); 

и да... я сейчас понимаю, что я - быдлокодер

gunja
()
Ответ на: комментарий от iero

в релизе - должно. в дебаг сборке у меня там NAN записан. какое поведение может добавить Qt - не знаю.

короче не обязано оно сегфолтить... уже не обязано

gunja
()
Ответ на: комментарий от iero

> при делении double на нуль

проверяю правда на виндах:

D:\>type test.cpp
#include <stdio.h>

int main() {
        long double a;
        int c;
        a = 1.;
        c=0;
        a /=c;
        printf("%g\n", a);
        return 0;
};


D:\>cl test.cpp
Microsoft (R) 32-bit C/C++ Optimizing Compiler Version 14.00.50727.762 for 80x86
Copyright (C) Microsoft Corporation.  All rights reserved.

test.cpp
Microsoft (R) Incremental Linker Version 8.00.50727.762
Copyright (C) Microsoft Corporation.  All rights reserved.

/out:test.exe
test.obj

D:\>test.exe
1.#INF

т.е. валиться-то и не обязаны. но это cpp

gunja
()
Ответ на: комментарий от gunja

там у куте еще и toDouble оказалось. но не дал ровным счетом никаког эффекта. буду прикручивать предложенные суммы.

iero
() автор топика
Ответ на: комментарий от iero

так... а корень в первой формуле откуда взялся? первая и вторая формулы по размерностям не сходятся. так и должно разве быть?

gunja
()
Ответ на: комментарий от iero

Deviats[currentindex]+=currentval*currentval;
...
Deviats[j]=sqrt(Deviats[j]/Counts[j]-Averages[j]*Averages[j]);

вроде сходятся?
это же среднеквадр. отклонение - корень из дисперсии.

iero
() автор топика
Ответ на: комментарий от iero

S=sqrt(M(x^2)-M(x)^2)

и

S=M(x^2-M(x)^2)

M-матожидание.

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

gunja
()

и я бы сказал, что формулы вообще какие-то отсебятины.

что видно в вики: http://ru.wikipedia.org/wiki/Дисперсия_случайной_величины

D[x] = M[ (x - M[x]) ^2 ] == M[X^2] - (M[x])^2

т.е. из первой формулы брать корень - не нужно. во второй формуле вы очень некрасиво раскрыли возведение в квадрат.

gunja
()
Ответ на: комментарий от gunja

>>т.е. из первой формулы брать корень - не нужно

Первая правильна - sqrt из дисперсии.

А вот вторая - вы правы - не сходится даже по размерности. Там нужно корень еще дописать.

gkrellm
()
Ответ на: комментарий от gkrellm

во второй, как я понял, в квадрат надо возводить отклонение от центра, но никак не величину в квадрат минус квадрат центра

gunja
()
Ответ на: комментарий от gunja

>>в квадрат надо возводить отклонение от центра, но никак не величину в квадрат минус квадрат центра

Чего гадать - на бумаге из первой выводится в две строки - там лишь корня не хватает во второй.

Неясно другое - может ли СКО выйти за выборку. Чего-то мне эксель даже считал вроде что может - хотя здравый смысл и вид формулы говорит об обратном.

gkrellm
()
Ответ на: комментарий от gkrellm

D[X] = M [ (x - M[x])^2 ] = M [x^2 - 2 x M[x] + (M[x])^2 ].
M[ M[x]] = M[x]
D[x] = M[x^2] - 2 M[x] * M[x] + M[x] * M[x] = M[x^2] - (M[x])^2

это вывод первой формулы. снова обращаю внимание, что
(x - M[x]) ^2 != x^2 - M[x]^2

т.е. вторая формула не верная. а вот выходить-то вроде не должно. можно набор данных, на котором выходило за выборку?

gunja
()
Ответ на: комментарий от gunja

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

iero
() автор топика
Вы не можете добавлять комментарии в эту тему. Тема перемещена в архив.