自分のブログ名

物理の空き地 by M.E_K

日々の学び、感じたことを書いております。ブログ移行中->https://physics-mek.com

移転しました。

リダイレクトします。

最小二乗法解くなら行列が楽って話

実験式がy=a+bxであるとき

最小二乗法の原理より,誤差y_i -(a+bx_i)の二乗の和を作ると,同じ精度を持つ測定の場合

\displaystyle{S=\sum_{i=1}^n {y_i - (a+bx_i)}=0}
\begin{align}
\therefore & \frac{\partial S}{\partial a} = -2\sum_{i=1}^{n}{\{y_i - (a+bx_i)}\}=0\\
& \frac{\partial S}{\partial b} = -2\sum_{i=1}^{n}{\{y_i - (a+bx_i)}\}x_i=0\end{align}
nはデータセットの個数である.
\begin{align} N&a + [x]b = [y]\\
 [x]&a+[x^2]b=[yx]\end{align}

行列式に置き換えて,
\begin{bmatrix}N & [x] \\
[x] & [x^2]
\end{bmatrix}\begin{bmatrix}a \\ b \end{bmatrix}
= \begin{bmatrix}[y] \\ [yx]\end{bmatrix}

\begin{bmatrix}a \\ b \end{bmatrix}
= \begin{bmatrix}N & [x] \\
[x] & [x^2]
\end{bmatrix}^{-1}\begin{bmatrix}[y] \\ [yx]\end{bmatrix}\cdots (1)
(1)を見ると,逆行列を計算することと,行列積を求めるだけでパラメータa,bが決定できる.

データ

f:id:Monster_K:20210917160528p:plain:w175:h280f:id:Monster_K:20210917160533p:plain:w320:h71

エクセルで逆行列

逆行列を計算する関数はMINVERSE関数.求まる逆行列の大きさに等しいセル範囲を選択した状態で下記のように数式を入力して、[Ctrl]+[Shift]+[Enter]を押すと結果が得られる.
f:id:Monster_K:20210917160530p:plain:w462:h116

エクセルで行列積

行列の積を計算する関数はMMULT関数.求まる逆行列の大きさに等しいセル範囲を選択した状態で下記のように数式を入力して、[Ctrl]+[Shift]+[Enter]を押すと結果が得られる.
f:id:Monster_K:20210921152149p:plain

結果

f:id:Monster_K:20210921153828p:plain
y = 0.66 + 0.50xとなった.

最確値a,bの確率誤差の計算

y_iに平均誤差\mux_iの誤差がないものとすれば,
\mu = \sqrt{\frac{\sum{\{y_i -(a+bx_i)\}^2}}{n-2}}=\sqrt{\frac{[y^2]-a[y]-b[yx]}{n-2}}
分母である自由度の個数はポイント数nからパラメータの個数を引いたものである.
ゆえに,y_iの確率誤差\sigma_y
\sigma_y =\sqrt{\frac{[y^2]-a[y]-b[yx]}{n-2}}

a,bの確率誤差\sigma_a,\sigma_bは,誤差の伝搬法より,
\sigma_a = \sqrt{\sigma_y^2 \sum_{i=1}^{n}(\frac{\partial a}{\partial y_i})^2}
\sigma_b = \sqrt{\sigma_y^2 \sum_{i=1}^{n}(\frac{\partial b}{\partial y_i})^2}

エクセルで計算するには

\sigma_y =\sqrt{\frac{[y^2]-a[y]-b[yx]}{n-2}}の計算は非常に単純なので,割愛.
\sigma_a = \sqrt{\sigma_y^2 \sum_{i=1}^{n}(\frac{\partial a}{\partial y_i})^2}
\sigma_b = \sqrt{\sigma_y^2 \sum_{i=1}^{n}(\frac{\partial b}{\partial y_i})^2}
の計算を行列で表すと,
\begin{bmatrix}\sigma_a&0\\
0&\sigma_b
\end{bmatrix}
=\sqrt{\frac{\sigma_y^{2}}{△}\begin{bmatrix}N & [x] \\
[x] & [x^2]\end{bmatrix}^{-1}}

MINVERSE関数は分数の△も含むので,\sigma_yの計算だけすればよく,エクセルで計算すると,
f:id:Monster_K:20211025113406p:plain

結果の対角成分が確率誤差になる.
y = 0.66\pm 0.09 + (0.50\pm0.01)x

まとめ

f:id:Monster_K:20211025124116p:plain
エクセルのシートを載せた.最小二乗法を連立方程式で解けばそれまでではあるが,エクセルだとどうしても計算ミスが怖い.
高次の計算こそ行列で済ませるととても早くなるので,ぜひ活用していただきたい.
「行列でやってみたい」「やってみてもうまくいかない」となった人はぜひ私に聞きに来てほしい.

C/C++ double型からint型へ変換

(int)で型を変換すると,切り下げになる.

#include<iostream>
int main() {
	double d = 10.811;
	int i;
	i = (int)d;#doubleからintへ変換している
	std::cout << i << std::endl;
}
#output
10

