Rのテキスト 知能情報処理演習１サポートサイト R Text

2013 年度 知能情報処理演習 1

R による科学技術計算

D1 クラス担当：北野

D2 クラス担当：谷口

平成 25 年 5 月 15 日

1kitano@ci.ritsumei.ac.jp

2taniguchi@ci.ritsumei.ac.jp

第 ⁷ 週  ^R の基礎 その ¹

Rの基礎的な使い方を練習する。マニュアル的な内容はここでは記載しないので、関連 資料¹²を参考にすること。また本演習のサポートサイト³も参照のこと。

コンソール環境

コンソールでは、プロンプトに関数（コマンド）を呼び出して実行したり、演算を行う ことができる。

課題₁ 起動時に表示される関数_demo()や_help()を実行せよ。

課題₂ 四則演算、累乗、剰余、数学関数（_sin(x)など）の演算を実行せよ。

オブジェクトと代入

Rでは、数値、配列、関数など処理の対象を格納するものをオブジェクトと呼ぶ。オブ ジェクトに代入を実行することで、オブジェクトを生成することができる（_Cで変数の宣 言と初期化を同時に行うのと同様）。また、オブジェクトに対し演算処理を施すことも出 来る。

課題₃ オブジェクト_xを数値、文字、_NULL（空＝何もなし）で初期化せよ。 課題₄ 初期化した_xに対し、課題₂で行った演算を実行せよ。

ベクトル

Rでは、複数のデータをひとまとまりのベクトルとして扱うことができる。_Cにおける 1^{次元配列まるごとが}1つのベクトルに対応する。上で扱った要素が₁つのデータ（_Cに おける変数）は、₁次元のベクトルと考える。

課題₅ ベクトルデータを生成する関数_c()を用いて、オブジェクト_xをベクトルとして 初期化せよ（要素数はいくつでもよい）。

課題₆ 初期化したベクトル_xに対し、課題₂で行った演算を実行し、その結果を確認せよ。

1RjpWiki: http://www.okada.jp.org/RWiki

2

上記サイトのリンク集からレベル別に応じた様々な資料を辿ることができる。

3https://sites.google.com/site/ciexercise/

ファイルへの保存とファイルからの実行

プログラムをコンソールのプロンプトに入力して実行するだけでなく、テキストファイ ルに編集して保存し、ファイルから実行することも出来る。

課題₇ 下の手順に従い、ファイルへの保存とファイルからの実行を確かめよ。

1. 「ファイル」メニューから「新しいスクリプト」を選択し、エディタを起動。 2. エディタでプログラムを編集。

3. 「ファイル」メニューでファイルに保存（例えば、ファイル名を”program.R”^）。 4. コンソールのプロンプトでsource("C:/program.R")^{4として、ファイルから}

実行。

4

作業ディレクトリに注意すること。

第 ⁸ 週  ^R の基礎 その ²

制御構造、関数、グラフ作成の方法について練習を行う。

処理の制御

条件分岐（_if-else文、_switch文）や反復（_for文）などの制御文を用いて、また、頻 繁に用いる処理は新しい関数を定義して、プログラムの制御を行える。

課題_{1 if-else}文を用い、適当な整数値を与えられたオブジェクト_xが、奇数の場合に

は_xの₂乗を、偶数の場合には_xの₂倍を計算する処理を実行せよ。 課題_{2 for}文を用いて、

∑₁₀₀₀

n=1 n¹ を計算するプログラムを作成せよ。

課題₃ 引数_xを₂乗する関数_square(x)を定義し、実行するプログラムを書け。

グラフ作成

Rではグラフ作成ツールも充実しており、様々なグラフを簡単に作成できる。

課題₄ ベクトルデータ_xを上で作成した関数_square()に適用し、オブジェクト_yに代入 せよ。得られた_x、_y座標を関数_plotを用いてグラフにせよ。

課題₅ 数学関数_cos(x)を_{[−2π, 2π]}の範囲で図示せよ。

第 ⁹ 週 行列と連立方程式

行列の操作、演算および連立方程式の解について練習する。

