* 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)

トップ   編集 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 検索 最終更新   ヘルプ   最終更新のRSS