<numeric>数値演算アルゴリズムひとめぐり

inner_product：内積

　2つのコンテナ：X[],Y[]に対し、内積：X[0]*Y[0] + X[1]*Y[1] + ……を求めます。

template<class InputIterator1, InputIterator2, class T>
T inner_product(InputIterator1 first1, InputIterator1 last1,
                InputIterator2 first2, T init);

template<class InputIterator1, InputIterator2, class T,
         class BinaryOperation1, class BinaryOperation2>
T inner_product(InputIterator1 first1, InputIterator1 last1,
                InputIterator2 first2, T init,
                BinaryOperation1 binary_op1, BinaryOperation2 binary_op2);

　後者はaccumulateと同様、加算（operatio+）、乗算（operatior*）に相当する二項関数オブジェクト：binary_op1、binary_op2を第5、6引数に与えます。

　accumulateとinner_productのサンプルとして最小二乗法による直線近似をやってみました。

　データの組がいくつか(X[N], Y[N])あり、Y[i]とf(X[i])の差の総和が最も小さくなるような直線f(x) = a*x + bのa, bを最小二乗法で求めます。（コ難しい偏微分方程式は端折って）a, bは以下の式から求まります:

fig-1

　上式に現れるΣx, Σyはaccumulate、Σx*x, Σx*y はinner_productを使って求めることができます。さらに、相関の強さを示す相関係数：r は、

fig-2

　平均値mx, myはそれぞれΣx/N, Σy/N ですから、X[], Y[]の各要素からmx, myを引いておけば上式のS(xy), S(xx), S(yy)は、それぞれ X[]とY[], X[]とX[], Y[]とY[]のinner_productです。

list-8　最小二乗法

#include <numeric>
#include <vector>
#include <random>
#include <iostream>
#include <algorithm>

using namespace std;

int main() {
  const size_t N = 100;

  vector<float> X;
  vector<float> Y;

  // データを作る: Y[i] = 5.0 * X[i] + 10 （+ 乱数）
  mt19937 gen;
  normal_distribution<float> dist(0.0f, 1.0f);
  for ( size_t i = 0; i < N; ++i ) {
    float xi = (float)i / N;
    X.emplace_back(xi);
    float yi = 5.0f * xi + 10.0f + dist(gen);
    Y.emplace_back(yi);
  }

  // X[],Y[]の総和 Σx, Σy
  float sx = accumulate(begin(X), end(X), 0.0f);
  float sy = accumulate(begin(Y), end(Y), 0.0f);
  // X[],Y[]の内積 Σxy
  float sxy = inner_product(begin(X), end(X), begin(Y), 0.0f);
  // X[]の二乗和 Σx^2
  float sxx = inner_product(begin(X), end(X), begin(X), 0.0f);

  float a = (N*sxy  - sx*sy ) / (N*sxx - sx*sx);
  float b = (sxx*sy - sxy*sx) / (N*sxx - sx*sx);

  // X[],Y[]の平均
  float mx = sx / N;
  float my = sy / N;

  // X[],Y[]の各要素から平均を引く: (x, y) → (x-mx, y-my)
  transform(begin(X), end(X), begin(X), [=](float x) { return x - mx;});  
  transform(begin(Y), end(Y), begin(Y), [=](float y) { return y - my;});  
  // X[]とY[]の内積 Σ(x-mx)(y-my)
  sxy = inner_product(begin(X), end(X), begin(Y), 0.0f);
  // X[]とX[]の内積 Σ(x-mx)^2
  sxx = inner_product(begin(X), end(X), begin(X), 0.0f);
  // Y[]とY[]の内積 Σ(y-my)^2
  float syy;
  syy = inner_product(begin(Y), end(Y), begin(Y), 0.0f);

  // 相関係数 r
  float r = sxy / (sqrtf(sxx)*sqrtf(syy));

  cout << "a = " << a << " b = " << b << " r = " << r << endl;

}