![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
このページでは、ブラウン運動の実験で得られたデータをScilabを用いて解析する方法の例を解説する。データ解析は実験課題の一部なので、解説はわざと不完全になっている。必要な情報は、ScilabのWeb pageや、マニュアルを参照して自力で探すこと。
Scilabのダウンロードページから、自分のプラットフォーム用のインストーラをダウンロードして、インストールする。
Linuxの場合、Distributionのパッケージシステムからインストール可能な事が多い。
現状では、Linux, Windows, MacOSX版が用意されているようである。
Scilabの基本文法はMatlabに良く似ている。Scilabの基本に関しては、Scilabのマニュアルでは分かりにくいので、例えばFinn's Scilab & Scicos Pageを参照すると良い。
不明な関数などについては、Scilabのプロンプトから、
help('command_name')
と打てば、ヘルプが表示される。
Scilabにおいて、//はコメントの開始を意味する。
brown.datというファイルに入っているデータを読み込むには、
d=read('brown.dat',-1,2);
とする。ここで、ファイル名の後の-1はファイルの最後の行まで読み込むことを意味し、2は2列のデータであることを意味する。詳しくはhelp('read')を参照。
t=d(:,1); v=d(:,2);
とすることで、tに時間データが、vに電圧データが代入される。
plot(t,v);
とすれば、時系列プロットが生成される。
平均と分散を求める方法は以下の通り
mean(v) //平均 variance(v) //分散
自己相関関数は、プログラム的に計算しなければならない。
c=[];tau=[]; dt=t(2)-t(1); //サンプリング間隔 for n=0:100, tau(n+1)=dt*n; c(n+1)=sum(v(1:$-n).*v(n+1:$))/(length(v)-n); // $は最後の要素を表す記号 end plot(tau,c);
Scilabにおいて、ヒストグラムをプロットするのは簡単である。以下のコマンドは、bin数20でvのヒストグラムをプロットする。
histplot(20,v)
しかし、histplot()コマンドは、Matlabのhistコマンドのようにヒストグラムの生データを返してはくれない。そこで、Matlabのhistと同様の働きをする関数を以下のように定義する。この間数は、返り値として、各binの度数nとbinの中央値xを返す。
function [n,x]=hist(nbin,v), vmax=max(v); vmin=min(v); binWidth=(vmax-vmin)/nbin; n=[]; x=[]; for i=1:nbin, binMin=vmin+(i-1)*binWidth; binMax=vmin+i*binWidth; x(i)=(binMax+binMin)/2; n(i)=sum((binMin<=v)&(v<binMax)); end n($)=n($)+sum(v==vmax); plot(x,n); endfunction
この関数を使うと、
[n,x]=hist(20,v);
のようにして、ヒストグラムの生データが得られる。
上記で得られたヒストグラムを、ガウシアンでFitしてみる。
まずは、Fitする関数型を定義する。pはFittingのパラメータである。
function y=FittingFunc(x,p), y=p(1)*exp(-(x-p(2)).^2/p(3)^2); endfunction
FittingのCriterion関数
function e=G(p,z), e=(z(2)-FittingFunc(z(1),p))^2; endfunction
これらの関数を用いて、以下のようにFittingを実行する。
p0=[2500;0;2.5]; //パラメータの初期値 (列ベクトルである必要がある) Z=[x';n']; //Fit対象のデータ [p,err]=datafit(G,Z,p0); //最適なpを探索
得られたpを用いて、結果をチェック。
nfit=FittingFunc(x,p); plot(x,n,x,nfit)