そのため,四捨五入がしたい場合は,double型に対して0.5を足してdouble型の状態で繰り上げておく必要がある.

#include<iostream>
int main() {
	double d = 10.811;
	int i;
	i = (int)(d +0.5);#double型に対して0.5を足してint型に変換する
	std::cout << i << std::endl;
}
#output
11

C/C++で二次関数のFitting

【やったこと】
C/C++で二次関数の最小二乗法を行った.計算では,逆行列・行列の積をやった.


【わかったこと】
C/C++にはPythonのような最小二乗法のライブラリはない
(おそらく,あるなら教えてほしい)
行列は簡単(手計算は地獄).ちなみにfor文でぶんまわすだけ.


【計算式】
二次関数の最小二乗法
y=a{x}^{2}+bx+c


 \begin{pmatrix}
a\\b\\c
\end{pmatrix} = 
{\begin{pmatrix}
[{x}^{4}]&[{x}^{3}]&[{x}^{2}]\\
[{x}^{3}]&[{x}^{2}]&[{x}]\\
[{x}^{2}]&[{x}]&N \end{pmatrix}}^{-1}
\begin{pmatrix}
[{x}^{2}y] \\ [xy]\\ [y]
\end{pmatrix}

【コード】

double A[3][3] = { {x4,x3,x2},{x3,x2,x},{x2,x,N} };
double inv_A[3][3], c[3], buf, term;
int i, j, k;//for文のための準備
int n = 3;//3×3行列用
//単位行列を作る
for (i = 0; i < n; i++) {
	for (j = 0; j < n; j++) {
		inv_A[i][j] = (i == j) ? 1.0 : 0.0;
	}
}
//逆行列を計算する
for (i = 0; i < n; i++) {
	buf = 1 / A[i][i];
	for (j = 0; j < n; j++) {
		A[i][j] *= buf;
		inv_A[i][j] *= buf;
	}
	for (j = 0; j < n; j++) {
		if (i != j) {
			buf = A[j][i];
			for (k = 0; k < n; k++) {
				A[j][k] -= A[i][k] * buf;
				inv_A[j][k] -= inv_A[i][k] * buf;
			}
		}
	}
}
//行列の積
for (i = 0; i < n; i++) {
	term = 0;
	for (k = 0; k < n; k++) {
		term = term + inv_A[i][k] * B[k];		
	}
	C[i] = term;
}
//結果の出力
for (i = 0; i < n; i++) {
		std::cout << C[i] << " ";
}
std::cout << std::endl;

【今回使ったデータセット
https://gitlab.com/-/snippets/2163095/raw/main/sample.txt?inline=false

【結果】

13.0000 22.0148 -10.0188

f:id:Monster_K:20210817220611p:plain

VSCodeでPythonコードが実行できない対処法

原因VSCodepython.exeのパスを読めていない.


対処法:設定ファイルを開きたい.以下の順に開いていき,jsonファイルを開く.
File->Preferences->Setting->OpenSetting

{
      "files.encoding": "shiftjis",
      "workbench.editorAssociations":{
      "*.ipynb"."jupyter-notebook"
    }
}

こんな感じのjsonファイルを開いて,以下を追加

"python.pythonPath":"C:\\Users\\name\\AppData\\Local\\Programs\\Python37\\python.exe"
//nameは自分のusernameを使う

IrfanView おすすめの設定

「IrfanView」定番の画像ビューワー - 窓の杜
IrfanViewpngやjpgを開くのにおすすめのソフトフェアである.
しかし,初期設定では画像を最後まで開くと,
こんな感じで毎回聞かれるためだるい.


そこで,画像送りをループ設定にするとよい.
Options->Properties/Setting->Browising/Editingを開き,Loop current folderにチェック.

C/C++ filesystem ディレクトリから指定拡張子をvectorに詰める

pythonのglob関数をC/C++で使いたい.
それができるといわれる"glob.h"なるものがあるらしい.しかし,それはLinuxのライブラリらしくWindowsでは使えそうになかった.(私の知識不足感は否めない)そこで,filesystemと呼ばれるライブラリでも同じようにglob的なことができるようだ.
本記事では,ディレクトリからpngファイルを取得し,vectorに詰めるコードを記す.

#include <iostream>
#include <string>
#include <vector>
#include <filesystem>

int main()
{
	std::string dir_name = "C:/Users/MEK/test";
	std::vector<std::string> files;
	std::filesystem::path dir_path(dir_name);
//例外処理
	if (!std::filesystem::is_directory(dir_path)) {
		std::cout << "No Dir" << std::endl;
		return -1;
	}
//pngファイルをvectorに詰める
	auto dir = std::filesystem::directory_iterator(dir_path);
	for (auto& p : dir) {
		std::string pathname = p.path().string();
		if (pathname.find(".png") != std::string::npos) {
			std::cout << p.path().string() << std::endl;
			files.push_back(pathname);
			std::cout << std::endl;
		}
	}
}

注意事項
VisualStudio2019では,プロジェクトプロパティの全般からC++標準をC++17にする必要がある.