べき乗則の近似式算出に関して

ある実験データがべき乗則に従うと考えられるので近似式を算出したいと考えております。

下記のコードで算出しようと試みておりますがどうも算出出来ないようです。

xdata=[43200,86400,345600,691200]
ydata=[1.29*10^-6,2.37*10^-6,3.36*10^-6,5.33*10^-6]

fun = @(x,xdata)x(1)*xdata.^x(2)
x0=[0,0]
x=lsqcurvefit(fun,x0,xdata,ydata)
times = linspace(10^4,10^7);
loglog(times,fun(x,times),'k-','LineWidth',2)
hold on
scatter(xdata,ydata)

エクセルで算出するとy=1*10^-8*x^0.46程度になりそうです。

一方でcurve fitting toolsを使用してべき乗則に近似すると

y=1.4*10^-9*x^0.62

となるようです。

どのようなコードを用いれば　べき乗則の近似曲線を求めることができるかお教えください。

またエクセルとMatlabで近似曲線の算出結果になぜ違いが出るのかお教えください。

Answer 1

[michio]

最適化結果が初期値（コード内の x0=[0,0]）に依存していますね。

Curve Fitting Toolbox では経験的に適当な（曖昧な表現ですが・・）初期値から係数の最適化、すなわちフィッティングを行うため、ご指摘の通り

x = [1.4e-9, 0.62]

という結果が得られていますが、また別の初期値（x0 = [0,0] や x0 = [0,1] など）で計算をすると違う結果になります。Excel 側の手法は分かりませんが、ちなみに Curve Fitting Toolbox で使用されている初期値は以下で確認できます（ x0 = [1.64e-09, 0.6165] です。）

他の方法としては y の値が非常に小さいので log をかけた値に対してフィッティングをするとロバスト性が向上するかもしれませんが、いずれにしてもある程度適切な（なんとなくオーダーが合っている）初期値を使ってやる必要があります。

xdata=[43200,86400,345600,691200];
ydata=[1.29*10^-6,2.37*10^-6,3.36*10^-6,5.33*10^-6];
ydatalog = log(ydata);
fun = @(x,xdata)x(2)*log(xdata)+log(x(1));

x0=[1e-9,1]; % この辺は結局データから適当に類推
x=lsqcurvefit(fun,x0,xdata,ydatalog)
times = linspace(10^4,10^7);
loglog(times,exp(fun(x,times)),'k-','LineWidth',2)
hold on
scatter(xdata,ydata)