Tips for MATLAB

お役立ちリンク

「メモ代わり」さんの「個人的MATLAB tips」：かゆいところに手がとどくことが書いてある！神経科学者っぽいです。

こまごましたこと

判別用

isnan NaN があるかどうか

isinf Inf or -inf があるかどうか

isinteger 整数値があるかどうか

デバッグ用

error('test') エラーメッセージが出て処理が止まる

warning('test') warning メッセージが出る。処理は止まらない

pause(n) n 秒処理が止まる

fprintf('model %d 終了 \n', model) コマンドにメッセージを出す

構造体

parameter.kf

parameter.kb

という感じで書くと、parameterと1つの変数に二種類の値を入れられる。便利。

並列化のメモ

といっても parfor を用いた並列化。Matlab Distributed Computing Serverは使ったことがありません。それでも数倍シミュレーションが早くなると感動します。

• Hyper-threading対応のCPUなら、4 core/cpuが8 core/cpu になりたくさんコアが使える。BIOSをいじろう。
• 全コアをparforでまわしちゃうと、ほかの仕事ができなくなる・・・
• 並列化計算の意味をよく考えて、並列計算にすべきところをParforにする。Parforで回るときは、完全にバラバラに回る。したがって、そのparforループ内は互いに独立に完結しないといけない。
  • 例えば、初期条件を数個ふってみて、それをparfor内で並列に計算させることは可能。なぜならそれらの初期条件の間で情報のやり取りは無いから。
  • あるループの情報が次のループに必要な場合とかは並列化無理。例えばオイラー法のように、逐次結果を更新していくような計算の場合は並列計算できない。
• Parforは、以下のように宣言する
  • parfor n = 1:1:10 って感じ。もちろん、n が整数であること。
• Parforのなかの変数、データの取り扱いに注意する。
  • parfor内の計算結果を出力するためには、sliced variableにする。
  • 例えば、
  • data = cell(10,1);
  • parfor n = 1:10
  • y = ode23(...)
  • data{n} = y;
  • end
  • な感じ。 n のところが sliced variable っぽい？
• セル配列にするとsliced variableにしやすい気がする。

if matlabpool('size') ~= 0

matlabpool close;

end

非線形Fitting（式が分かっている場合）

例として、解離定数の式。

1. [Kd, r, J, COVB, mse] = nlinfit(xdata, ydata, @(Kd, x) x./(Kd(1)+x), Kd0)
2. % Kd パラメーター; r 残差; J ヤコビアン; COVB 共分散行列, mse 誤差項の推定分散
3. ci = nlparci(Kd, r, 'Jacobian', J)
4. % ci 95%信頼区間（デフォルト値）;

非線形Fitting（式が分かっていない場合。例えば微分方程式の解析解がないときなど。）

局所的な最適化でよければ fminsearch を使う。

このとき、CellDesignerで出力した.mファイルをode functionファイルとして用いる場合は、functionのなかにfunctionを使う入れ子関数にして、パラメーターをふり、自分で最適化関数を作って最小化する。

なかなかめんどくさい・・・

一例

function fitpara = optimization_test()

% 最適化のプログラム例

% 入れ子関数を使って、解析解の無い微分方程式の最適化をする。

% 構成としては、入れ子関数を使っている。

% optimization_test 実行ファイル。functionファイル

% -model.m ODEを含むモデル部分。functionファイル。

% -costfcn.m コスト関数。function。

% 実験例

% In vitro の ERK リン酸化のパラメーターを最適化により求める。

% 逐次反応モデル

% [npERK]' = -k1*[npERK]

% [pYERK]' = k1*[npERK]-k2*[pYERK]

% [pTpYERK]' = k2*[pYERK]

% ちなみに、このモデルは解析解があるのでfminsearchしなくても実は最適化できる・・・

% 実験結果（column 1, npERK; column 2, pYERK; column 3, pTpYERK）

init = ...

[2.20096 0 0

2.117027895 0.067465619 0.016466485

2.058372739 0.142102058 0.000485203

1.542959495 0.627330246 0.030670259

1.050852802 1.045774858 0.10433234

0.490494103 1.344476356 0.365989541

0.17618356 1.30477363 0.72000281

0.055793023 0.839737568 1.305429409];

