近似曲線をMATLABで書く!!

本記事はMATLAB R2021aを用いています.

# Excelの最小二乗法による近似曲線の種類
Excelの近似曲線は主に以下がある.
![excel.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1635383/ea54d56d-0994-e354-9cbf-dd73ce3e7af4.png)

* 線形近似
* 多項式近似
* 対数近似
* 指数近似
* 累乗近似

今回はMATLABでこれらを描いていく.

# MATLABでデータを作る
今回は近似を行うデータもMATLABで作っていく. もし, Excelのdataを用いたい場合はMATLABの[readtable関数](https://jp.mathworks.com/help/matlab/ref/readtable.html)を使えば, 簡単にデータをインポートできる.

今回は. 近似のターゲットとなる関数に[randn関数](https://jp.mathworks.com/help/matlab/ref/randn.html?s_tid=doc_ta)による誤差を足すことによってデータを作った.

# MATLABでの近似曲線の実装
ここから, MATLABで近似曲線を作っていく. 主に参考となるリンクは以下となる.

**近似曲線を計算する**

* [polyfit](https://jp.mathworks.com/help/matlab/ref/polyfit.html?s_tid=doc_ta)
* [polyval](https://jp.mathworks.com/help/matlab/ref/polyval.html)
* [プログラムによる近似](https://jp.mathworks.com/help/matlab/data_analysis/programmatic-fitting.html)

**グラフを描く**

* [plot](https://jp.mathworks.com/help/matlab/ref/plot.html?s_tid=doc_ta)
* [semilogx](https://jp.mathworks.com/help/matlab/ref/semilogx.html?s_tid=doc_ta)
* [semilogy](https://jp.mathworks.com/help/matlab/ref/semilogy.html)

**グラフのいろいろ**

* [xlim](https://jp.mathworks.com/help/matlab/ref/xlim.html?s_tid=doc_ta)
* [ylim](https://jp.mathworks.com/help/matlab/ref/ylim.html)
* [grid](https://jp.mathworks.com/help/matlab/ref/grid.html)
* [xlabel](https://jp.mathworks.com/help/matlab/ref/xlabel.html)
* [ylabel](https://jp.mathworks.com/help/matlab/ref/ylabel.html)
* [legend](https://jp.mathworks.com/help/matlab/ref/legend.html?s_tid=doc_ta)
* [title](https://jp.mathworks.com/help/matlab/ref/title.html?s_tid=doc_ta)
* [MATLABで作れるグラフの種類](https://jp.mathworks.com/help/matlab/creating_plots/types-of-matlab-plots.html)

## 最小二乗法による線形近似
今回のターゲット関数は, 
$$ y = 1-x $$
として, 以下のデータを作った.

``` code:matlab
data_num = 100;
data_x = 10*rand(1,data_num);
data_y = 1 - data_x + 0.5*randn(1,data_num);
```

![linear_regression_data.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1635383/ce710af1-dbfd-3ca5-2a59-4c26b999f303.png)

MATLABで最小二乗法による近似は[polyfit関数](https://jp.mathworks.com/help/matlab/ref/polyfit.html?s_tid=doc_ta)と[polyval関数](https://jp.mathworks.com/help/matlab/ref/polyval.html?s_tid=doc_ta)で行うことができる.

``` code:matlab
x = linspace(0,10,100);
p = polyfit(data_x,data_y,1);
y = polyval(p,x);
```

![linear_regression_plot.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1635383/2c99bd9b-531a-afbb-1936-e9921b0d7a17.png)

ちなみに, この近似曲線の式は,
$$ y = p(1)x + p(2) $$
となっている.

## 最小二乗法による多項式近似
ここでのターゲット関数は,
$$ y = 2 + x - \frac{1}{2}x^2 - x^3 $$
とした. そして, 以下のデータを作った.

```code:matlab
order = 3;
data_num = 100;
a = [2, 1, -1/2, -1];
data_x = 10*rand(1,data_num);
data_y = zeros(1,data_num);
for i=1:order+1
    data_y = data_y + a(i)*data_x.^(i-1);
end
data_y = data_y + 5*randn(1,data_num);
```

![polynomial_regression_data.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1635383/ebdbeeba-753b-4324-0f4f-649069ca3145.png)

MATLABでは, 多項式近似は線形近似とほぼ何も変わらない. [polyfit関数](https://jp.mathworks.com/help/matlab/ref/polyfit.html?s_tid=doc_ta)は,

```matlab
p = polyfit(xデータ, yデータ,多項式の次数);
```

となっているため, 今回は,

```code:matlab
x = 0:10/99:10;
p = polyfit(data_x,data_y,3);
y = polyval(p,x);
```

により, 以下のような近似曲線が得られる.

![polynomial_regression_plot.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1635383/629bf31d-2742-606d-0099-9b3c7db31804.png)

ちなみに, 近似曲線の式は, 
$$ y = p(1)x^3+p(2)x^2+p(3)x+p(4) $$
となる.

## 最小二乗法による対数近似
ここで得たい, 近似曲線は,

$$ y = a + b\log x $$

であり, この$a$と$b$を求める.

まず, ターゲット関数は,

$$ y = 1 -\log x $$

として以下のデータを作成した.

```code:matlab
data_num = 100;
a = [1,-1];
data_x = 10*rand(1,data_num);
data_y = a(1) + a(2)*log(data_x) + 0.1*randn(1,data_num);
```

![log_approx_data.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1635383/13cbe643-7f75-8d02-b15e-d640115d6a5e.png)


そして, このような非線形の近似曲線を求める場合は, [MATLABプログラムでの近似](https://jp.mathworks.com/help/matlab/data_analysis/programmatic-fitting.html#:~:text=%E5%B8%B8%E3%81%AB%E9%87%8D%E8%A6%81%E3%81%A7%E3%81%99%E3%80%82-,%E9%9D%9E%E5%A4%9A%E9%A0%85%E5%BC%8F%E9%A0%85%E3%82%92%E3%82%82%E3%81%A4%E7%B7%9A%E5%BD%A2%E3%83%A2%E3%83%87%E3%83%AB,-MATLAB%20Online%20%E3%81%A7)を用いて, 以下のように書ける.

```code:matlab
x = 0.01:10/99:10;
X = [ones(data_num,1), log(data_x')];
p = X\data_y';
y = [ones(size(x')), log(x')]*p;
```

このプログラム中にあるバックスラッシュは, [線形方程式](https://jp.mathworks.com/help/matlab/math/systems-of-linear-equations.html)の最小二乗解を求める. この手順を踏むことで, 以下の近似曲線を得られる.

![log_approx_plot.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1635383/4db75561-a7b5-b46d-7cd7-54b7e2d54a24.png)

また, 対数近似により得られた曲線は片対数グラフにおいて直線になる.

![log_approx_semix.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1635383/bce28ab3-4453-e99c-a9d2-f90e3bbd2472.png)

## 最小二乗法による指数近似
指数近似では,
$$ y = ae^{bx} $$
の$a$および$b$を求める. この式の両辺について対数をとると,
$$ \log y = \log a + bx $$
となり, 線形近似に変換することができる. 今回は, この性質を用いて近似曲線を求める.

今回は, ターゲット関数を,
$$ y = 10e^{-\frac{1}{4}x} $$
として, 以下のようなデータを作った.

```code:matlab
data_num = 100;
a = [10,-1/4];
data_x = 10*rand(1,data_num);
data_y = a(1).*exp(a(2).*data_x) + 0.1*randn(1,data_num);
```

![exp_approx_data.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1635383/5af4197c-87d9-5647-cce6-73d9c706651a.png)

ここでは線形近似を行うため, [polyfit関数](https://jp.mathworks.com/help/matlab/ref/polyfit.html?s_tid=doc_ta)を用いるが, 式変形に従い,

```code:matlab
x = 0:10/99:10;
p = polyfit(data_x,log(data_y),1); % data_yではなくlog(data_y)にする!!
y = exp(p(2)+p(1).*x);
```

として, 以下の近似曲線を得る.

![exp_approx_plot.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1635383/37e74b11-75f6-db8b-0847-b0f0b56d257d.png)

また, 指数近似についても片対数グラフにおいて直線となる.

![exp_approx_semiy.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1635383/0c8a85a6-392a-d326-5151-479382439dc4.png)

## 最小二乗法による累乗近似
累乗近似では,

$$ y=ax^b $$

の$a$および$b$を求める. こちらについても両辺の対数をとり,

$$ \log y = \log a + b\log x $$

となる. これを用いて対数近似と同様の作業を行ってもよいが, 今回は, この式を

$$ y' = s + b x' $$

($y' = \log y$, $s = \log a$, $x'=\log x$)とみて, [polyfit関数](https://jp.mathworks.com/help/matlab/ref/polyfit.html?s_tid=doc_ta)での線形近似を行う.

ここでのターゲット関数は,

$$ y = 10x^{-\frac{1}{4}} $$

として以下のようにデータを作った.

```code:matlab
data_num = 100;
a = [10,-1/4];
data_x = 10*rand(1,data_num);
data_y = a(1).*data_x.^a(2) + 0.1*randn(1,data_num);
```

![exponent_approx_data.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1635383/d3a6f618-91ea-3b83-9d14-147d5ccdec02.png)

そして, [polyfit関数](https://jp.mathworks.com/help/matlab/ref/polyfit.html?s_tid=doc_ta)を用いて, 近似曲線を求める.

```code:matlab
x = 0.01:10/99:10;
p = polyfit(log(data_x),log(data_y),1); % data_xとdata_yをlog(data_x)とlog(data_y)とする
y = exp(p(2)).*x.^p(1);
```

![exponent_approx_plot.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1635383/b817ea18-cc9c-ffa0-aeed-37fa06999f52.png)

# おわりに
コロナ禍により, 試験が減りレポート課題が増えてきているように感じます(公平性の担保が簡単だからなのかな..)
レポート課題が増えると, 大学の先生はたくさんの似たようなレポートを見ることになると思います(たぶん..)
こんなとき, 目で見て理解できるグラフなどがあると, それだけであなたのレポートは目立つでしょう!! 
[MATLABはたくさんのグラフを作れる](https://jp.mathworks.com/help/matlab/creating_plots/types-of-matlab-plots.html)とても優秀なツールなので, みんなでグラフを作っていい点もぎ取りましょう!!(作りすぎ&載せすぎには注意です.)