mysample.c 114514

作者：小鳥遊時宿

/***********************************************************************

mysample.c：

情報科学基礎実験I

[2] 正規母集団と統計的推測

課題プログラム（標本抽出実験プログラム）

概要：

母集団データから，与えられた標本数，繰り返し回数で

標本の無作為抽出を行い，統計量と信頼区間の計算を行う．

コマンド： ./mysample population_file sample_size repeat_count

./mysample population_file < urand_data

第1引数(population_file)：母集団データ(一変量データ)ファイル名

第2引数(sample_size )：標本サイズ

第3引数(repeat_count )：標本抽出の繰り返し回数

母集団データ(一変量データ)のファイル形式：

1行目はコメント行

2行目は3行目以降から入力するデータ数

3行目以降は空白文字で区切られた指定数個のデータ値

---------- 2010年度4月から ------------------------------------------

第2引数,第3引数が無い場合:

標準入力から一変量データファイル形式でsample_size個の一様乱数を読み込み，

それらから標本番号データを作成して標本を1回だけ抽出する．

読み込まれる一様乱数 u[i] (i=0,1,2,...,sample_size-1) は整数値で

以下の条件を満たすことが必要

0<=u[i]<=U : ただし，Uは母集団のサイズ以上の整数値

標本番号は mod(u[i], 母集団のサイズ) で計算される．

--------------------------------------------------------------------

***********************************************************************/

#include <stdio.h>

#include <stdlib.h>

#include <time.h>

#include <math.h>

#include <gsl/gsl_cdf.h> /* GNU Scientific Library */

/****************************** 定数 ******************************/

#define NMAX 10000 /* 一変量データの最大データ数 */

#define MMAX 100 /* 標本データの最大データ数 */

#define IMAX 10000 /* 標本抽出の繰り返し最大数 */

#define SMEAN_FILE "smean.out" /* 標本平均データの保存ファイル */

#define SVAR_FILE "svar.out" /* 標本分散データの保存ファイル */

#define ZSTAT_FILE "zstat.out" /* Z統計量データの保存ファイル */

#define CHI2_FILE "chisq.out" /* カイ二乗統計量データの保存ファイル*/

#define TSTAT_FILE "tstat.out" /* t統計量データの保存ファイル */

#define RSAMP_FILE "rsamp.out"/* 無作為抽出インデックスファイル */

#define ALPHA 0.05 /* α値 */

/****************************** 関数 ******************************/

int univar_readdata(char *, double *);

void shuffle(int *, int);

void randsamp(int, int, int *);

double univar_mean(double *, int);

double univar_var_p(double *, int);

double sample_var(double *, int);

double z_alpha(double);

double tdist_alpha(double, double);

double chisq_alpha(double, double);

void univar_savedata(double *, int, char *, char *);

void save_smean(double *, double *, int, int, int,

double, double, double, double, char *);

void save_svar(double *, int, int, int, double, double,

double, double, char *);

/**********************************************************************

main(): メイン関数

**********************************************************************/

int main(int argc, char **argv)

