CentralLimitTheorem.cpp

#include <iostream>
#include <vector>
#include <random>
#include <iomanip>
#include <algorithm>
#include <cmath>

int main() {
  const int probSize = 2560000; // 配列サイズ
  const int sampleNum = 1000000; // 実験サンプル数
  const int normalization = 1000; // 正規化量
  
  std::mt19937 random; // mersenne twister
  std::uniform_real_distribution<double> dist(-1.0, 1.0);
  std::vector<double> prob(probSize), gaussian(probSize);
  std::vector<double>::iterator p = prob.begin() + probSize/2;
  std::vector<double>::iterator g = gaussian.begin() + probSize/2;
  
  int n;
  std::cin >> n;
  
  // sampleNum個のサンプルについてSを求める
  for (int i = 0; i < sampleNum; i++) {
    long double S = 0;
    for (int j = 0; j < n; j++) {
      S += dist(random);
    }
    int index = S / sqrt(n / 3.0) * normalization;
    p[std::min(probSize/2-1, std::max(-probSize/2, index))] ++;
  }
  
  // 累積和による計算
  for (int i = 0; i < prob.size() - 1; i++) {
    prob[i+1] += prob[i];
  }
  
  // 正規分布の計算
  for (int i = 0; i < probSize - 1; i++) {
    double x = (double)(i - probSize / 2) / normalization;
    gaussian[i+1] = gaussian[i]
      + exp(-x*x/2) / sqrt(2.0 * acos(-1.0)) / normalization;
  }
  
  // 出力
  std::cout << std::fixed;
  std::setprecision(5);
  double errorAverage = 0;
  for (int i = - 2 * normalization; i <= 2 * normalization; i++) {
    std::cout << (double)i / normalization << "\t";
    std::cout << p[i] / sampleNum - g[i] << std::endl;
    errorAverage += pow(p[i] / sampleNum - g[i], 2.0) / (2 * normalization + 1);
  }
  std::cout << errorAverage << std::endl;
  return 0;
}