* Scilabを用いたブラウン運動データ解析 [#b5d36144]
このページでは、ブラウン運動の実験で得られたデータをScilabを用いて解析する方法の例を解説する。データ解析は実験課題の一部なので、解説はわざと不完全になっている。必要な情報は、Scilabの[[Web page>http://www.scilab.org/]]や、[[マニュアル>http://www.scilab.org/download/index_download.php?page=documentation]]を参照して自力で探すこと。
#contents
** インストール [#w9646c91]
Scilabの[[ダウンロードページ>http://www.scilab.org/download/index_download.php]]から、自分のプラットフォーム用のインストーラをダウンロードして、インストールする。
Linuxの場合、Distributionのパッケージシステムからインストール可能な事が多い。~
現状では、Linux, Windows, MacOSX版が用意されているようである。
** Scilabの基本 [#efe731f4]
Scilabの基本文法はMatlabに良く似ている。Scilabの基本に関しては、Scilabのマニュアルでは分かりにくいので、例えば[[Finn's Scilab & Scicos Page>http://home.hit.no/~finnh/scilab_scicos/]]を参照すると良い。
不明な関数などについては、Scilabのプロンプトから、
help('command_name')
と打てば、ヘルプが表示される。
Scilabにおいて、//はコメントの開始を意味する。
** データの読み込み [#h4162120]
brown.datというファイルに入っているデータを読み込むには、
d=read('brown.dat',-1,2);
とする。ここで、ファイル名の後の-1はファイルの最後の行まで読み込むことを意味し、2は2列のデータであることを意味する。詳しくはhelp('read')を参照。
t=d(:,1);
v=d(:,2);
とすることで、tに時間データが、vに電圧データが代入される。
plot(t,v);
とすれば、時系列プロットが生成される。
** 平均、分散、自己相関関数 [#ef9062ea]
平均と分散を求める方法は以下の通り
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);
c(n+1)=sum(v(1:$-n).*v(n+1:$))/(length(v)-n); // $は最後の要素を表す記号
end
plot(tau,c);
** ヒストグラム [#k8aecd1c]
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);
のようにして、ヒストグラムの生データが得られる。
** フィッティング [#u0dac302]
上記で得られたヒストグラムを、ガウシアンで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)
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
