/*
 *   母平均の信頼区間の被覆確率をシミュレーションで近似する。
 *      母集団分布：自由度 p のカイ２乗分布
 *      母平均は p
 *   正規近似（中心極限定理）と、t 分布による近似とどっちが良いか知りたい。
 */

//#include <iostream.h>
#include <iostream>
using namespace std;

#include "random.c"
#include "Distributions.c"

// 標本数
#define  N  100

const long L = 10000;      // シミュレーションの繰り返し数
const int p = 5;          // 母平均（カイ２乗分布の自由度）
const double a = 0.05;    // 1 - 信頼係数
const unsigned long init = 13; // 擬似乱数の初期値

struct Stats
{
  int num;         // 標本数
  double xbar;   // 標本平均
  double sd;     // 不偏分散の平方根
};

// x : データ, n : データ数
void make_Stats(Stats &, double *x, int n);

// j[0] : 信頼区間の下限、j[1] : 信頼区間の上限、 z : 両側 100 a ％点
void interval(Stats &, double j[2], double z);

int main()
{
  init_genrand(init);
  long l;
  Stats t; 
  double x[N];
  double j_n[2], j_t[2];        // 信頼区間の下限と上限
                                // （_n は正規近似、_t は ｔ分布による近似）
  double count_n = 0, count_t = 0;  // 信頼区間が母平均を含んだ回数
  double p_n, p_t ;             // 被覆確率のシミュレーションによる推定値

  double z_n = n_inv(1 - a/2);
  double z_t  = sqrt( F_inv(1 - a, 1, N - 1) ); 
        // Ｔ〜自由度(n-1)のｔ分布のとき、T^2 は、自由度 1, (n-1) のＦ分布
  
  for(l = 0; l < L; l++){
    m_chi_rand(5, x, N);
    
    make_Stats(t, x, N);           // 標本平均、標本標準偏差の計算

    interval(t, j_n, z_n);         // 信頼区間の上限と下限の計算
    if(j_n[0] <= p && j_n[1] >= p)
      count_n += 1;                // 母平均を含んでいたら 1 加える

    interval(t, j_t, z_t);
    if(j_t[0] <= p && j_t[1] >= p)
      count_t += 1;
  }
  p_n = count_n/L;                 // 被覆確率のシミュレーション推定値
  p_t = count_t/L;

  cout << "標本数 = " << N
       << ", 自由度 = " << p
       << ", シミュレーションの繰り返し数 = " << L << endl;
  
  cout << "正規近似による被覆確率　：" << p_n
       << " ± " << sqrt(p_n*(1-p_n)/L)*z_n << endl;// ここは正規近似で、95%信頼区間

  cout << "ｔ分布近似による被覆確率：" << p_t
       << " ± " << sqrt(p_t*(1-p_t)/L)*z_n << endl; // ここも正規近似で、95%信頼区間
}

void make_Stats(Stats &t, double *x, int n)
{
  t.num = n;
  t.xbar = 0;
  t.sd = 0;
  int i;
  for(i = 0; i < n; i++){
    t.xbar += x[i];
    t.sd += x[i]*x[i];
  }
  t.xbar /= n;
  t.sd = sqrt( (t.sd - n*t.xbar*t.xbar)/(n-1) );
}

void interval(Stats &t, double j[2], double z)
{
  double h = z * t.sd / sqrt(t.num);
  j[0] = t.xbar - h;
  j[1] = t.xbar + h;
}