lsqcurvefit または lsqnonlin の用のコード生成

`lsqcurvefit` または `lsqnonlin` の用のコード生成

この例では、非線形最小二乗用の C コードを生成する方法を示します。

最小二乗用のデータとモデル

この例では、ベクトル xdata が 100 個のデータ点を示し、ベクトル ydata が対応する測定値を示します。xdata と ydata のモデル化された関係は以下になります。

${ydata}_{i} = a_{1} + a_{2} \exp (- a_{3} {xdata}_{i}) + ε_{i} .$

問題のデータを生成します。

rng(5489,'twister') % For reproducibility
xdata = -2*log(rand(100,1));
ydata = (ones(100,1) + .1*randn(100,1)) + (3*ones(100,1)+...
    0.5*randn(100,1)).*exp((-(2*ones(100,1)+...
    0.5*randn(100,1))).*xdata);

このコードは、平均 2 の指数分布から 100 個の独立したサンプルによって xdata を生成します。コードは、a = [1;3;2] を使用して定義方程式から ydata を生成します。これには、標準偏差 [0.1;0.5;0.5] をもつ正規偏差を追加して摂動が与えられています。

`lsqcurvefit` 用のコード生成の解決

ソルバーアプローチ

ゴールは、データに最も合うモデルのパラメーター ${\hat{a}}_{i}$ 、i = 1, 2, 3 を見つけることです。

lsqcurvefit を使用してパラメーターをデータに当てはめるには、近似関数を定義する必要があります。lsqcurvefit では、近似関数がパラメーターベクトル a とデータ xdata を取得して応答の予測を返します。この予測は、ノイズのない完全なモデルを使用した ydata と一致する必要があります。無名関数として近似関数 predicted を定義します。

predicted = @(a,xdata) a(1)*ones(100,1)+a(2)*exp(-a(3)*xdata);

モデルをデータに当てはめるためにlsqcurvefit はパラメーターの初期推定 a0 を必要とします。

a0 = [2;2;2];

lsqcurvefit を呼び出して、最適適合パラメーター ${\hat{a}}_{i}$ を求めます。

[ahat,resnorm,residual,exitflag,output,lambda,jacobian] =...
   lsqcurvefit(predicted,a0,xdata,ydata);

コード生成アプローチ

コード生成を使用して同じ問題を解くには、次の手順を実行します。

以前のステップ (データの生成、近似関数の作成、初期点の作成、および lsqcurvefit の呼び出し) をすべて組み込んだ関数を記述します。

function [x,res] = solvelsqcurve
rng(5489,'twister') % For reproducibility
xdata = -2*log(rand(100,1));
ydata = (ones(100,1) + .1*randn(100,1)) + (3*ones(100,1)+...
    0.5*randn(100,1)).*exp((-(2*ones(100,1)+...
    0.5*randn(100,1))).*xdata);
predicted = @(a,xdata) a(1)*ones(100,1)+a(2)*exp(-a(3)*xdata);
options = optimoptions('lsqcurvefit','Algorithm','levenberg-marquardt','Display','off');
a0 = [2;2;2];
lb = [];
ub = [];
[x,res] = lsqcurvefit(predicted,a0,xdata,ydata,lb,ub,options);
end

コード生成用の構成を作成します。この場合は、'mex' を使用します。
```
cfg = coder.config('mex');
```
関数 solvelsqcurve 用のコードを生成します。
```
codegen -config cfg solvelsqcurve
```
solvelsqcurve_mex.mexw64 のような名前の生成されたファイルを実行して、生成されたコードをテストします。
```
[x,res] = solvelsqcurve_mex
```
```
x =

    1.0169
    3.1444
    2.1596


res =

    7.4101
```

`lsqnonlin` 用のコード生成の解決

ソルバーアプローチ

ゴールは、データに最も合うモデルのパラメーター ${\hat{a}}_{i}$ 、i = 1, 2, 3 を見つけることです。

lsqnonlin を使用してパラメーターをデータに当てはめるには、近似関数を定義する必要があります。lsqnonlin では、近似関数がパラメーターベクトル a、データ xdata、およびデータ ydata を取得します。近似関数は、応答の予測とデータ ydata の差を返します。この差は、ノイズのない完全なモデルでは 0 になるはずです。無名関数として近似関数 predicted を定義します。

predicted = @(a)(a(1)*ones(100,1)+a(2)*exp(-a(3)*xdata) - ydata)

モデルをデータに当てはめるためにlsqnonlin はパラメーターの初期推定 a0 を必要とします。

a0 = [2;2;2];

lsqnonlin を呼び出して、最適適合パラメーター ${\hat{a}}_{i}$ を求めます。

[ahat,resnorm,residual,exitflag,output,lambda,jacobian] =...
   lsqnonlin(predicted,a0);

コード生成アプローチ

コード生成を使用して同じ問題を解くには、次の手順を実行します。

以前のステップ (データの生成、近似関数の作成、初期点の作成、および lsqnonlin の呼び出し) をすべて組み込んだ関数を記述します。

function [x,res] = solvelsqnon
rng(5489,'twister') % For reproducibility
xdata = -2*log(rand(100,1));
ydata = (ones(100,1) + .1*randn(100,1)) + (3*ones(100,1)+...
    0.5*randn(100,1)).*exp((-(2*ones(100,1)+...
    0.5*randn(100,1))).*xdata);
predicted = @(a) (a(1)*ones(100,1)+a(2)*exp(-a(3)*xdata) - ydata);
options = optimoptions('lsqnonlin','Algorithm','levenberg-marquardt','Display','off');
a0 = [2;2;2];
lb = [];
ub = [];
[x,res] = lsqnonlin(predicted,a0,lb,ub,options);
end

コード生成用の構成を作成します。この場合は、'mex' を使用します。
```
cfg = coder.config('mex');
```
関数 solvelsqnon 用のコードを生成します。
```
codegen -config cfg solvelsqnon
```
solvelsqnon_mex.mexw64 のような名前の生成されたファイルを実行して、生成されたコードをテストします。
```
[x,res] = solvelsqnon_mex
```
```
x =

    1.0169
    3.1444
    2.1596


res =

    7.4101
```
解は solvelsqcurve_mex によって生成されるものと同じになります。これは、ソルバーの基礎となるアルゴリズムが同じだからです。そのため、都合の良い方のソルバーを使用できます。