% 実験結果の反応時間の値 (単位は秒）

tspan = [0;50;100;300;600;1200;1800;3600];

% 初期値パラメーター

x0 = [init(1,1) init(1,2) init(1,3)];

% 最適化処理のための初期値設定

para = [0.01 0.01]; % この初期値近傍から最適化パラメーターを探索する。

% 最適化処理のためのオプション。自分で調べて。

options = optimset('Display','OFF','PlotFcns',@optimplotfval,'TolX',1e-3,'TolFun',1e-3,'MaxIter',250);

% 最適化処理。

[fitpara, norm] = fminsearch(@costfcn, para, options);

% コスト関数。errorの値を最小にするようにサーチする。

function error = costfcn(para)

k1 = para(1);

k2 = para(2);

% ここで ode 関数を使う。

[t, fittedData] = ode23(@model,tspan,x0);

% 入れ子関数。関数の中に関数を入れる。これだと、ode関数の中の変数を最適化できる。

function xdot = model(time,fittedData)

xdot = zeros(3,1);

% 初期値の代入

npERK = fittedData(1);

pYERK = fittedData(2);

pTpYERK = fittedData(3);

%npERK

xdot(1) = -k1*npERK;

%pYERK

xdot(2) = k1*npERK-k2*pYERK;

%pTpYERK

xdot(3) = k2*pYERK;

end

% 目的関数。この関数を最小化するようにパラメーターをフィッティングする。

% つまり最小二乗法。

% sum一回だけでは3行それぞれの和が出るので、

% 二回かけることで、その三行分のまとめた和が出る。

error = sum(sum((fittedData(:,1:3)-init(:,1:3)).^2));

end

% グラフ化

figure(1)

plot(tspan, init(:,1),'bo')

hold on

plot(tspan, init(:,2),'go')

plot(tspan, init(:,3),'ro')

plot(tspan, fittedData(:,1),'b-')

plot(tspan, fittedData(:,2),'g-')

plot(tspan, fittedData(:,3),'r-')

hold off

xlabel('Time (sec)');

ylabel('ERK concentration [/\muM]');

end

グラフ用のoptions （自分用）

% Figureウィンドウの大きさ指定。結構便利。

% gcf は current figure のプロパティを指定するときに使う

set(gcf,'Position',[50,300,800,300]);

%================================================

% plot

%================================================

% グラフ化

h = plot(X_axis, Y_axis);

% 一般

h = plot(x, y, 'linestyle_marker_color','linewidhth',2)

% ラインスタイル ('-', '--', ':'); など

% マーカタイプ ('x', '*', 'o','+'); など

% カラー識別子('c', 'm', 'y', 'k', 'r', 'g', 'b', 'w');

h = plot(x, y, 'go--'); % 緑の破線、丸マーカつき

% 線幅、色

set(h,'linewidth',2,'color','blue');

% X, Y 軸

xlabel('Time(min)');

ylabel('kcat-cRaf [/\muM/sec]');

xlim([0 1]);

ylim([0 2]);

axis([0 12 -0.5 1]);

% gca は current axis のプロパティを指定するときに使う

set(gca,'XScale','Log'); % X 軸片ログ

set(gca,'YScale','Log'); % Y 軸片ログ

% タイトル

title(['ERK-' char(C(c)) ']);

% レジェンド

legend(h,'First','Second','Third');

%================================================

% hist

%================================================

histではなくBar関数を使うと色や重なりなどを変更しやすくなる。

%================================================

% pcolor (Heat mapを作るときに使う。)

%================================================

% pcolor 表示

% X_axis (１行X列ベクトル), Y_axis (Y行１列ベクトル), C(X,Y) 行列。

% C の次元をちゃんとX,Yに合うようにすること。

h = pcolor(X_axis, Y_axis, C);

% ラインスタイルをなしにする。ほぼ必須。

set(h,'linestyle','none');

% グラフの横にカラーバーつける。必須

colorbar;

% カラーバーの上限、下限を決める

caxis([0, 0.74]);

% 新しいMatlabでは色が変なのでjetにする

colormap('jet')

% 通常は左下が原点になるが、axis ijをいれると左上が原点になる。

axis ij;

% pcolorで行列データを格子状に表現する。こちらを参考にさせていただきました。

A = [1,0,0;0,1,0;0,0.5,0.7]; % 視覚化対象の行列

[x,y] = size(A); % 行列の大きさを測る

pcolor([A,zeros(x,1);zeros(1,y+1)]); % 行列サイズを１行１列を大きくしてプロット

Saveするときのデータに日付を入れる

データをsaveするときに日付を入れてoverwriteを防ぐ

save(['allData' datestr(now, 'yyyymmdd_HHMMSS')], '-v7.3');

で、allData20121223_210517.matというファイル名で保存される。

2012年12月23日21時5分17秒って感じ。

データ型

Matlabとか以前の話ですが・・・

• Zeros: デフォルトで double 。
• int8: -128 ~ 127。 めちゃ軽くなる。

CellDesignerから出力された.mファイルを使ってシミュレーションする

• 下の方にあるいらない行を削る
  • function z = pow ....以降削る。
  • ファイル名と関数名（mファイルの一行目）を一致させる。
• ode関数でまわす
  • 関数名のみタイプすると初期値が得られる（引数nargin）
  • [y, t] = ode23(@function_name, [0 100], x0)
  • % x0は初期値
• ode23やode45で遅い場合は、微分方程式がstiffである可能性が高い（陽解法）。ode15s （陰解法）とかにすると爆速になったりする。計算誤差が大きくなるみたいだけど、今のところ誤差はほとんど感じない。