{

int i, j;

int sample_size, repeat_count, num_p;

int rsmp[MMAX];

int z_count, t_count, chisq_count;

double pdata[NMAX];

double sdata[MMAX];

double pmean, pvar;

double smean[IMAX], svar[IMAX], z[IMAX], chi2[IMAX], t[IMAX];

double za, chi2a_l, chi2a_r, ta;

FILE *fp;

char buff[BUFSIZ];

/* 引数のチェック */

if((argc != 4 && argc != 2)) {

fprintf(stderr, "\n Usage: mysample population_file sample_size repeat_count\n");

fprintf(stderr, " mysample population_file < urand_data\n\n");

return EXIT_FAILURE;

}

/* 標本サイズ */

if(argc == 4) {

sample_size= atoi(argv[2]);

if(sample_size <= 0) sample_size = 1;

if(sample_size > MMAX) {

fprintf(stderr, "\n ### ERROR: mysample: sample_size (%d) is too big.\n\n",

sample_size);

return EXIT_FAILURE;

}

else

sample_size = 0;

/* 繰り返し回数 */

if(argc == 4) {

repeat_count= atoi(argv[3]);

if(repeat_count <=0) repeat_count = 1;

if(repeat_count > IMAX) {

fprintf(stderr, "\n ### ERROR: mysample: repeat_count (%d) is too big.\n\n",

repeat_count);

return EXIT_FAILURE;

}

else

repeat_count = 1;

/* 母集団データファイルの読み込み */

for(i=0; i<NMAX; i++) pdata[i]=0.0;

num_p=univar_readdata(argv[1], pdata); /* 母集団サイズ */

if(num_p <= 0) return EXIT_FAILURE;

/* 母平均と母分散の計算 */

pmean = univar_mean(pdata, num_p);

pvar = univar_var_p(pdata, num_p);

/* 母平均と母分散の表示 */

printf("[ Population Parameter ]\n");

printf(" population mean : %12.5e\n", pmean);

printf(" population variance: %12.5e\n", pvar);

if(argc == 4) {

#ifdef DEBUG

/* rand()の初期化(デバッグ用) */

srand((unsigned)1);

#else

/* rand()の初期化(現在時刻で設定) */

srand((unsigned)time(NULL));

#endif

}

else { /* 一様乱数データの読み込み */

fgets(buff, BUFSIZ, stdin); /* 1行目はコメント */

fscanf(stdin,"%d", &sample_size); /* 2行目はデータ数 */

if(sample_size <= 0 || sample_size > MMAX) {

fprintf(stderr, " ### ERROR: mysample: sample_size(=%d) is invalid.\n\n", sample_size);

return EXIT_FAILURE;

}

for(i=0; i<sample_size; i++) {

if(fscanf(stdin,"%d", &j) == EOF) {

fprintf(stderr, "\n ### Error: mysample: %d-th data could not read.\n\n", i+1);

return EXIT_FAILURE;

}

rsmp[i] = (j % num_p); /* (一様乱数値)を(母集団サイズ)で割ったときの余り */

}

/* 信頼区間のパーセント点の計算 */

za = z_alpha(0.5*ALPHA);

chi2a_r = chisq_alpha(0.5*ALPHA, sample_size-1);

chi2a_l = chisq_alpha(1-0.5*ALPHA, sample_size-1);

ta = tdist_alpha(0.5*ALPHA, sample_size-1);

/* 無作為抽出インデックスファイルを開く */

if((fp=fopen(RSAMP_FILE, "w"))==NULL){

fprintf(stderr, " ### ERROR: mysample: Cannot open %s.\n\n", RSAMP_FILE);

return EXIT_FAILURE;

}

fprintf(fp, "Indexes Selected by Random Samplings (M=%6d, R=%6d)\n",

sample_size, repeat_count);

fprintf(fp, "%d\n", sample_size*repeat_count);

/* 標本抽出の繰り返し */

printf("[ Random Sampling ]\n");

printf(" Sampling Size: %d\n", sample_size);

printf(" Repeat Count : %d\n", repeat_count);

z_count = 0; t_count = 0; chisq_count = 0;

for(i=0; i<repeat_count; i++) {

/* 乱数列の生成 */

if(argc == 4) {

randsamp(num_p, sample_size, rsmp);

}

/* 標本抽出 */

for(j=0; j<sample_size; j++) {

sdata[j] = pdata[ rsmp[j] ];

/* debug */

/*printf("sdata[%4d]=pdata[%6d]=%12.5le\n", j, rsmp[j], sdata[j]);*/

fprintf(fp, "%6d ", rsmp[j]);

if( (j % 10) == 9 ) fprintf(fp, "\n");

}

if((j % 10) != 0) fprintf(fp, "\n");

/* 統計量の計算 */

smean[i] = univar_mean(sdata, sample_size); /* 標本平均 */

svar[i] = sample_var(sdata, sample_size); /* 標本分散 */

/*svar[i] *= (double)(num_p-sample_size)/(num_p-1); */ /* 標本分散の有限母集団修正 */

z[i] = (smean[i]-pmean)/sqrt(pvar/sample_size); /* 標本平均のZ統計量 */

chi2[i] = (sample_size-1)*svar[i]/pvar; /* 標本分散のCHI2統計量 */

t[i] = (smean[i]-pmean)/sqrt(svar[i]/sample_size); /* 標本平均のt統計量 */

/* 信頼区間内に母平均(=0)，母分散(=1)が入る数を数える */

if(smean[i]-za*sqrt(pvar/sample_size) <= pmean &&

smean[i]+za*sqrt(pvar/sample_size) >= pmean ) {

z_count++;

}

if(smean[i]-ta*sqrt(svar[i]/sample_size) <= pmean &&

smean[i]+ta*sqrt(svar[i]/sample_size) >= pmean ) {

t_count++;

}

if((sample_size-1)*svar[i]/chi2a_r <= pvar &&

(sample_size-1)*svar[i]/chi2a_l >= pvar ) {

chisq_count++;

}

/* debug

printf("[%4d] %12.5le, %12.5le, %12.5le, %12.5le, %12.5le\n",

i+1, smean[i], svar[i], z[i], chi2[i], t[i]);

printf(" Sampling (%6d) \r", i+1); fflush(stdout); */

}

/* printf(" Sampling (%6d) ...done.\n", i); */

/* 無作為抽出インデックスファイルを閉じる */

fclose(fp);

/* 信頼区間のパーセント点の表示 */

printf("[ Percent Point ]\n");

printf(" alpha : %f\n", ALPHA);

printf(" Z(alpha/2) : %f\n", za);

printf(" t(alpha/2, %3d) : %f\n", sample_size-1, ta);

printf(" Chisq(alpha/2, %3d) : %f\n", sample_size-1, chi2a_r);

printf(" Chisq(1-alpha/2, %3d): %f\n", sample_size-1, chi2a_l);

/* 推定される信頼区間内に母平均，母分散が含まれる標本の数の表示 */

printf("[ # of Samplings Actually Including P-mean within Its %2.0f%% CI ]\n",

(1-ALPHA)*100.);

printf(" for known p-var : %d samplings (%5.1f %% )\n",

z_count, (double)z_count/repeat_count*100. );

printf(" for unknown p-var : %d samplings (%5.1f %% )\n",

t_count, (double)t_count/repeat_count*100. );

printf("[ # of Samplings Actually Including P-var within Its %2.0f%% CI ]\n",

(1-ALPHA)*100.);

printf(" for unknown p-mean : %d samplings (%5.1f %% )\n",

chisq_count, (double)chisq_count/repeat_count*100. );

/* 標本データの出力(repeat_count == 1のときのみ) */

if(repeat_count == 1) {

printf("[ Sampled Data ]\n");

for(i=0; i<sample_size; i++) {

printf("%12.5e ", sdata[i]);

if(i % 5 == 4)

printf("\n");

}

if(sample_size % 5 != 0) printf("\n");

printf("[ Statistic ]\n");

printf(" sample mean : %12.5e\n", smean[0]);

printf(" sample variance: %12.5e\n", svar[0]);

printf(" Z-statistic : %12.5e\n", z[0]);

printf(" chisq-statistic: %12.5e\n", chi2[0]);

printf(" t-statistic : %12.5e\n", t[0]);

}

/* 信頼区間，標本分布データの保存*/

save_smean(smean, svar, repeat_count, sample_size, num_p, pmean, pvar, za, ta, SMEAN_FILE);

save_svar(svar, repeat_count, sample_size, num_p, pmean, pvar, chi2a_l, chi2a_r, SVAR_FILE);

univar_savedata(z, repeat_count, "Z-statistic", ZSTAT_FILE);

univar_savedata(chi2, repeat_count, "Chi-square-statistic", CHI2_FILE);

univar_savedata(t, repeat_count, "t-statistic", TSTAT_FILE);

/* 正常終了 */

return EXIT_SUCCESS;

}

/**********************************************************************

univar_readdata():

ファイル fn から一変量データを読み込み，一次元配列dataに格納する．

戻り値: 読み込んだデータの数

エラーが発生した場合は，終了するか，負の整数値を返す．

**********************************************************************/

int univar_readdata(char *fn, double *data)

{

FILE *fp;

int i, datasize;

char buff[BUFSIZ];

if((fp=fopen(fn, "r"))==NULL){

fprintf(stderr, "\n ### ERROR: univar_readdata(): Could not open %s.\n\n", fn);

return -1;

}

printf("[ Data File ]\n");

printf(" File Name: %s\n", fn);

/* 1行目：コメント */

if(fgets(buff, BUFSIZ, fp) == NULL) {

fprintf(stderr, "\n ### ERROR: univar_readdata(): Could not read 1st line.\n\n");

return -1;

}

printf(" Comment : %s", buff);

/* 2行目：データ数 */

if(fscanf(fp,"%d", &datasize) != 1) {

fprintf(stderr, "\n ### ERROR: univar_readdata(): Could not read data size.\n\n");

return -1;

}

if(datasize <= 0 || datasize > NMAX) {

fprintf(stderr,

"\n ### ERROR: univar_readdata(): Invalid data size (=%d) was given.\n\n", datasize);

return -1;

}

printf(" Data Size: %d\n", datasize);

/* 3行目以降 */

for(i=0; i<datasize; i++)

if(fscanf(fp, "%lf", &data[i]) == EOF) break;

fclose(fp);

if(i < datasize) {

fprintf(stderr, "\n ### WARNING: univar_readdata(): \n");

fprintf(stderr, " Read data size (%d) is less than the specified one (%d).\n\n", i, datasize);

}

/* 読み込んだデータ数を返す */

return i;

}

/**********************************************************************

shuffle():

ランダム順列の生成

rperm[0],...,rperm[n-1]をランダムに切り混ぜる

戻り値: なし

**********************************************************************/

void shuffle(int *rperm, int n)

{

int i, j, tmp;

double r;

for(i=0;i<n;i++)rperm[i]=i;

for(i=n-1;i>n:i--){

r=((double)rand()/((double)RAND_MAX+1.0));

j=(int)((i+1)*r);

tmp=rperm[i];

rperm[i]=rperm[j];

rperm[j]=tmp;

}

//////////////////////////////////////////////////

///// シャッフル(実験テキスト3.2.3節を参照) /////

//////////////////////////////////////////////////

}

/**********************************************************************

randsamp():

0,1,...,n-1からランダムにm個選ぶ

手順 :

(1) 0,1,...,n-1のランダム順列rperm[0],...,rperm[n-1] を作成

(2) rperm[0],...,rperm[n-1] から m個(m<=n)を添え字の小さい方から

順に無作為に選ぶ

(3) m個のrpermの値を rsmp[0],...,rsmp[m-1]に格納する

使用する一様乱数生成器の質が良ければ，

(1)のランダム順列は冗長な処理で不要

戻り値: なし

**********************************************************************/

void randsamp(int n, int m, int *rsmp)

{

int i, j;

int rperm[NMAX];

double r;

/* 0...n-1のランダム順列の生成 */

for(i=0; i<n; i++) rperm[i]=i;

shuffle(rperm, n);

ランダム順列rperm[0],...,rperm[n-1] からm個を無作為に選び,

その値をrsmp[0],...,rsmp[m-1]に格納する

j=m;

for(i=0;i<j;i++)

r=(double)rand()/((double)RAND_MAX+1.0);

if(r<(double)j/(double)(n-i)){

rsmp[m-j]=rperm[i];

j=j-1;

if(j<=0)return;

//////////////////////////////////////////////////

///// 未完成部分 (実験テキスト3.2.3節を参照) /////

///// テキスト3.2.3節(3) の変数m,M が /////

///// この関数のどの変数と対応するか? /////

//////////////////////////////////////////////////

}

/**********************************************************************

univar_mean():

一変量データdata[0], ..., data[n-1] の算術平均を計算

戻り値:算術平均値

**********************************************************************/

double univar_mean(double *data, int n)

{

int i;

double sum;

sum = 0.0;

for(i=0;i<n;i++){

sum=sum+data[i];

}

//////////////////////////////////////////////////

///// 平均部分 (univar.cに同じ関数がある) /////

//////////////////////////////////////////////////

return sum/n;

}

/**********************************************************************

univar_var_p():

一変量データdata[0], ..., data[n-1]の母分散を計算

戻り値:母分散値

**********************************************************************/

double univar_var_p(double *data, int n)

{

int i;

double sum, err, mean;

mean = univar_mean(data, n); /* 平均値の計算 */

sum=0;

for(i=0;i<n;i++){

err = data[i] - mean;

sum=sum+err*err;

}

//////////////////////////////////////////////////

///// 母分散 (univar.cに同じ関数がある) /////

//////////////////////////////////////////////////

return sum/n;

}

/**********************************************************************

sample_var():

一変量データdata[0], ..., data[n-1]の標本分散を計算

戻り値:標本分散値

**********************************************************************/

double sample_var(double *data, int n)

{

int i;

double sum, err, mean;

mean = univar_mean(data, n); /* 平均値の計算 */

sum=0;

for(i=0;i<n;i++){

err = data[i] - mean;

sum=sum+err*err;

}

return sum/(n-1);

//////////////////////////////////////////////////

///// 未完成部分 univar_var_p() とほぼ同じだが /////

///// 1箇所だけ異なる /////

//////////////////////////////////////////////////

}

/**********************************************************************

z_alpha():

標準正規分布の 100αパーセント点

戻り値:標準正規分布の片側・上側確率100alpha[%]に対応するz値を返す

**********************************************************************/

double z_alpha(double alpha)

{

return gsl_cdf_ugaussian_Qinv(alpha);

}

/**********************************************************************

tdist_alpha():

t分布の 100αパーセント点

戻り値:自由度dfのt分布の片側・上側確率100alpha[%]に対応するt値を返す

**********************************************************************/

double tdist_alpha(double alpha, double df)

{

return gsl_cdf_tdist_Qinv(alpha, df);

}

/**********************************************************************

chisq_alpha():

カイ二乗分布の 100αパーセント点

戻り値:自由度dfのカイ二乗分布の片側・上側確率100alpha[%]を返す

**********************************************************************/

double chisq_alpha(double alpha, double df)

{

return gsl_cdf_chisq_Qinv(alpha, df);

}

/**********************************************************************

univar_savedata():

一変量データdata[0],...,data[n-1]を

univarへの入力ファイル形式で保存する

fn は保存ファイル名

commentは1行目に記述するコメント文字列

戻り値: なし

**********************************************************************/

void univar_savedata(double *data, int n, char *comment, char *fn)

{

FILE *fp;

int i;

if((fp=fopen(fn, "w"))==NULL){

fprintf(stderr, " ### ERROR: univar_savedata(): Cannot open %s.\n\n", fn);

return;

}

/* 1行目 */

fprintf(fp, "%s\n", comment);

/* 2行目 */

fprintf(fp, "%d\n", n);

/* 3行目以降 */

for(i=0; i<n; i++)

fprintf(fp, "%12.5e\n", data[i]);

fclose(fp);

}

/**********************************************************************

save_smean():

標本平均と信頼区間の保存

smean[0],...,smean[n-1]の標本平均データと，

正規母集団を仮定する場合の母平均の信頼区間(母分散が既知と未知)を

ファイル名 fnで保存する

svar : 標本分散データ: svar[0],...,svar[n-1]

n : 標本データ数

m : 標本のサイズ

num_p: 母集団のサイズ

pmean: 母平均

pvar : 母分散

za : Z値の100αパーセント点

ta : 自由度m-1のt分布の100αパーセント点

戻り値: なし

**********************************************************************/

void save_smean(double *smean, double *svar, int n, int m, int num_p,

double pmean, double pvar, double za, double ta,

char *fn)

{

FILE *fp;

int i;

double zlim, tlim;

/* 母分散pvarが既知の場合の区間推定

smean-zlim <= pmean <= smean+zlim */

zlim =za*sqrt(pvar/m); ///// 未完成部分(実験テキスト2.4.2節) /////

if((fp=fopen(fn, "w"))==NULL){

fprintf(stderr, " ### ERROR: save_smean(): Cannot open %s.\n\n", fn);

return;

}

/* コメント行 */

fprintf(fp, "# %s: sample means (s.mean) and their confidence intervals\n", fn);

fprintf(fp, "# population size : %d\n", num_p);

fprintf(fp, "# population mean (p.mean): %12.5e\n", pmean);

fprintf(fp, "# population variance : %12.5e\n", pvar);

fprintf(fp, "# sample size : %d\n", m);

fprintf(fp, "# repeat count : %d\n", n);

fprintf(fp, "#\n");

fprintf(fp, "# confidence coefficient : %f\n", 1-ALPHA);

fprintf(fp, "# confidence interval:\n");

fprintf(fp, "# if p.var is known : s.mean - z_lim <= p.mean <= s.mean + z_lim\n");

fprintf(fp, "# if p.var is unknown: s.mean - t_lim <= p.mean <= s.mean + t_lim\n");

fprintf(fp, "#\n");

fprintf(fp, "# s.mean, z_lim, t_lim\n");

/* データ行 */

for(i=0; i<n; i++) {

/* 母分散pvarが未知の場合の区間推定

smean-tlim <= pmean <= smean+tlim */

tlim =ta*sqrt(svar[i]/m); ///// 未完成部分(実験テキスト2.4.3節) /////

fprintf(fp, "%12.5e %12.5e %12.5e\n", smean[i], zlim, tlim);

}

fclose(fp);

}

/**********************************************************************

save_svar():

標本分散と信頼区間の保存

svar[0],...,svar[n-1]の標本分散データと

正規母集団を仮定する場合の母分散の信頼区間を

ファイル名 fnで保存する

n : 標本データ数

m : 標本のサイズ

num_p : 母集団のサイズ

pmean : 母平均

pvar : 母分散

chi2a_l: 自由度m-1のカイ二乗分布の100(1-α/2)パーセント点

chi2a_r: 自由度m-1のカイ二乗分布の100α/2パーセント点

戻り値: なし

**********************************************************************/

void save_svar(double *svar, int n, int m, int num_p, double pmean,

double pvar, double chi2a_l, double chi2a_r, char *fn)

{

FILE *fp;

int i;

double lim_low, lim_up;

if((fp=fopen(fn, "w"))==NULL){

fprintf(stderr, " ### ERROR: save_svar(): Cannot open %s.\n\n", fn);

return;

}

/* コメント行 */

fprintf(fp, "# %s: sample variances (s.var) and their confidence intervals\n", fn);

fprintf(fp, "# population size : %d\n", num_p);

fprintf(fp, "# population mean : %12.5e\n", pmean);

fprintf(fp, "# population variance(p.var) : %12.5e\n", pvar);

fprintf(fp, "# sample size : %d\n", m);

fprintf(fp, "# repeat count : %d\n", n);

fprintf(fp, "#\n");

fprintf(fp, "# confidence coefficient: %f\n", 1-ALPHA);

fprintf(fp, "# confidence interval : lim_low <= p.var <= lim_up\n");

fprintf(fp, "#\n");

fprintf(fp, "# s.var, lim_low, lim_up\n");

/* データ行 */

for(i=0; i<n; i++) {

/* 母分散pvarの区間推定

lim_low <= pvar <= lim_up */

lim_low = (m-1)*svar[i]/chi2a_r;//// 未完成部分(実験テキスト2.4.4節) ///// /* 下側信頼限界 */

lim_up = (m-1)*svar[i]/chi2a_l;///// 未完成部分(実験テキスト2.4.4節) ///// /* 上側信頼限界 */

fprintf(fp, "%12.5e %12.5e %12.5e\n",

svar[i], lim_low, lim_up);

}

fclose(fp);

}