Взвешенная дисперсия и взвешенное стандартное отклонение в C++

Привет, я пытаюсь рассчитать взвешенную дисперсию и взвешенное стандартное отклонение ряда целых чисел или чисел с плавающей запятой. Я нашел эти ссылки:

http://math.tutorvista.com/statistics/standard-deviation.html#weighted-standard-deviation

http://www.itl.nist.gov/div898/software/dataplot/refman2/ch2/weightsd.pdf (предупреждающий pdf)

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

template <class T>
inline float    Mean( T samples[], int count )
{
    float   mean = 0.0f;

    if( count >= 1 )
    {
        for( int i = 0; i < count; i++ )
            mean += samples[i];

        mean /= (float) count;
    }

    return mean;
}

template <class T>
inline float    Variance( T samples[], int count )
{
    float   variance = 0.0f;

    if( count > 1 )
    {
        float   mean = 0.0f;

        for( int i = 0; i < count; i++ )
            mean += samples[i];

        mean /= (float) count;

        for( int i = 0; i < count; i++ )
        {
            float   sum = (float) samples[i] - mean;

            variance += sum*sum;
        }

        variance /= (float) count - 1.0f;
    }

    return variance;
}

template <class T>
inline float    StdDev( T samples[], int count )
{
    return sqrtf( Variance( samples, count ) );
}

template <class T>
inline float    VarianceWeighted( T samples[], T weights[], int count )
{
    float   varianceWeighted = 0.0f;

    if( count > 1 )
    {
        float   sumWeights = 0.0f, meanWeighted = 0.0f;
        int     numNonzero = 0;

        for( int i = 0; i < count; i++ )
        {
            meanWeighted += samples[i]*weights[i];
            sumWeights += weights[i];

            if( ((float) weights[i]) != 0.0f ) numNonzero++;
        }

        if( sumWeights != 0.0f && numNonzero > 1 )
        {
            meanWeighted /= sumWeights;

            for( int i = 0; i < count; i++ )
            {
                float   sum = samples[i] - meanWeighted;

                varianceWeighted += weights[i]*sum*sum;
            }

            varianceWeighted *= ((float) numNonzero)/((float) count*(numNonzero - 1.0f)*sumWeights);    // this should be right but isn't?!
        }
    }

    return varianceWeighted;
}

template <class T>
inline float    StdDevWeighted( T samples[], T weights[], int count )
{
    return sqrtf( VarianceWeighted( samples, weights, count ) );
}

Прецедент:

int     samples[] = { 2, 3, 5, 7, 11, 13, 17, 19, 23 };

printf( "%.2f\n", StdDev( samples, 9 ) );

int     weights[] = { 1, 1, 0, 0, 4, 1, 2, 1, 0 };

printf( "%.2f\n", StdDevWeighted( samples, weights, 9 ) );

Результат:

7.46
1.94

Должно быть:

7.46
5.82

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

http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Weighted_incremental_algorithm

template <class T>
inline float    VarianceWeighted( T samples[], T weights[], int count )
{
    float   varianceWeighted = 0.0f;

    if( count > 1 )
    {
        float   sumWeights = 0.0f, meanWeighted = 0.0f, m2 = 0.0f;

        for( int i = 0; i < count; i++ )
        {
            float   temp = weights[i] + sumWeights,
                    delta = samples[i] - meanWeighted,
                    r = delta*weights[i]/temp;

            meanWeighted += r;
            m2 += sumWeights*delta*r;   // Alternatively, m2 += weights[i] * delta * (samples[i]−meanWeighted)
            sumWeights = temp;
        }

        varianceWeighted = (m2/sumWeights)*((float) count/(count - 1));
    }

    return varianceWeighted;
}

Результат:

7.46
5.64

Я также пытался посмотреть на boost и esutil, но они не очень помогли:

http://www.boost.org/doc/libs/1_48_0/boost/accumulators/statistics/weighted_variance.hpp http://esutil.googlecode.com/svn-history/r269/trunk/esutil/stat/util.py

Мне не нужна целая библиотека статистики, и, что более важно, я хочу понять реализацию.

Может ли кто-нибудь опубликовать функции для их правильного расчета?

Бонусные баллы, если ваши функции могут сделать это за один проход.

Кроме того, кто-нибудь знает, дает ли взвешенная дисперсия тот же результат, что и обычная дисперсия с повторяющимися значениями? Например, будет ли дисперсия выборок [] = {1, 2, 3, 3} такой же, как взвешенная дисперсия выборок [] = {1, 2, 3}, весов [] = {1, 1, 2} ?

Обновление: вот таблица документов Google, которую я создал для изучения проблемы. К сожалению, мои ответы далеки от PDF-файла NIST. Я думаю, что проблема в шаге unbias, но я не вижу, как это исправить.

https://docs.google.com/spreadsheet/ccc?key=0ApzPh5nRin0ldGNNYjhCUTlWTks2TGJrZW4wQUcyZnc&usp=sharing

Результатом является взвешенная дисперсия 3,77, которая является квадратом взвешенного стандартного отклонения 1,94, полученного в моем коде на С++.

Я пытаюсь установить октаву в своей установке Mac OS X, чтобы я мог запускать их функцию var() с весами, но установка ее с помощью brew занимает целую вечность. Сейчас я увлечен бритьем яков.