行列の作成と操作

課題₁ ベクトルデータx = (1, 2, 3, 4, 5, 6)^{を作成し、関数}matrix()^{を用いて、} 2 × 3^行列、3 × 2^{行列を作成せよ。}

課題₂ 正方行列を₂つ作成し、それらの和、差、積、スカラー積（これらは演算子に注 意）、を求めよ。

課題₃ 上で作成した行列に対し、行列式、転置行列、逆行列を求めよ。

課題₄ 関数_eigen()を用いて、固有値・固有ベクトルを求めよ。関数_eigen()は、固有

値と固有ベクトルの₂つの要素（ベクトル）をリスト⁵として返す。リスト内の複数 のベクトルのうち、₁つを取り出すには、その名前＝ラベルを指定すればよい。例え ば、固有値の場合、固有値に付けられたラベル_$valuesを用い、

> y <- eigen(x)

> y$values

とする。

課題₅ 正方行列_Aを作成し、その固有ベクトルを並べてできる行列_Pと_P⁻¹を求め、行 列_Aの対角化_P⁻¹_APを行え。

連立方程式の解法

A^を行列、x^、b^{をベクトルとする時、}Ax = bの形に表される連立方程式の解_xは、_R に用意された関数で簡単に求めることが出来る。

課題₆ 行列_A、ベクトル_bを定め、関数_solve()を用いて解を求めよ。

52つ以上のベクトル一纏めにしたもの

第 ¹⁰ 週 回帰分析

測定データなどの数値の組の間の関係性をモデル式に当てはめて理解することを回帰分 析という。モデル式＝回帰式は、実データと回帰式の間の誤差が出来る限り小さくなるよ うに決定する。誤差を二乗誤差で定義する場合、この回帰式決定の方法として最小二乗法 が用いられる。

線形回帰分析と最小二乗法

データの組_{(xi, yi})(i = 1, 2, · · · , n)^{に対し、モデル式}y = f (x) = ax + b^でx_i^とy_i^の関 係が説明できるか⁶を検討するために、_{f (xi)}と_yiの差が出来る限り小さくなるように係 数_aと_bを求めたい。今、

G =













 1 x1

1 x₂ ... ^... 1 xn















, a= (b, a)^T, y= (y₁, y₂, · · · , yn^{)T (1)}

と定義し、_{f (xi)}と_yiの差を二乗誤差J = (y − Ga)^T(y − Ga)^{で定め、Jを最小にするよ} うな係数を求める。_Jを_aの関数_J(a)と見れば、_J(a)が最小になるのは、^∂J

∂a ^{= 0を満た}

す時である（_Jを_aの各成分で偏微分したものが₀）。これをみたす_aは、_G^T_{Ga = G}^T_y と表される（_G^T_Gに逆行列が存在すれば、_{a= (G}^T_G)⁻¹_G^T_yと書ける）。_a、_bは、成分 表示すると、

a = ⁿ

∑n

i=1^xiyi−^∑ni=1^xi∑ni=1^yi

n^∑n_i=1x²_{i − (}^∑n_i=1x_i)^{2 (2)} b =

∑_n

i=1^x2i

∑_n

i=1^yi−^∑ni=1^xiyi

∑_n

i=1^xi

n^∑n_i=1x²_{i − (}^∑n_i=1xi)^{2 (3)} (4) と表せる（簡単な計算なので導出してみるとよい）。

組み込みデータセット

• 組み込みデータセット一覧の表示_{: data()}

