{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 2022年度プログラミング演習A・B\n",
"\n",
"# 第15回レポート課題の解説"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 演習1\n",
"\n",
"次の $11$ 個の点に対して、正規方程式を解いて最小二乗法による線形回帰を行い、それらの点および求めた回帰直線（範囲は $0\\leq x\\leq 1$）を描画してください。\n",
"\n",
"$$\n",
"(0,1.09),\\quad (0.1,1.17),\\quad (0.2,1.42),\\quad (0.3,1.54),\n",
"$$\n",
"\n",
"$$\n",
"(0.4,1.71),\\quad (0.5,1.99),\\quad (0.6,2.26),\\quad (0.7,2.44),\n",
"$$\n",
"\n",
"$$\n",
"(0.8,2.55),\\quad (0.9,2.84),\\quad (1,3.04)\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%plot --format svg\n",
"\n",
"data = [0,1.09; 0.1,1.17; 0.2,1.42; 0.3,1.54; 0.4,1.71; 0.5,1.99; 0.6,2.26; 0.7,2.44; 0.8,2.55; 0.9,2.84; 1,3.04]; %点を並べた行列\n",
"n = size(data,1); %点の個数\n",
"\n",
"plot(data(:,1),data(:,2),\".\") %点の描画\n",
"grid on\n",
"hold on\n",
"\n",
"A = [ones(n,1),data(:,1)];\n",
"b = data(:,2);\n",
"c = (A'*A)\\(A'*b); %正規方程式を解く\n",
"%実は「c = A\\b」でもよい（https://jp.mathworks.com/help/matlab/math/systems-of-linear-equations.htmlを参照）\n",
"\n",
"x = 0:0.01:1;\n",
"y = c(1)+c(2)*x;\n",
"plot(x,y) %回帰直線の描画"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 演習2\n",
"\n",
"演習1と同じ $11$ 個の点に対して、統計量を用いた計算で最小二乗法による線形回帰を行い、それらの点および求めた回帰直線（範囲は $0\\leq x\\leq 1$）を描画してください。\n",
"\n",
"ただし、統計量の計算には平均は`mean(ベクトル)`、分散は`var(ベクトル)`、共分散は`cov(ベクトル1,ベクトル2)`を使用してよいです。"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%plot --format svg\n",
"\n",
"data = [0,1.09; 0.1,1.17; 0.2,1.42; 0.3,1.54; 0.4,1.71; 0.5,1.99; 0.6,2.26; 0.7,2.44; 0.8,2.55; 0.9,2.84; 1,3.04]; %点を並べた行列\n",
"\n",
"plot(data(:,1),data(:,2),\".\") %点の描画\n",
"grid on\n",
"hold on\n",
"\n",
"c2 = cov(data(:,1),data(:,2))/var(data(:,1)); %統計量を用いた計算\n",
"c1 = mean(data(:,2))-c2*mean(data(:,1));\n",
"\n",
"x = 0:0.01:1;\n",
"y = c1+c2*x;\n",
"plot(x,y) %回帰直線の描画"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 演習3（オプション）\n",
"\n",
"節点\n",
"\n",
"$$\n",
"x_i=-1+\\frac{2(i-1)}{n-1}\\quad (i=1,\\ldots,n)\n",
"$$\n",
"\n",
"を用いて関数\n",
"\n",
"$$\n",
"f(x)=\\frac{1}{1+25x^2}\n",
"$$\n",
"\n",
"を $m-1$ 次以下の多項式\n",
"\n",
"$$\n",
"g(x)=\\sum_{j=1}^m c_j x^{j-1}\n",
"$$\n",
"\n",
"で近似することを考えます。\n",
"\n",
"$n=30,60$、$m=10,20$ を組み合わせた計 $4$ 個の場合について近似多項式を求め、元の関数 $f(x)$ とともに $-1\\leq x\\leq 1$ の範囲でグラフを描画してください。"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%plot --format svg\n",
"\n",
"function value = f(x) %関数f(x)\n",
" value = 1./(1+25*x.^2);\n",
"end\n",
"\n",
"x = -1:0.01:1; %x座標の配列\n",
"y = f(x);\n",
"plot(x,y,\"DisplayName\",\"f(x)\") %f(x)のグラフを描画\n",
"grid on\n",
"hold on\n",
"\n",
"%2つのfor文でnとmの値を変化させる\n",
"for n = [30,60]\n",
" node = (-1:2/(n-1):1)'; %節点の配列（縦ベクトル）\n",
" for m = [10,20]\n",
" A = zeros(n,m);\n",
" for j = 1:m\n",
" A(:,j) = node.^(j-1); %行列Aの各列を計算\n",
" end\n",
" b = f(node); %縦ベクトルbを定める\n",
" c = (A'*A)\\(A'*b); %正規方程式を解いて係数ベクトルcを求める\n",
" y = 0;\n",
" for j = 1:m\n",
" y = y+c(j)*x.^(j-1); %近似多項式を計算\n",
" end\n",
" plot(x,y,\"DisplayName\",[\"n=\",num2str(n),\", m=\",num2str(m)]) %近似多項式のグラフを描画\n",
" end\n",
"end\n",
"\n",
"legend"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Octave",
"language": "octave",
"name": "octave"
},
"language_info": {
"file_extension": ".m",
"help_links": [
{
"text": "GNU Octave",
"url": "https://www.gnu.org/software/octave/support.html"
},
{
"text": "Octave Kernel",
"url": "https://github.com/Calysto/octave_kernel"
},
{
"text": "MetaKernel Magics",
"url": "https://metakernel.readthedocs.io/en/latest/source/README.html"
}
],
"mimetype": "text/x-octave",
"name": "octave",
"version": "5.2.0"
}
},
"nbformat": 4,
"nbformat_minor": 4
}