Преобразование однопроходной формулы StandardDev из C в Python

avatar
Lance Gibson
8 августа 2021 в 19:45
69
1
-1

Я пытаюсь преобразовать следующий однопроходный код стандартного отклонения в C в Python:

double std_dev(double a[], int n) {
    if(n == 0)
        return 0.0;
    int i = 0;
    double meanSum = a[0];
    double stdDevSum = 0.0;
    for(i = 1; i < n; ++i) {
        double stepSum = a[i] - meanSum;
        double stepMean = ((i - 1) * stepSum) / i;
        meanSum += stepMean;
        stdDevSum += stepMean * stepSum;
    }
    // for poulation variance: return sqrt(stdDevSum / n);
    return sqrt(stdDevSum / (n));

Вот что у меня есть на Python:

def std_dev(a,n):
    if n == 0:
        return 0.0
    i = 0
    meanSum = float(a[0])
    stdDevSum = float(0.0)
    for i in range(1,n,1):
        stepSum = float(float(a[i]) - meanSum)
        stepMean = float(((i - 1)*stepSum)/i)
        meanSum += stepMean
        stdDevSum += stepMean*stepSum
        print(stdDevSum)
    value = float(sqrt(stdDevSum/(n)))
    print(value)

Однако я не получаю правильного результата для стандартного отклонения совокупности. Например, программа возвращает стандартное отклонение набора [10,20,500,40,50] как 175,33, тогда как онлайн-калькулятор или ручной расчет возвращает 188,53. Как объяснить разницу?

Спасибо!

Алгоритм C Источник: https://www.strchr.com/standard_deviation_in_one_pass

Источник
Eric Postpischil
8 августа 2021 в 19:59
0

Каково наименьшее количество чисел, для которых вы обнаружили, что это не возвращает правильных результатов? Кто они такие? Прежде чем опубликовать здесь, вы распечатали значения каждой переменной и увидели, чем они отличаются между реализациями C и Python? На какой итерации происходит первое отклонение и в какой переменной? Отредактируйте вопрос, чтобы предоставить минимальный воспроизводимый пример.

Eric Postpischil
8 августа 2021 в 20:01
0

Почему код Python имеет sqrt(stdDevSum/(n), а код C имеет sqrt(stdDevSum / (n - 1)?

Adrian Mole
8 августа 2021 в 20:02
0

@Eric - я тоже спрашивал об этом, но потом увидел комментарий в конце кода C.

Eric Postpischil
8 августа 2021 в 20:03
0

@AdrianMole: Тем не менее, код должен быть аналогичным. Код C с комментарием о том, что может быть получен другой результат, не эквивалентен коду Python, который дает другой результат.

Adrian Mole
8 августа 2021 в 20:04
0

@Eric Эрик, я согласен - ОП должен отредактировать, чтобы не было версии «образцовой сигмы», а только кода «популяционной сигмы».

Adrian Mole
9 августа 2021 в 06:14
0

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

Ответы (1)

avatar
Adrian Mole
8 августа 2021 в 20:02
1

Синтаксис вашего выражения range() неверен. Третий аргумент должен быть действительным приращением (т. е. 1), а не выражением, включающим индекс цикла. [ Ссылка ]

Кроме того, вы должны использовать math.sqrt() и вам необходимо import math:

import math 

def std_dev(a,n):
    if n == 0:
        return 0.0
#   i = 0 # This line is unnecessary
    meanSum = float(a[0])
    stdDevSum = float(0.0)
    for i in range(1,n,1): # Note the third argument!
#   for i in range(1, n):  # Alternative version - the default for "step" is 1
        stepSum = float(float(a[i]) - meanSum)
        stepMean = float(((i - 1)*stepSum)/i)
        meanSum += stepMean
        stdDevSum += stepMean*stepSum
        print(stdDevSum)
    value = float(math.sqrt(stdDevSum/(n)))
    print(value)

test = [1.  2.  3.  2.  1.]
std_dev(test, 5)

Вывод:

0.0
2.0
2.0
2.75
0.7416198487095663

Используя те же тестовые данные с вашим кодом C (или его версией, использующей делитель n, а не (n - 1)), следующим образом:

int main()
{
    double test[] = { 1.  2.  3.  2.  1. };
    double answer = std_dev(test, 5);
    printf("%lg\n", answer);
    return 0;
}

дает 0.74162 в качестве вывода.

Lance Gibson
9 августа 2021 в 04:30
0

Извините за путаницу в отношении выборки и населения. Нам нужна SD популяция, а не выборка. Кроме того, я исправил проблему с диапазоном(). Верный вопрос, почему эта версия не возвращает полностью точную SD. В вашем примере [1,2,3,2,1] стандартное отклонение было достаточно маленьким, чтобы результат был разумным. Однако возьмите набор [10,20,500,40,50], и вычисленное программой SD станет сильно отличаться от того, что вы вычислили бы вручную или с помощью онлайн-калькулятора. Как мы можем объяснить эту неточность в программе?

Adrian Mole
9 августа 2021 в 05:45
0

Хм. Но и код C, и код Python дают один и тот же ответ (175.339).

Lance Gibson
9 августа 2021 в 05:52
0

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