• 個々のデータセットの表示_{: "}データセット名_"（例えば_cars）

演習課題

課題₁ データセット_carsは、自動車の速度_speedとその速度でブレーキをかけた場合の 停止するまでの距離_distの組を表すデータである（サンプル数₅₀）。このデータ セットについて、_speedを_x、_distを_yとして、まず散布図を作成せよ。次に回帰 式を求め、その回帰直線をデータの散布図に重ねて表示せよ。

6

より正確には、_xiで_yiを説明できるか

課題₂ データセット_anscombeは複数の組を格納したデータセットである。例えば、_x1と y1^、x2^とy2が対応している。まず、対応する系列の組の散布図を作成せよ。次に、 それぞれの組に対し回帰式を求め、回帰直線をデータともに図示してみよ。さらに ４つ組で比較してみよ。

課題₃ 適当なデータセットを選び、回帰分析を行え（回帰分析に適さないデータもある ことに注意）。

第 ¹¹ 週 非線形方程式の解法

非線形関数（_xの₁次式で表せない関数）_{f (x)}に対し、非線形方程式_{f (x) = 0}を満た す解_α（_{f (α) = 0}）を数値的に求める問題扱う。

ニュートン法

非線形方程式の代表的な解法であるニュートン法は以下の手続きで反復的に解を求める 方法である。_{f (x) = 0}の真の解_αの正確な値は分かっていないが、_αに近い値_x0は分かっ ているとする。この_x0を初期値とし、_n回目の反復により得られる_αに対するよりよい 近似解を_xnとするとき、_{n + 1}回目の近似解は次のように与えられる。

x_n+1 = x_n^{―f (xn)}

f^′(xn)⁽⁵⁾

ここで、_f^′_(xn)は_{f (x)}の_xについての導関数に_xnを代入したものである。これは、近似 解_xnが与えられているとき、

1. x = x_n^におけるy = f (x)^{の接線の方程式}y = f^′(x_n)(x^―x_n) + f (x_n)^を求め、 2. ^{この接線と}y = 0^との交点x_n−f^{f (x}_′(xⁿ⁾

n) ^を得、

3. この求めた交点を次の近似解とする、

として理解できる。具体例を図に示す。ここでは_{f (x) = x}³₋³

2^x2−32^{x + 1、初期値x0 = 1}

としている。反復により、_{x1 =} ¹

3^{, x2 =} 117^{59,と真の解 1}2に近づいて行くことがわかる。 

-1.5 -1 -0.5 0 0.5 1 1.5

-0.5 0 0.5 1 1.5 2 2.5

y = x³ –(3/2)x² –(3/2)x +1 y = –(3/2)x +1/2 y = –(13/6)x +59/54

x₀ = 1 x1 = 1/3

x2 = 59/117

図 _1: ニュートン法の具体例

反復の停止条件は_ϵを十分に小さな値（₁₀⁻⁶など）とし

• |xn+1− xn| < ϵ

• |^xn+1_xn^−xn| < ϵ

• |f(xn)| < ϵ

など適当なものを設定する。注意すべきは、方程式_{f (x) = 0}が₂つ以上の解を持つ場合 でも、１つの初期値から始める反復では、１つの数値解しか得ることができない。

ニュートン法のアルゴリズム

1. ^{適当な初期値}x₀^を選ぶ。ϵを十分小さな正の値とし、停止条件を設定する。 2. ^上式(5)^によりx_n^からx_n+1^{を求める。}

3. 停止条件を満たせば終了、そうでなければ_nをインクリメントして₂へ。

演習課題

課題₁ ニュートン法を実行するプログラムを作成せよ（解の系列を表示するように作成す ること）。また、そのプログラムを用いて、_{f (x) = x}³₋³

2^x2−32^{x + 1とし、f (x) = 0}

の解を求めよ。

課題_{2 f (x) = e}^−x_{+ x}²_{− 6x}とし、_{f (x) = 0}の数値解を求めよ。 課題_{3 f (x) = xe}^−x/2_{− e}^−xとし、_{f (x) = 0}の数値解を求めよ。

第 ¹² 週 微分方程式の解法

未知の変数_xが満たす関係式は方程式と呼ばれ、その方程式を解くことで変数_xの値を 求めることができた。未知の関数_y(x)およびその導関数（_y^′_{, y}^′′_{, · · ·}）が満たす関係式は 微分方程式と呼ばれ、これを解くことで関数_y(x)を決定することができる。常微分方程 式を数値的に解く方法を扱う。

オイラー法

変数_xの関数_yについての微分方程式と初期条件 dy

dx ^{= y}

′= f (x, y), y(a) = A (6)

を代表的な数値的解法であるオイラー法で解く手順は以下の通りである。

y(x)の解を仮定し（まだ得られていないので）、任意の_xのまわりでテイラー展開す ると、

y(x + ∆x) = y(x) + y^′(x)∆x +¹ 2^y

′′_(x)(∆x)2₊¹

6^y

′′′_(x)(∆x)3_{+ · · · (7)}

と書ける。元の問題₍₆₎は連続値で表されているるが、計算機で扱う場合には、離散化し て扱わなければならない。_∆xは変数_xを離散化した時の間隔（＝刻み幅、分点間隔）であ る。つまり、変数_xを分点間隔_∆xで離散化した場合の_i番目の分点_xiは、_{xi= a + i∆x} と書ける。また、_xiにおける_y(xi)の“ 数値解 ” は通常誤差を含むので、真の値_y(xi)と 区別して_Yiと書くことにする。

∆xが十分に小さいとすると、式₍₇₎において、

• ^{右辺第3項以降は、第1、2}項に比べ無視出来るほどに小さくなるので無視、

• ^{式(6)より、y′ = f (x, y)を右辺第2項に代入、}

• x = xi^{とすると、x + ∆x = x}i+1^、

となることから、式₍₆₎の微分方程式は、

Y_i+1= Y_i+ ∆xf (x_i, Y_i) (8)

と表すことができる。式₍₈₎を用い、_{x0 = a}、_{Y0 = y(x0}) = y(a) = A^{を初期値として、逐} 次_Yiを求めていけば良い。

オイラー法を微分方程式とその初期条件_{dy/dx = y, y(0) = 1}に適用した例を示す。 この微分方程式は解析的に解くことができ、_{y = e}^xがその解である（普通は解析解はわか らないことが多い）。図の実線が解析解、点列が数値解_Yiを表している。誤差を含んだ_Yi を用いて_Yi+1を計算するので、_iが増加するにつれ_e^xiと_Yiの誤差が大きくなっていくこ とがわかる。なるべく誤差を小さくするには_∆xを小さくとればよい。

1 1.5 2 2.5 3

0 0.2 0.4 0.6 0.8 1

y

x

dy/dx = y (y = e^-x)

(x₁, Y₁) = (0.2, 1.2)

(x₂, Y₂) = (0.4, 1.44)

(x3, Y3) = (0.6, 1.728) (x₀, Y₀) = (0, 1), Δx = 0.2

図_2: オイラー法の具体例

オイラー法のアルゴリズム

1. i = 0^、x₀ = a^、Y₀ = A^{と初期化、反復回数}N^と∆x^{を定める。} 2. x_i+1= xi + ∆x^、Y_i+1= Y_i+ ∆xf (x_i, Y_i)^を計算

3. i^{をインクリメントし}i = Nならば終了、そうでなければ₂へ。

演習課題

課題₁ オイラー法を実行するプログラムを作成せよ。さらに、微分方程式

dy

dx ^{= yを初期}

条件_{y(0) = 1}とし、_{0 ≤ x ≤ 1}の範囲で解を求めよ。また、刻み幅_∆xを変え、解 析解_{y = e}^xとの誤差を比較せよ₍グラフを作成するとよい₎。

課題₂ 微分方程式dy/dx = y(1 − ky), ^{y(0) = 0.1について、定数kを適当に設定して} 解を求めよ。また、数値解をグラフにし、_y(x)の特徴について考察せよ。

課題₃ 微分方程式dy/dx = −y + cos x, ^{y(0) = 0について、定数kを適当に設定して解} を求めよ。また、数値解をグラフにし、_y(x)の特徴について考察せよ。

第 ¹³ 週 乱数を用いる計算 ーモンテカルロ法ー

計算機で生成する疑似乱数（以下、単に乱数と呼ぶ）を用いることで様々な計算、シミュ レーションを行うことができる。乱数を利用する計算アルゴリズムはモンテカルロ法と呼 ばれ、その代表的な計算をいくつか扱う。

数値積分

2次元平面において面積を求めたい領域_S、領域_Sを含む正方形の領域を領域_Aを考え る。領域_A内に生成した一様乱数の個数を_NA、そのうち領域_S内に含まれたものの個数 を_NSとする。領域_Aの面積_Aと領域_Sの面積_Sの比_{A : S}は、_{NA: NS}に近似する。こ れを利用して、半径₁の円の面積あるいは円周率_πの近似値を求めることができる。

1. xy^{平面において、}(0, 0)^、(1, 0)^、(1, 1)^、(0, 1)^{を頂点とする正方形}A^{と原点を中} 心とする半径₁の円の第₁象限にある_1/4円を領域_Sを定める。

2. r_x^とr_y^{は互いに独立な}[0, 1]の範囲の一様乱数とし、これらを用いてランダムな点 (r_x, r_y)^をN_A^{個生成する。}

3. ^{このうち、領域}S^{内に入る点（}r_x²+ r²_{y ≤ 1}^{を満たす点）の個数}N_S^{をカウントする。} 4. ^面積S^はN_S/N_A^{、円周率は}4N_S/N_Aで求めることができる。

任意の確率分布に従う乱数の生成：棄却法

一様分布や正規分布など標準的な確率分布ではない任意の確率分布に従う乱数の生成は、 次の関数を用意する。

• p(x)：生成したい乱数が従う確率密度関数（数式で記述出来るものはもちろん、数

値的にしか表現できないものでも可）

• q(x)：乱数の生成が簡単な確率密度関数（一様分布、正規分布など）

• k^{：定義域において} p(x) ≤ kq(x)^{が成り立つ定数} これらを用い、₁個の乱数を生成する実行手順は以下の通り。

1. q(x)^{に従い、乱数}x^を生成。

2. [0, kq(x)]^{の範囲で一様乱数}u^を生成。

3. u ≤ p(x)^ならxを採択（乱数として使用）、そうでなければ棄却（使用しない）。

p(x)^としてp(x) = (x² _{− x +} ¹₄)/Z^{（ただし、定義域は}[0, 1]^；Z ^{は規格化定数、}Z =

∫1

0^(x2− x +¹₄^{)dx）、q(x)として一様分布、k = 1/4とした例を示す。}

Histogram of sample

sample

Frequency

0.0 0.2 0.4 0.6 0.8 1.0

01000200030004000

図3: p(x) = (x²_{− x +}¹₄)/Zに従う乱数のヒストグラム

演習課題

課題₁ モンテカルロ法による数値積分の方法で、円周率_πの近似値を求めよ。乱数の個 数を変えた場合に、近似値がどのようになるか確かめよ。

課題₂ モンテカルロ法による数値積分の方法で、^x

2

a^{2 +} y² b^{2 +z}

2

c^{2 ≤ 1（a、b、cは正の数）を}

満たす領域の体積を求めよ。_a、_b、_cは各自で設定すること。

課題₃ 棄却法を用いて、p(x) = x(1 − x)/Z, (0 ≤ x ≤ 1)を確率密度関数とする乱数を 生成し、そのヒストグラムを作成せよ。

第 ¹⁴ 週 検定

データの統計的性質に関する仮説の確からしさを定量的に確かめる仮説検定について、 実際にデータを扱うことにより、その考え方を理解する。

仮説検定

観測対象（母集団）の統計的性質について立てた予想（仮説_H）に対し、その観測対象 から実際に得られたデータの集合（標本）が予想から著しく離れた場合、その予想の正し さを疑うことになる。この「著しく」を定量的に評価し、予想の正しさを見積もることが 仮説検定である。

 コイントスにより得られる、表、裏それぞれの確率を_p、_{1 − p}とするとき、そのコイ ンが公正（＝表裏の確率が同じ）かどうかについての検定を行う。そのため、コイントス を₁₀₀回行い、表₆₀回、裏₄₀回の結果を得た。

• ^仮説H0は、コインの表が出る確率は_{p = 1/2}である。

• コイントスで表が出る回数_yは確率変数であり、仮説_H0に従う場合、その確率分布 は二項分布B(100, 1/2)^{で表すことができる。（ただし、}B(n, p) =_nC_kp^k_{(1 − p)}^n−k）

• ^{二項分布B(n, p)はnp}が十分大きいとき、ガウス分布N (np, np(1 − p))^{で近似でき} る。つまり、仮説_H0の下でy ∼ N(50, 25)^である。

• ˆ^{p = y/nと定義すると、}p ∼ N(p, p(1 − p)/n)^{ˆ 。あるいは、Z = (ˆ}p − p)/^√p(1 − p)/n とすると、Z ∼ N(0, 1)^。

• コイントスの結果から、_{Z = 2}。この値が起こる確率は_0.05以下、つまり、確率_0.05 以下の稀な現象が起こった（仮説_H0の下では通常起こりえない）。生じた事実（コイ ントスの結果）は覆すようがないので、仮説_H0が誤り（有意水準_0.05の時）。よっ て仮説_H0は棄却。

課題₁ 適当な値の_pを設定し、コイントスをシミュレートして、仮説検定を行え。

標本平均の分布

X ∼ N(µ, σ²⁾とする。この確率分布から_n個サンプルした標本データ_{X1, X2}, · · · , Xⁿ の平均（標本平均）_{X =}¯ ¹

n

∑n

i=1^Xi^{、分散7v2 =} _n−1^{1 ∑n}i=1^(Xi− ¯^{X)2も確率変数である。} n^{が十分大きければ、}_{v ∼ σ}^であり、_{X ∼ N(µ, σ}^{¯ 2}/n)^{となる。あるいは}t = √^X−µ¯

v²/n^に対

し、t ∼ N(0, 1)^となる。

課題_{2 X ∼ N(µ, σ}²₎とする。この確率分布から_n個サンプルした標本データ_{X1, X2}, · · · , Xn

の平均（標本平均）_{X =}^{¯ 1}

n

∑n

i=1^{Xi、分散8v2 =} _n−1^{1 ∑n}i=1^(Xi− ¯^{X)2も確率変数で}

7

不偏分散

8

不偏分散

ある。_nが十 分大きけ れば、_{v ∼ σ}であ り、_{X ∼ N(µ, σ}¯ ²_/n)とな る。あるい は t = √^X−µ¯

v²/n^に対し、t ∼ N(0, 1)。サンプリングを繰り返して_tの分布をヒストグラム に表せ。_nの値を変え、特に_nが小さい場合（₁桁の値）に、ヒストグラムと_{N (0, 1)} を比較せよ。

t _検定

標本に対し、その母平均_µについての検定を行う場合、通常その分散（母分散）_σ²につ いては未知である場合がほとんどである。標本の数_nが十分多ければ、標本から求めた分 散（不偏分散）_vを_σと見なしても問題ない。しかし、課題₂で見たように、_nが小さい 場合には、_tは標準正規分布_{N (0, 1)}からずれることがわかる。つまり、

X−µ¯

√σ²/n ^̸= X−µ¯

√v²/n^。

この場合、標準正規分布ではなく、_{t =}

X−µ¯

√v²/n が従う分布（スチューデントの_t分布と呼 ぶ）を用いて検定しなければならない。_t分布を用いた検定を_t検定（検定すべきは母平 均_µ）と呼ぶ⁹。

課題₃ サポートサイト¹⁰からテキストファイル”t-test.dat”をダウンロードし、格納さ れたデータの母平均について検定せよ。ファイルからの読み込みには

x <- read.table("^{ファイルへのパス}/^{ファイル名}")

として、_xにリストとして格納することができる。

9t^{分布の形、}t検定のより詳細な内容については、₃回生開講の「統計データ解析」にて学ぶ。

10https://sites.google.com/site/ciexercise/

発展課題 ¹

サポートサイトに掲載された課題のうち₁つを選び、実施せよ。

発展課題 2

本演習で扱った手法（あるいはその発展的手法でも可）を用いて、各自で用意したデー タ、あるいはテーマについて解析し、考察せよ。