2019年度数理科学科の情報処理演習Ⅱの補講課題やってみた
はじめに
今日、友達のもとに情報処理演習Ⅱの単位取得ラインに届かなかったので、課題の内容がよければ単位をあげますよというメールが届いたらしく、テスト前にもかかわらず頼ってきたのですが、明日のテストは大したものじゃないのでとりあえず片づけとくかという感じでやることにしました。
問題の内容
問題は2問あり、
- $e^x$の近似値をマクローリン展開を用いて求めるC言語の関数を作成する。
- $x$のデータ列と$y$のデータ列があり、MATLABを用いてそれの回帰直線の切片と傾き求める関数を作成し、その関数に散布図と回帰直線を同時に書かせるプログラムを入れる。
というような問題だった。ということでサクッと解いていこうと思います。
問1の考えたプログラム
まず考え方としてはe^xのマクローリン展開は\displaystyle\sum^{\infty}_{n=0}\frac{x^n}{n!}で近似できます。今回はnというマクローリン展開の項数があらかじめ指定されているということで、プログラムでの計算は\displaystyle\sum^{n}_{k=0}\frac{x^k}{k!}となります。ここで、\displaystyle\sum^{n}_{k=0}\frac{x^k}{k!}=\frac{x^n}{n!}+\sum^{n-1}_{k=0}\frac{x^k}{k!}よりn項のマクローリン展開の値を求めるとき、n-1項のマクローリン展開の値に\dfrac{x^n}{n!}を足せばよいということになり、これをうまく用いて関数を再起呼び出しを行って関数を作成していきます(授業で関数の再起呼び出しは習ったので多分こういったやり方が望まれているのではないかと予想しました)。再起呼び出しを行う際無限ループは回避しないといけないのでn=0の時は1を返して終了という形にしました。また、n項目の値である\dfrac{x^n}{n!}のx^nはxn、n!はmという変数で置きました。ということで、関数の名前をe_xということにしたので、関数からreturnするのは
return xn / m + e_x(x, n - 1);
ということになります。
次に関数の引数、戻り値の型に関して考えていきます。まず引数に関して、マクローリン展開の項数は必ず正数になるのでint型、e^xのxは特に指定は無いということで、小数も考えられるためdouble型にしました。また、戻り値に関してもe^xの値が返ってくるということでdouble型としました。以上が再起呼び出しなどを用いて作成した関数のプログラムが以下になります。このプログラムはVisualStudio2019の環境では問題なく動きました。が、学校の課題提出の環境で同様に動作するかは不明なので提出する際は参考程度にしてください。
#include
double e_x(double x, int n) {
if (n == 0) {
return 1;
}
else {
int i, m = 1;
double xn = 1;
for (i = 0; i < n; i++) {
m = m * (i + 1);
xn = xn * x;
}
return xn / m + e_x(x, n - 1);
}
}
このあとにmain関数でこの関数を実行するとしっかり出てくると思いますが、nの値が大きすぎるとdouble型の大きさでは足りなくなりinfと出てきてしまうので注意してください。これの回避方法はしりません。
問2の考えたプログラム
次はMATLABになりますが、これに関しては大変楽なことに回帰直線(ax+b)のaとbの値を計算してくれる関数が存在しているのでそれを使ってしまえば終わりなのですが、多分そんなものはこの課題で求められていないということで、今回は最小二乗法というものを使って求めていきます。参考は高校数学の美しい物語 最小二乗法の行列表現(単回帰,多変数,多項式)です。
この最小二乗法というのは、データ(x_1,y_1),\cdots ,(x_n,y_n)に対してA=\begin{pmatrix}1 & x_1\\ 1& x_2\\ \vdots & \vdots \\1 &x_n \end{pmatrix}, \overrightarrow{b}=\begin{pmatrix}y_1\\y_2\\ \vdots \\ y_n\end{pmatrix}, \overrightarrow{x}=\begin{pmatrix}\beta \\\alpha\end{pmatrix}とおくと、回帰直線(\alpha x+\beta)の\alphaと\betaに関して\overrightarrow{x}=(A^{\mathrm{T}}A)^{-1}A^{\mathrm{T}}\overrightarrow{b}で求めることができるというものです。これをそのままMATLABで書いていきます。
事前に与えられている情報に関しては
x=1:0.1:5;
y=2*x+1+0.5*randn(1,41);
というxとyなのですが、既にデータの数が41と決まっているのでそれを用いて関数を作成します。関数名はmyLine、引数はxとyのデータ、戻り値は回帰直線の傾きaと切片bということで関数側のプログラムは次のようになりました。
function [a,b] = myLine(x,y)
X=horzcat(ones(41,1),x.');
Y=(X.'*X)\(X.'*y.');
a=Y(2);
b=Y(1);
plot(x,y,'.');
hold on;
fplot(@(k)Y(2)*k+Y(1),[0,6]);
end
あまりうまく出来たようには見えませんがとりあえず求められているものに関しては正しく返せているので補講課題としてはこのようなものでいいのではないのかなと思いました。このプログラムは上の数学的操作をそのままコードに直したものなので特に補足事項はありません。流れを簡単に解説すると、最初に上の解説のAを作成し、xを求め、a,bに返し、最後に散布図をプロットした状態で回帰直線も表示するようにしています。fplotの下の行にhold off;と入れたほうがよかったかもしれませんが、まあ問題ないかなと勝手に思っているだけなので、提出する際は入れておいたほうが無難なんじゃないかなって思います。
終わりに
以上で数理科学科2019年度情報処理演習Ⅱの補講課題のプログラムに関して作り方に関して説明してきたんですが、とりあえずこの程度のプログラムなんて大したことないので頑張って作成しましょう。来年度に再履修したらまためんどくさい課題が待ってるのでやる気のある人だけ再履修の選択肢をとるのが賢明だと思います。そして、このページの存在が教授達にしられてしまうと問題になる可能性があるので、教授達にはばれないようにしてください。
コメント入力