Past, Present and Future
Past, Present and Future
コンテスト中に自分のブログから該当するエントリを探す作業をよくしていることに気付いたので、インデックスを作成しておくことにした。
Google Code Jam 2012 Round 1A
- Kingdom Rush
- std::sort()のコンパレータオブジェクトについて考察
Understanding How std::string::substr works
- Understanding How std::string::substr works
- std::string::substr()について考察
Recalling Bayes' Theorem
定理を素早く思い出したりするために自分にとっての「定型」の問題を用意しておくことがままあります。
例えばベイズの定理の場合。「囚人問題」や「モンティホール問題」が特に有名なので、それを「定型」としている方もいらっしゃるのではないでしょうか。
私が「定型」としている問題の一つは道具としてのベイズ統計に掲載されていた以下の問題です(第2章 3. 壷の問題を考える。p. 51)。
二つの壷a、bがある。壷aには赤玉が3個、白玉が2個入っている。壷bには赤玉が8個、白玉が4個入っている。壷aと壷bが選ばれる割合は1:2とする。どちらかの壷から玉1個を取り出したとき、それが赤玉であった。その赤玉が壷aから選ばれている確率を求めよ。
道具としてのベイズ統計
「壷aと壷bが選ばれる割合は1:2とする」の一節はどのように解釈するものなのだろう...? などと疑問に思い続けてはいますが、まあ壷は目の届かない場所にでも置いてあるんだろうななどとといいように解釈しつつ、長年付き合っています。
さて、壷aを選ぶ事象をA、どちらかの壷から玉1個を取り出したときにそれが赤玉である事象をRとします。問題が求めているのは、「どちらかの壷から玉1個を取り出したとき、それが赤玉であった。その赤玉が壷aから選ばれている確率」です。これは条件付き確率Pr(A | R)となります。
さて、ここでベイズの定理を書き出してみましょう。
同時確率Pr(A, R)は条件付き確率Pr(A | R)と周辺確率Pr(R)の積として計算できるのでした。同様に、同時確率Pr(R, A)は条件付き確率Pr(R | A)と周辺確率Pr(A)の積として計算できるのでした。
同時確率Pr(A, R)とPr(R, A)は等しいので、以下の等式が成り立ちます。
両辺をPr(R)で割って、
ベイズの定理を導くことができました。ここまでは簡単ですね。では、それぞれの確率に値を入れて計算してみましょう。
...うん、ここで少し立ち止まってみましょう。単純に確率の値を入れる前に、それぞれの確率がベイズの定理ではどのように呼ばれるかを復習しておきます。そのほうが効率的です。
左辺Pr(A | R)は求める確率です。
式を見れば一目瞭然、条件付き確率です。「どちらかの壷から玉1個を取り出したら赤玉であった」ときに「選んだ壷が壷aであった」確率です。ベイズの定理ではこれを事後確率(Posterior Probability)と呼びます。
右辺Pr(A)は「壷aが選ばれる」確率です。
ベイズの定理ではこれを事前確率(Prior Probability)と呼びます。右辺の中では一番重要です。壷の中がどのような状態になっているか知らなくても、その選ばれ方に関して何かしら事前に知っている情報だからです(どうです、重要そうでしょう?)。
右辺Pr(R | A)も条件付き確率です。
ベイズの定理ではこれを尤度(Likelihood。「ゆうど」と読みます)と呼びます。
この場合の尤度は、「壷aが選ばれた」ときに「取り出した玉が赤玉であることはどのくらいであるか」を示しますが、あまりこのように記載されることはありません。どちらかといえば、「壷aが選ばれた」ときに「取り出した玉が赤玉であることはどのくらいもっとも(尤も)らしいか」といったように記載されます。
最後に、右辺のPr(R)は規格化定数(Normalize Constant。証拠(Evidence)とも呼ばれます)などと呼ばれています。
この項はベイズの定理を考える上であまり活躍することはありません(これはベイズの定理を学んで行けば行くほどはっきりしてきます)。
さて、では実際に計算をしてみましょう。まずは事前確率Pr(A)です。「壷aと壷bが選ばれる割合は1:2」なのでした。そのため、壷aが選ばれるのは1 / 3ですね。
次に尤度Pr(R | A)を考えてみましょう。尤度は「壷aが選ばれた」ときに「取り出した玉が赤玉であることはどのくらいもっともらしいか」を示しているのでした。ちょっと何言っているかよく分かりませんよね。
こういうときは極端な例を考えてみるのも一つの手です。例えば、壷aの中には赤玉しか入っていないとしましょう。赤玉の数はなんでもよいです。流れに沿って赤玉が5個入っていると考えてもよいかもしれません。このとき、尤度はどうなるでしょう?
尤度は「壷aが選ばれた」ときに「取り出した玉が赤玉であることはどのくらいもっともらしいか」と考えるのでした。今は極端な例を考えているので、「壷aが選ばれる」と「取り出す玉は全て赤玉である!」と考えられます。先ほどの書き方をすれば、「壷aが選ばれた」ときに「取り出した玉が赤玉であることはたいへんもっともらしい」と考えておきましょう。
では壷aの中には白玉しか入っていない例を考えてみましょう。この場合の白玉の数もなんでもよいです。もちろん、白玉が5個は言っていると考えてもよいです。
先ほどと同様に考えると、「壷aが選ばれた」とき「取り出す玉は全て白玉である!」と考えられます。こちらは「壷aが選ばれた」ときに「取り出した玉が赤玉であることは全くもってもっともらしくない!」とでも考えられます。
どうでしょう、少し「もっともだ」という言い回しに慣れたでしょうか?
では元の例を考えましょう。「壷aには赤玉が3個、白玉が2個入っている」例です。
「壷aを選ばれた」ときに「取り出した玉が赤玉であることはどのくらいもっともらしいか」。5個のうち3個は赤玉ですから、「60%の確率でもっともらしい」と言うことができそうです。初めのうちはもっとカジュアルに、「赤玉を取り出すのは6割くらいでもっともらしいかな...」などと表現してもよいかもしれません。
さて、事前分布、尤度に関して考えました。後は規格化定数について考えれば計算できますね...ですがこのエントリではここでPr(B | R)について考えます。つまり、「どちらかの壷から玉1個を取り出したとき、それが赤玉であった。そのとき、その赤玉が壷bから選ばれている確率は如何ほどか」を考えます。もう少しお付き合いください。
Pr(B | R)はベイズの定理を用いて、以下のように計算できます。
AとBが入れ替わっただけです。簡単ですね。あまり重要視されていなかった分母に至っては変わりません。
事前分布Pr(B)や尤度Pr(R | B)も簡単に求めることができます。
事前分布も尤度も66.7%でした(これらが同じ値であることに意味はありません)。ここでまた例え話をしましょう。「壷aと壷bは同じ割合で選ばれる」としたらどうでしょう? 「どちらかの壷から玉1個を取り出したとき、それが赤玉であった」ら、「赤玉が出る尤度がより高い壷から取り出したのでは?」と考えるのがより自然ではないでしょうか? 尤度という概念を導入すると、何か人間の自然な思考で考えられるような気がしてきませんか。
もちろん尤度だけでは判断できません。例えば壷aが全て赤玉で詰まっていたとしても、つまり最高潮にもっともらしかったとしても、「壷aと壷bが選ばれる割合は1:100」であったら、壷aから玉自体が取り出される機会がそのものが大変少ないのです。
尤度と事前分布の積を取るのはそのためです。言い換えれば、「ところで全体からするとどの割合でもっともらしいの?」という問いの答えなんですね。ちょっと条件付き確率の考え方に似てますね。
では早速計算してみましょう。
積の値を比較してみたいですね。そのような場合は通分すればよいのでした。
その結果、9:20の割合であることが分かりました。
この割合は「どちらかの壷から玉1個を取り出したら赤玉であった」ときに「選んだ壷が壷aであったか壷bであったか」の割合でもあります。そのため、元の問題の答えもここから計算することができます。
少しトリッキーでしたが、規格化定数を使うことなしに事後確率、つまり「どちらかの壷から玉1個を取り出したら赤玉であった」ときに「その赤玉が壷aから選ばれている」確率を計算することができました。
どうやら事後確率は31%であるようです。事前確率、つまり「事前に情報がないときに壷を選んだら壷aである」確率が33.3%であったので、「どちらかの壷から玉1個を取り出したら赤玉であった」ことで確率が少し落ちたことになります。元々壷bから玉を取り出すことが多い上に、壷bのほうが少し尤度が高かったため、壷aから赤玉を取り出す確率が下がった、と考えるとよいでしょう。
まとめ
今回はベイズの定理について記載しました。
ほとんどの書籍は式を導いた後にすぐ具体的な例を解説していますが、自分にとっては
- どの項が「事後確率」「事前確率」「尤度」「規格化定数」を表しているか
- またその項の重要度はどうか
を頭に入れた後に読み進めたほうが理解が早かったように記憶しています。
どの分野を学ぶときもそうですが、自分の得意なパターンや例題を見つけるのは大事なように思います。理解できるなと思えるのは本当に嬉しいことですし、効率も上がるように思います。
- 作者: 涌井良幸
- 出版社/メーカー: 日本実業出版社
- 発売日: 2009/11/19
- メディア: 単行本
- 購入: 10人 クリック: 71回
- この商品を含むブログ (18件) を見る
- 作者: J.アルバート
- 出版社/メーカー: 丸善出版
- 発売日: 2012/04/05
- メディア: 単行本（ソフトカバー）
- この商品を含むブログ (1件) を見る
SRM 611 Div II
SRM 611 Div II
http://community.topcoder.com/stat?c=round_overview&er=5&rd=15844
SRM 611 Div Iに参加。戦績は惨憺たるものだった。
SRM終了後に@nico_shindanninさんのTopCoderでプログラムしてみた 第2060回(SRM611 直後放送)を視聴。そこで解説されていた素因数、最小公倍数(LCM)の解説があまりにも明瞭で素晴らしかったため、氏の許可を得てブログにまとめることにした。
さて、Div I 250とDiv II 500は類似した問題であった。Div I 250はDiv II 500の知見を利用すると効率よく解くことができる。そのため、このエントリでは生放送での解説から得られる知見を解説し、Div II 500を解くことを目標とする。
LCMSetEasy
http://community.topcoder.com/stat?c=problem_statement&pm=13040&rd=15844
- 問題について考える前にまず、約数、最大公約数、最小公倍数について考察してみよう
- 生放送で氏はこれを以下のように図示した
- 続いて24という数字も考えてみよう。24を素因数分解すると以下となる
- これは以下のように図示することができる
- さて、この数字は600を割り切る。氏の手法を用いた場合、このように図示することができる
- ものすごく分かりやすい...次に70という数字を考えてみよう
- 70は600を割り切らない。これは以下のように図示することができる。7がはみ出るのだ
- 以下は18を考えた例だ。18も600を割り切らないので、はみ出る
- 逆にいえば、はみ出る要素がなければその値は約数であると言えるし、はみ出てしまうのであれば約数ではないと言えることができる
- 次に最大公約数を考えてみよう。計算機を使って最大公約数を求める場合はユークリッドの互除法を用いるのが一般的だが、氏の手法でも考えてみよう
- 24と15の最大公約数を考えてみよう。これははそれぞれ、以下のように図示することができる
- 最大公約数を考えるには重なった部分のみ抽出すればよい。論理積を取る感覚だ
- もう一例挙げておこう。24と900の最大公約数は12だ
- 同様に、複数の数値の最大公約数も論理積を取る感覚で考えることができる
- 15、24、900の最大公約数は3だ
- この図を眺めていると、複数の数値の最大公約数は以下のように考えることができることが分かる
- また、式を以下のように変形していくと、端から数値を処理すればよいということが分かる
- これも図からも直感的に得ることができる
- さて、最小公倍数も同じように図示することが出来る。先ほどとは違い、論理和を取る感覚だ
- 最小公倍数を算出する以下の式と見比べてみると楽しい
- 複数の数値の最小公倍数も論理和を取る感覚で考えることができる
- 例えば、24、105、900の最小公倍数は1800だ
- 最大公約数で考えたときと同様に、複数の数値の最小公倍数は以下のように考えることが出来る
- 式変形を行うと、最小公倍数に関しても端から数値を処理すればよいということが分かる
- ここまでの知識をもとに、SRMの問題を解いてみよう。ある整数xと整数の集合Sが与えられるのであった。例として、xを600、集合Sの数値をS = {8, 12, 15, 18, 105, 300}とする
- 集合Sの部分集合の最小公倍数が数値xとなるか判定するにはどうすればよいだろうか?
- 一先ず、全ての数のLCMを取ってみよう
- ...なんとなく600に近い形にはなっているように見えるが、3が一つ、また7がはみ出してしまっている
- はみ出してしまうのは約数ではないからであった。なので、これらの値は除外してしまってよい。これは、値が整数xを割り切るかどうかで判断できる
- この例の場合、18、105は除外してしまってよい
- 改めて、はみ出さなかった数のみ、つまり約数のみの最小公倍数を考えてみよう
- 最小公倍数の数は600となった!
- そのため、以下の部分集合は最小公倍数がxとなるものの一つであると言える
- 最小公倍数が600とならない場合も考えてみよう。先ほどの集合から300を抜いた場合を考えてみる
- 先ほどと同様にはみ出す、つまり約数ではない数を除外して、最小公倍数を求めよう。その値は300となり、600には足りないことが分かる
これを実装すると以下のようになる(システムテストをパスした実装)。
class LCMSetEasy { public: std::string include(std::vector<int> S, int x) { long long lcm = 1; for (const auto& i : S) if (x % i == 0) lcm = lcm * i / std::__gcd<long long>(lcm, i); return lcm == x ? "Possible" : "Impossible"; }; };
まとめ
- この図示手法は大変簡潔で分かりやすい。最大公約数や最小公倍数を考える際、有用なアプローチの一つに成りえる
- 氏は解説の文書化を快諾くださった。筆者の力量不足であまりよい文章に落とし込むことが出来なかったが、感謝したい
おまけ
@kojingharangさんが大変有用なつぶやきをされていたので抜粋させていただく。
Understanding How cbind/rbind works in R
Rを使うようになってきた最大の動機はR自身で柔軟にデータを作成できることにあります。例えば、パッケージユーザーのための機械学習(1)：決定木 - 銀座で働くデータサイエンティストのブログでは以下のようなXORパターンのデータを用いて決定木の解説を行っています。
「少しランダムが入った市松模様」とでもいうのでしょうか。上述のものは「シンプルなXORパターン」と称されていました。少し点群の数を増やすと市松模様を基調としている様子がはっきりしてきます。
自分も是非、こういうデータを短時間でゴリゴリ作れるようになりたいと思いました。先のエントリでは、ありがたいことにデータを作成した際のRの実装を掲載してくださっていました。
# シンプルなXORパターン
> p11<-cbind(rnorm(n=25,mean=1,sd=0.5),rnorm(n=25,mean=1,sd=0.5))
> p12<-cbind(rnorm(n=25,mean=-1,sd=0.5),rnorm(n=25,mean=1,sd=0.5))
> p13<-cbind(rnorm(n=25,mean=-1,sd=0.5),rnorm(n=25,mean=-1,sd=0.5))
> p14<-cbind(rnorm(n=25,mean=1,sd=0.5),rnorm(n=25,mean=-1,sd=0.5))
> t<-as.factor(c(rep(0,50),rep(1,50)))
> d1<-as.data.frame(cbind(rbind(p11,p13,p12,p14),t))
> names(d1)<-c("x","y","label")...これで出来るのか、という印象。初見では何をやっているか全く検討もつきませんでした。
部分的に実行してコードリーディングを行った結果、cbindとrbindという関数が何回も登場することに気付きました。
> help(cbind)
Combine R Objects by Rows or Columns
Description:
Take a sequence of vector, matrix or data frames arguments and
combine by _c_olumns or _r_ows, respectively. These are generic
functions with methods for other R classes.
Usage:
cbind(..., deparse.level = 1)
rbind(..., deparse.level = 1)
:なるほど、cbindもrbindもベクトルや行列を結合する関数のようです。cとrはそれぞれ、ColumnとRowを表すとのこと。要素数の少ないベクトルや行列で実験を行った結果、なかなか一貫性のある設計であることに気付きました。
まずベクトルのcbindとrbindから考えてみることにします。ベクトルaとbを以下のように置きます。
Rでもベクトルを定義しましょう。
> a = c(1, 2, 3) > b = c(4, 5, 6) > a [1] 1 2 3 > b [1] 4 5 6
これらのベクトルをcbindしたものは行列となります。これをMcとします。Mcは以下のような、2つのベクトルを結合した行列となります。
> cbind(a, b)
a b
[1,] 1 4
[2,] 2 5
[3,] 3 6同様に、ベクトルをrbindしたものをMrとします。Mrの結果は以下のようになりました。
> rbind(a, b) [,1] [,2] [,3] a 1 2 3 b 4 5 6
次に行列のcbind/rbindを考えてみましょう。以下のような行列SとTを考えます。
Rでも行列を定義しましょう。
> S = matrix(c(1, 4, 2, 5, 3, 6), 2, 3)
> T = matrix(c(7, 10, 8, 11, 9, 12), 2, 3)
> S
[,1] [,2] [,3]
[1,] 1 2 3
[2,] 4 5 6
> T
[,1] [,2] [,3]
[1,] 7 8 9
[2,] 10 11 12ベクトルの場合と同様に、行列をcbindしたものをMcとします。Mcは以下のような行列となります。
> cbind(S, T)
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1 2 3 7 8 9
[2,] 4 5 6 10 11 12行列をrbindしたものをMrとします。Mrは以下のような行列となります。
> rbind(S, T)
[,1] [,2] [,3]
[1,] 1 2 3
[2,] 4 5 6
[3,] 7 8 9
[4,] 10 11 12さて、どうでしょう? とても一貫していますね。引数がベクトルである場合でも行列である場合でもcbind(a, b)は以下のような行列を作成するのです。
同様にrbind(a, b)は以下のような行列を作成します。
ベクトルのcbind/rbindには少し注意が必要です - ベクトルが列ベクトルとして扱われるか行ベクトルとして扱われるかは関数によって異なります。また、cbindで以下のように合成することは出来ません(このような結合を行うにはc(a, b)とします)。
さて、ここまでのcbind/rbindの知識を元に、今一度XORパターンの実装を読んでみましょう(再掲)。
# シンプルなXORパターン
> p11<-cbind(rnorm(n=25,mean=1,sd=0.5),rnorm(n=25,mean=1,sd=0.5))
> p12<-cbind(rnorm(n=25,mean=-1,sd=0.5),rnorm(n=25,mean=1,sd=0.5))
> p13<-cbind(rnorm(n=25,mean=-1,sd=0.5),rnorm(n=25,mean=-1,sd=0.5))
> p14<-cbind(rnorm(n=25,mean=1,sd=0.5),rnorm(n=25,mean=-1,sd=0.5))
> t<-as.factor(c(rep(0,50),rep(1,50)))
> d1<-as.data.frame(cbind(rbind(p11,p13,p12,p14),t))
> names(d1)<-c("x","y","label")さて、もう読めるようになりましたね...と簡単に行けばよいのですが、やはりまだ難解ですね。
この実装では結合した行列から最終的にデータフレームを生成していますが、行列を作成することを目標としてみましょう。as.data.frame(...)の中身を作成するところまで、と考えてください。
少し分解しましょう。まず、p11〜p14までのデータを考えてみます。これらはX座標、Y座標共に標準偏差0.5の正規分布に従ったランダムな点群です。X座標、Y座標のそれぞれのベクトルとして用意し、cbindを用いて合成をしています。以下のようなイメージです。
meanを変えることで分布の中心を変えながら、p12、p13、p14を同様に作成します。二つ目の添字が象限を表しています。例えば、p13は第三象限に分布した点群です。そのため、ほとんどの値は負の値となります。
さて、ベクトルを合成した結果は行列となるのでした。これをrbindで合成します。市松模様とするために少し工夫を施し、p11、p13、p12、p14の順番で合成しましょう。
これとは別に0と1のラベルからなるベクトルを作成します。先ほど少し触れましたが、ベクトルとベクトルを直列に合成するにはc()を用います。
最後にcbindを用いて行列とベクトルを合成します。
これで目的の行列を作成することができました。データフレームに変換しなくても行列の段階でプロットすることも可能です。
plot(M[,1:2], col=M[,3], pch=16, xlab='x', ylab='y')
まとめ
ベクトルや行列等を合成して新たな行列を作成する場合はcbind/rbindを使うと便利です。Rにはまだ慣れていないためどんな場合でもデータフレームにまとめようとする傾向にありますが、今後はcbind/rbindをどんどん活用していきたいと考えています。
もちろん好みによるとは思いますが、個人的にはmatrixを使うよりもrbindを使ったほうが直感的に行列が記述出来るよう感じています。
> S = rbind(c(1, 2, 3), c(4, 5, 6))
> S
[,1] [,2] [,3]
[1,] 1 2 3
[2,] 4 5 6以下にmatrixを使った行列の作成を再掲します。是非比較してみてください。
> S = matrix(c(1, 4, 2, 5, 3, 6), 2, 3)
How matrix and vector are multiplied in R?
構文的な一貫性をなかなか感じられないのでRを敬遠していました。しかし時勢には敵わず、使用頻度が増えてきました。
どうせやるのであれば手が増えたほうがよいので、スニペットを多数作成し検証しています。気付いたことをブログに少しずつまとめようかと思っています。
さて、今日はvector型とmatrix型のオブジェクトの乗算を検証しました。行列とベクトルの演算は空気のように扱えないとストレスが貯まるので、どの処理系を学ぶときにも比較的早い段階で検証することにしています。
例えばベクトルの向きを調べることは重要です。vector型が列ベクトルとして解釈されるのであるか、行ベクトルとして解釈されるかで話は大きく違ってきます。
具体例を挙げてみましょう。3 × 3の行列をM、要素数3のベクトルをvとします。行列Mとベクトルvを乗じたとき、ベクトルが列ベクトルであれば、結果も列ベクトルになります。列ベクトルは行列の右側からかけます。
この演算方法を基調とする処理系をPre-multiplyな処理系などと呼びます。右辺がベクトルであることを簡潔に示すためにv'を導入しました。v'は以下のように計算されます。
これに対しベクトルが行ベクトルであれば行列の左側からかけないといけませんし、結果も行ベクトルになります。
この演算方法を基調としている処理系をPost-multiplyな処理系と呼びます。v'の要素は以下のように計算します。そのまま記述すると水平方向に長くなってしまうため右辺のベクトルを転置していることに注意してください。
Pre-multiplyな処理系とPost-multiplyな処理系で同様の結果を得るにはベクトルの方向だけでなく、行列を転置しなければいけないことに注意してください。
例えば同時座標系での移動行列を考えてみましょう。簡単に、移動行列とはあるベクトルの要素を平行移動する行列である、と考えてください。まずPre-multiplyな処理系での移動行列を考えてみることにします。
この行列は二次元ベクトル(x, y)を(Tx, Ty)だけ平行移動します。
さて、移動行列を適用することによって(x, y)が(x + Tx, y + Ty)に射影されました。これをPost-multiplyな処理系でもやってみましょう。Post-multiplyな処理系では行ベクトルを用いますが、要素の結果は等しくならなければなりません。
試しにPre-multiplyな処理系の移動行列をPost-multiplyな処理系で使ってみましょう。
何かめちゃくちゃなことになってしまいました。Pre-multiplyな処理系の行列をPost-multiplyな処理系で用いるには転置をしなければならなかったのです。
今回はうまくいっているようですね。つまり、Post-multiplyな処理系での移動行列は以下のようなものなのでした。Pre-multiplyな処理系のものと見比べて転置されていることを確認するのもよいでしょう。
ここまでで書いたように、処理系でベクトルが列ベクトルで扱われるか行ベクトルで扱われるかで行列の意味合いが全く異ります。そのため、ベクトルが列ベクトルなのか行ベクトルなのか予め調べておくことはとても重要です。言葉を変えれば、処理系でベクトルが列ベクトルで扱われるか行ベクトルで扱われるかを調べることによって、その処理系がPre-multiplyを基調として処理系かPost-multiplyを基調とした処理系かを見極めることができます。
さて、Rのvector型オブジェクトは以下のように記述します。
> v <- c(1, 2, 3) > v [1] 1 2 3
同様に、Rのmatrix型オブジェクトは以下のように記述します。
> M <- matrix(c(1, 0, 0, 0, 1, 0, 3, 5, 1), 3, 3)
> M
[,1] [,2] [,3]
[1,] 1 0 3
[2,] 0 1 5
[3,] 0 0 1Rはvector型とmatrix型の*演算子による乗算を許していますが、これは望む乗算ではありませんでした。調べてみたところ、%*%演算子を使えば望む行列ベクトル演算をするようです。早速やってみましょう。Rの処理系はPre-multiplyな処理系であると仮定し、(x, y) = (0, 0)を(Tx, Ty) = (3, 5)で移動してみましょう。
> M
[,1] [,2] [,3]
[1,] 1 0 3
[2,] 0 1 5
[3,] 0 0 1
> v
[1] 0 0 1
> M %*% v
[,1]
[1,] 3
[2,] 5
[3,] 1数式で表すと以下のようになります。
さて、うまく移動出来たようですね。どうやらRの行列は列ベクトルを基調としていると言えそうです。また、戻り値がmatrix型であることにも注意しましょう。
念のため、Post-multiplyな処理系での演算も試してみましょう。数式で表すと以下です。
> M <- t(M)
> M
[,1] [,2] [,3]
[1,] 1 0 0
[2,] 0 1 0
[3,] 3 5 1
> v %*% M
[,1] [,2] [,3]
[1,] 3 5 1驚くことに、これもPost-multiply系であるかのように動作してしまいました。この結果から考えるに、Rで行列ベクトル演算を行う際、ベクトルは以下のように扱われることが分かります。
- 行列 %*% ベクトル
- ベクトルは列ベクトルであると解釈される
- ベクトル %*% 行列
- ベクトルは行ベクトルであると解釈される
とても柔軟な処理であるようにも見えますが少し注意が必要だな、と改めて感じました。
まとめ
- 行列とベクトルの乗算を行うには%*%演算子を用いる
- Rでのベクトルは文脈によって行ベクトルとも列ベクトルとも扱われる
- RはPre-multiply前提の処理系でもPost-multiply前提の処理系でもなさそうだ
ベクトルの向きは、例えばLaTeXでベクトルを記述する際にも影響します。結果の通りRの処理系は予想以上に柔軟であり決定打には欠けますが、今後ベクトルをLaTeXで表記するときは列ベクトルを基調とし、行列演算はPre-multiplyな処理系を前提にしようと考えています。その理由はmatrixオブジェクトを生成するのに列ベクトルを並べるように要素を記述しなければならなかったことにあります。もっとも、Rの言語全体の一貫性のなさはよく話題になるようですので、最終的には列ベクトル + Pre-multiplyな処理系であるように考えるのが好みであることも少なからず影響しています。
