2016年5月23日月曜日

Bashスクリプトで中心極限定理

中心極限定理実験 Bashスクリプト。学生さんの(スクリプトの書き方も含めた) 勉強用。
# Japanese is written in UTF-8.
# シャープはコメントです。

##勉強用制御
No1=0
No2=0
No3=0


##色々設定
bin=50 # bin面体サイコロ
dice_test=10000 # サイコロのチェックのために振る回数
one_test=10 # one_test回振って平均を求める
test=1000 # 平均を求める試行回数

##出力ファイル
dice_check="dice_check_bin${bin}_test${diec_test}"
dice_check_txt="${dice_check}.txt"
dice_check_ps="${dice_check}.ps"

CLT_check="CLT_check_bin${bin}_test${test}"
CLT_check_txt="${CLT_check}.txt"
CLT_check_ps="${CLT_check}.ps"

rm -f $dice_check_txt
rm -f $dice_check_ps
rm -f $CLT_check_txt
rm -f $CLT_check_ps


#さいころ関数
dice_func(){
 echo "$RANDOM" $bin  | awk '{print int($1/32767*$2)}'
}

#サイコロ関数のテスト
echo "Dice is" `dice_func 0`"."



####
if (($No1))  ;
then
####
 #配列の初期化1
 array=()
 for((i=0;i do
  array[${i}]=0
 done
 
 
 #サイコロ関数の正しさをtest回ふって確認。
 for((i=0;i do
  rand=`dice_func 0`
  array[${rand}]=$((${array[${rand}]}+1))
 done
 
 #結果をファイルに書き込み
 for((i=0;i do
  echo "$i ${array[${i}]}" >> ${dice_check_txt}
 done
 echo ${dice_check_txt} "is made."
###
fi #No1
###


####
if (($No2)) ;
then
####

 #配列の初期化2
 for((i=0;i do
  array[${i}]=0
 done
 
 
 #「サイコロ関数をone_test回ふってその平均を調べる。」試行をtest回実施。
 for((i=0;i do
  s=0
  for((j=0;j  do
   rand=`dice_func 0`
   s=$((${s}+${rand}))
  done
  ave=`echo $s $one_test | awk '{print int($1/$2)}'`
  array[${ave}]=$((${array[${ave}]}+1))
  echo $i "kai end."
 done
 
 #結果をファイルに書き込み
 for((i=0;i do
  echo "$i ${array[${i}]}" >> ${CLT_check_txt}
 done
 echo ${CLT_check_txt} "is made." 
 
###
fi #No2
###



####
if (($No3)) ;
then
####
 #初期値を探す
 n0=`awk 'BEGIN{max_n=0}{if($2>max_n) max_n=$2}END{print max_n}' ${CLT_check_txt}`
 m0=`awk 'BEGIN{max_m=0;max_n=0}{if($2>max_n){max_n=$2;max_m=$1}}END{print max_m}' ${CLT_check_txt}`
 
 #Gnuplotで作図
gnuplot<
set ylabel "Number"
set xlabel "Value"
set te po co solid

set title "Dice Check : Bin ${bin} Test ${dice_test}"
set ou "${dice_check_ps}"
pl [][0:] "${dice_check_txt}" w histeps lw 2
set ou


Gauss(x)=n*exp(-0.5*((x-m)/s)**2)

n=${n0}
m=${m0}
s=m*0.1

fit Gauss(x) "${CLT_check_txt}" via n,m,s

set title "CLT Check : Ave. of ${one_test} Bin ${bin} (Test ${dice_test})"
set ou "${CLT_check_ps}"
pl "${CLT_check_txt}" w histeps lw 2, Gauss(x)
set ou

exit
EOF

###
fi #No3
###