密度トレイスを計算するための関数を書いてみた

とりあえずこんな感じ(どんな感じだ)で密度トレイスを出力できるようにした。

/Users/yasuhisa/cpp% ./a.out
2.72898e-07
0.000121768
0.00827342
0.100135
0.290664
0.319885
0.203549
0.064411
0.0112637
0.00164454

正規乱数を発生させて、そこから計算させてます。データ一つに対して、インスタンス一つ生成しているのがどうなのかというところだけど、まあ最初だからとりあえずできることが重要ということにしておく。boostとかで適当にメルセンヌツイスターやっていたりなど。テンプレートしてみようかと思ったけど、意味分からずコピペではきつそうだったからやめておいた。

#include <functional>
#include <algorithm>
#include <vector>
#include <iterator>
#include <iostream>
#include <boost/foreach.hpp>
#include <cmath>
#include <boost/random.hpp>

using namespace std;
using namespace boost;

class Kernel
{
  double X;
public:
  Kernel(double x);
  double density(double x, double h);
  void show();
  operator double () { return X; }
};

Kernel::Kernel(double x)
{
  X = x;
}

double Kernel::density(double x, double h)
{
  const static double PI = 4*atan(1.0);
  return exp( - pow((x-X)/h,2) / 2) / sqrt(2*PI);
}

void Kernel::show()
{
  cout << X << endl;
}

double sum_of_kernel_density(vector<Kernel> vec, double x, double h);
const uint32_t SEED = 123;
int main()
{
  mt19937 gen(SEED);
  normal_distribution<> dst(5.0, 1.0);

  variate_generator<
  mt19937, normal_distribution<>
	> r_norm(gen, dst);
  
  vector<Kernel> nums;

  for( int i = 0; i < 100; ++i )
	{
	  Kernel k(r_norm());
	  nums.push_back( k );
	}

  double h;
  h = 0.7;
  for(int i=0;i<10;++i)
	{
	  cout << sum_of_kernel_density(nums, i, h) << endl;
	}

  return 0;
}

double sum_of_kernel_density(vector<Kernel> vec, double x, double h)
{
  vector<Kernel>::iterator it;
  it = vec.begin();
  double sum;
  sum = 0.0;
  while( it != vec.end() )
	{
	  sum = sum + (*it).density(x, h);
	  ++it;
	}
  return sum / (vec.size() * h) ;
}