шаблоны — Взвешенная дисперсия и взвешенное стандартное отклонение в переполнении стека

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

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 Docs, которую я настроил для изучения проблемы. К сожалению, мои ответы нигде не близки к NIST pdf. Я думаю, что проблема находится в шаге непредвзятости, но я не вижу, как это исправить.

https://docs.google.com/spreadsheet/ccc?key=0ApzPh5nRin0ldGNNYjhCUTlWTks2TGJrZW4wQUcyZnc&УСП = обмен

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

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

float mean(uint16_t* x, uint16_t n) {
uint16_t sum_xi = 0;
int i;
for (i = 0; i < n; i++) {
sum_xi += x[i];
}
return (float) sum_xi / n;
}

/**
* http://www.itl.nist.gov/div898/software/dataplot/refman2/ch2/weigmean.pdf
*/
float weighted_mean(uint16_t* x, uint16_t* w, uint16_t n) {
int sum_wixi = 0;
int sum_wi = 0;
int i;
for (i = 0; i < n; i++) {
sum_wixi += w[i] * x[i];
sum_wi += w[i];
}
return (float) sum_wixi / (float) sum_wi;
}

float variance(uint16_t* x, uint16_t n) {
float mean_x = mean(x, n);
float dist, dist2;
float sum_dist2 = 0;

int i;
for (i = 0; i < n; i++) {
dist = x[i] - mean_x;
dist2 = dist * dist;
sum_dist2 += dist2;
}

return sum_dist2 / (n - 1);
}

/**
* http://www.itl.nist.gov/div898/software/dataplot/refman2/ch2/weighvar.pdf
*/
float weighted_variance(uint16_t* x, uint16_t* w, uint16_t n) {
float xw = weighted_mean(x, w, n);
float dist, dist2;
float sum_wi_times_dist2 = 0;
int sum_wi = 0;
int n_prime = 0;

int i;
for (i = 0; i < n; i++) {
dist = x[i] - xw;
dist2 = dist * dist;
sum_wi_times_dist2 += w[i] * dist2;
sum_wi += w[i];

if (w[i] > 0)
n_prime++;
}

if (n_prime > 1) {
return sum_wi_times_dist2 / ((float) ((n_prime - 1) * sum_wi) / n_prime);
} else {
return 0.0f;
}
}

/**
* http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Weighted_incremental_algorithm
*/
float weighted_incremental_variance(uint16_t* x, uint16_t* w, uint16_t n) {
uint16_t sumweight = 0;
float mean = 0;
float M2 = 0;
int n_prime = 0;

uint16_t temp;
float delta;
float R;

int i;
for (i = 0; i < n; i++) {
if (w[i] == 0)
continue;

temp = w[i] + sumweight;
delta = x[i] - mean;
R = delta * w[i] / temp;
mean += R;
M2 += sumweight * delta * R;
sumweight = temp;

n_prime++;
}

if (n_prime > 1) {
float variance_n = M2 / sumweight;
return variance_n * n_prime / (n_prime - 1);
} else {
return 0.0f;
}
}

void main(void) {
uint16_t n = 9;
uint16_t x[] = { 2, 3, 5, 7, 11, 13, 17, 19, 23 };
uint16_t w[] = { 1, 1, 0, 0,  4,  1,  2,  1,  0 };

printf("%f\n", weighted_variance(x, w, n)); /* outputs: 33.900002 */
printf("%f\n", weighted_incremental_variance(x, w, n)); /* outputs: 33.900005 */
}

5.82

template <class T>
inline float    VarianceWeighted( T samples[], T weights[], int count )
{
...
varianceWeighted *= ((float) numNonzero)/((float) count*(numNonzero - 1.0f)*sumWeights);    // this should be right but isn't?!
...
}

template <class T>
inline float    VarianceWeighted( T samples[], T weights[], int count )
{
...
varianceWeighted *= ((float) numNonzero)/((float) (numNonzero - 1.0f)*sumWeights);  // removed count term
...
}