BrownScilab - Ando Lab Wiki

Scilabを用いたブラウン運動データ解析

このページでは、ブラウン運動の実験で得られたデータをScilabを用いて解析する方法の例を解説する。データ解析は実験課題の一部なので、解説はわざと不完全になっている。必要な情報は、ScilabのWeb pageや、マニュアルを参照して自力で探すこと。

インストール

Scilabのダウンロードページから、自分のプラットフォーム用のインストーラをダウンロードして、インストールする。 Linuxの場合、Distributionのパッケージシステムからインストール可能な事が多い。
現状では、Linux, Windows, MacOSX版が用意されているようである。

Scilabの基本

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)

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2013-07-09 (火) 14:50:24 (2316d)