29

値がアプリオリに制限されていない、多数の数値セットの幾何平均を計算する必要があります。素朴な方法は

double geometric_mean(std::vector<double> const&data) // failure
{
  auto product = 1.0;
  for(auto x:data) product *= x;
  return std::pow(product,1.0/data.size());
}

ただし、これは累積のアンダーフローまたはオーバーフローのために失敗する可能性がありますproduct(注:long doubleこの問題を実際に回避するわけではありません)。したがって、次のオプションは対数を合計することです。

double geometric_mean(std::vector<double> const&data)
{
  auto sumlog = 0.0;
  for(auto x:data) sum_log += std::log(x);
  return std::exp(sum_log/data.size());
}

これは機能しますが、すべての要素を呼び出すstd::log()ため、遅くなる可能性があります。それを避けることはできますか？productたとえば、累算の指数と仮数 (と同等) を別々に追跡することによって?

4

7 に答える 7

12

「分割指数と仮数」ソリューション：

double geometric_mean(std::vector<double> const & data)
{
    double m = 1.0;
    long long ex = 0;
    double invN = 1.0 / data.size();

    for (double x : data)
    {
        int i;
        double f1 = std::frexp(x,&i);
        m*=f1;
        ex+=i;
    }

    return std::pow( std::numeric_limits<double>::radix,ex * invN) * std::pow(m,invN);
}

exがオーバーフローする可能性があることを懸念している場合は、 の代わりに double として定義し、すべてのステップでlong long乗算するinvNことができますが、このアプローチでは多くの精度が失われる可能性があります。

編集大きな入力の場合、計算をいくつかのバケットに分割できます。

double geometric_mean(std::vector<double> const & data)
{
    long long ex = 0;
    auto do_bucket = [&data,&ex](int first,int last) -> double
    {
        double ans = 1.0;
        for ( ;first != last;++first)
        {
            int i;
            ans *= std::frexp(data[first],&i);
            ex+=i;
        }
        return ans;
    };

    const int bucket_size = -std::log2( std::numeric_limits<double>::min() );
    std::size_t buckets = data.size() / bucket_size;

    double invN = 1.0 / data.size();
    double m = 1.0;

    for (std::size_t i = 0;i < buckets;++i)
        m *= std::pow( do_bucket(i * bucket_size,(i+1) * bucket_size),invN );

    m*= std::pow( do_bucket( buckets * bucket_size, data.size() ),invN );

    return std::pow( std::numeric_limits<double>::radix,ex * invN ) * m;
}
于 2013-11-14T15:54:54.223 に答える