/* 課題１解答例 */

/* ［数値処理］→［自動的に数値で出力］*/
if numer#false then numer:false else numer:true;

/*浮動小数点の桁数=7に指定 */
fpprec : 7;

/* 浮動小数点表示桁数を7に指定*/
fpprintprec:7;

/* (I) */
Center(n):=(
load(distrib),
/* 正規乱数を発生させるためにdistribパッケージを読み込む */
StdSampleMean: makelist( mean (random_normal(0,1,n))/(1/sqrt(n)), 10000),
/* n個の標準正規乱数の標準化平均，10000個のリスト */
load(descriptive),
/* ヒストグラムを描くためにdescriptiveパッケージを読み込む*/
load(draw), /* 描画パッケージ読み込み */
draw2d(
histogram_description(StdSampleMean, frequency=relative),
/* StdSampleMeanのヒストグラム，相対度数表示*/
explicit(pdf_normal(x,0,1),x, -4, 4),
/* 標準正規分布確率密度関数グラフ */
yrange=[0,0.4], xrange=[-4,4]
)
);

kill(all);

/* (II) */
load(distrib);
data0 : random_normal(0,1,20);
data1 : random_normal(1,1,20);
load(stats) /*推測統計パッケージstats読み込み*/;

/* (1) */
test_mean(data0);
print("t値=",
t : mean(data0)/std1(data0)*sqrt(20));
print("p値=", 
2 * (1-cdf_student_t(abs(t),19)));
/* 信頼区間 */
print("信頼区間=",
[ mean(data0)-quantile_student_t(0.975,19)*std1(data0)/sqrt(20),
  mean(data0)-quantile_student_t(0.025,19)*std1(data0)/sqrt(20) ]);
print("平均推定値=",
mean(data0));

/* (2) */
test_mean(data1, conflevel=0.99);
print("t値=",
t : mean(data1)/std1(data1)*sqrt(20) );
print("p値=", 
2 * (1-cdf_student_t(abs(t),19)) );
print("信頼区間=",
[ mean(data1)-quantile_student_t(0.995,19)*std1(data1)/sqrt(20),
  mean(data1)-quantile_student_t(0.005,19)*std1(data1)/sqrt(20) ] );
print("平均推定値=",
mean(data1));

/* (3) */
test_mean(data1,conflevel=0.99,alternative=less);
print("t値=",
t : mean(data1)/std1(data1)*sqrt(20));
print("p値=", 
cdf_student_t(abs(t),19));
print("信頼上限=",
mean(data1)-quantile_student_t(0.01,19)*std1(data1)/sqrt(20) );
print("平均推定値=",
mean(data1));

/* (4) */
test_means_difference(data0,data1,varequal=true); 
s : sqrt((19*var1(data0)+19*var1(data1))/(20+20-2));
print("t値=",
t : (mean(data0)-mean(data1))/(s*sqrt(1/20+1/20)));
print("p値=", 
2 * (1-cdf_student_t(abs(t),20+20-2)));
print("信頼区間=",
[ mean(data0)-mean(data1)-quantile_student_t(0.975,20+20-2)*s*sqrt(1/20+1/20),
  mean(data0)-mean(data1)-quantile_student_t(0.025,20+20-2)*s*sqrt(1/20+1/20)]
);
print("平均差推定値=",
mean(data0)-mean(data1));

kill(all);