この記事はQiita 数学 Advent Calendar 2015の5日目の記事です。（後付け）
http://qiita.com/Akai_Banana/items/b328fe0116d248127a36
にインスパイアされました。
といっても、東ロボに勝とうということではなく、リンク先で解かれていなかった問題をRで解きます。
具体的には、この問題です。
(3)
後日、このクラスでハンドボール投げの記録を取りなおした。次に示したA～Dは、最初にとった記録から今回の記録への変化の分析結果を記述したものである。a～dの各々が今回取りなおしたデータの箱ひげ図となる場合に、⓪～③の組み合わせのうち分析結果と箱ひげ図が矛盾するものは、[カ]、[キ]である。
選択肢：
⓪ A-a
① B-b
② C-c
③ D-d
A : どの生徒の記録も下がった。
B : どの生徒の記録も伸びた。
C : 最初にとったデータで上位1/3に入るすべての生徒の記録が伸びた。
D : 最初にとったデータで上位1/3に入るすべての生徒の記録は伸び、下位1/3に入るすべての生徒の記録は下がった
※図はリンク先を参考にしてください。
道具の準備
and <- function(a,b){
a&b
}
or <- function(a,b){
a|b
}
forall <- function(...){
Reduce(and, ...)
}
exists <- function(...){
Reduce(or, ...)
}
createData <- function(value, num){
if(length(value) > 0){
c(rep(value[1], num[1]), createData(value[-1], num[-1]))
}
}
これ、定義した後で思ったんですが、本当はforallは値域、命題１、命題２、みたいな引数を取るのが適切な気がします。まあ、大目に見てください。ちなみに、orやexistsに至っては定義しても使いませんが、大目に見てください。
最後のCreateDataの定義は、問題分の図を簡単に入力するためのおまじないです。
東ロボ君の場合は、命題をきちんと自分で読み込むのですが、今回はそこまでは目指していないので、とりあえず手動で問題のグラフを読み込みます。
データの準備
data <- createData(1:9 * 5 + 2.5, c(1,4,6,11,9,4,3,1,1))
qx <- c(quantile(data))
qa <- c(7.5, 22.5, 22.5, 27.5, 47.5)
qb <- c(12.5, 22.5, 27.5, 32.5, 52.5)
qc <- c(12.5, 17.5, 22.5, 27.5, 42.5)
qd <- c(2.5, 12.5, 22.5, 32.5, 52.5)
qa～qdは、読み込んだデータです。qxはquantileを入れていて、まさに箱髭図で表現している各種の値が入っています。
解く
さあ、解くぞ解くぞ。
Prop1 <- forall(qx >= qa)
Prop2 <- forall(qx <= qb)
Prop3 <- forall(qx[4:5] <= qc[4:5])
Prop4 <- forall(qx[4:5] <= qd[4:5]) & forall(qx[1:2] >= qd[1:2])
東ロボ君と比べるとちょっと出来が悪いシャイなので、解いても答えを教えてくれません。
> Prop1
[1] FALSE
> Prop2
[1] TRUE
> Prop3
[1] FALSE
> Prop4
[1] TRUE
ということで0と2が答えということになるのですが、確かにその通りですね。
なぜこれで解けるのか
一応、まじめな解説もつけてみます。
A:すべての生徒の記録が下がったとします。最小値を取った人を$x$とするとき、一回目の計測値を$f_1(x)$、二回目の計測値を$f_2(x)$とすると、少なくとも元の最小値よりも小さくなりますから$f_1(x)>f_2(x)$です。そのため、最小値は$f_2(x)$以下になります。※$f_2(x)$とは限りません！
最大値についても同様に考えると、最大値を取った人を$y$とすれば、任意の人$z$に対して$f_1(y)\geq f_1(z)$ですから、$f_2(z)$はそれより小さくなるので、任意の$z$について$f_1(y)>f_2(z)$となります。
同じような考察を行えば、四分位点の各値について、一回目の値>二回目の値であることが分かります。これを計算するのが、Prop1に代入した式の右辺ですね。($qx >= qa$は、TRUEやFALSEのリストになり、これに対してandを繰り返し適用しています。)
B:に関しても、全く上と同じ議論が成り立ちます。
C,Dに関しては、どうでしょうか？
Cの場合、二回目の方が最大値が大きいのは明らかですね。では四分位点はどうなのかというと、上位1/3の範囲に収まっている方の四分位点は、やはりさきほどの議論と同様にして二回目の方が大きいはずです。
そこで、AやBでは四分位数全体に適用した命題を、CやDでは最大値と大きい方の四分位点に限定して適用します。
ということを、Rに計算させることで結果が得られるのでした。
まだ納得できない部分がある方へ
このProp1～4の右辺、なぜ等号を含むのでしょうか？
それは、xx2.5/xx7.5という数値に丸めているからです。
もう少し丁寧に言うと、この問題では具体的な分布が与えられておらず、四分位数がどの数の間にあるか、という条件しか与えられていません。そこで、例えば「0と5の間にある」ということを、2.5という数で表わしています。上の式で計算できるのは、必要条件となります。言い換えると、Prop1～Prop4は、主張が矛盾しないとすれば必ず成り立たないといけない式となります。
とはいえ、一般論としては、もう少し丁寧に議論をすることで、本質的に同じ式で議論ができることが分かる…はずです。それは、また時間があれば。
そういえば
全く余談で、高校数学では$a \geq 0$などと書いた時に、$a=0$となる場合が実際にあることを示さないといけない、とかいう話があったことを思い出しました。値域の議論をしているのならばそうなのでしょうが、$\geq$という記号は、基本的には$a>0$ or $a=0$の場合に真となる記号、という意味だと思いたいですね。
他にも他にも…その１：論理
他の問題もRで解こうとしてみます。他はだいたい力づくで解けるのですが、少しプログラミングらしい道具が必要そうな問題を選んで、遊んでみます。
[1]
(1)命題「($p_1$かつ$p_2$)⇒($q_1$かつ$q_2$)」の対偶は[ア]である。
選択肢：
⓪($\bar{p_1}$または$\bar{p_2}$)⇒($\bar{q_1}$または$\bar{q_2}$)
①($\bar{q_1}$または$\bar{q_2}$)⇒($\bar{p_1}$または$\bar{p_2}$)
②($\bar{q_1}$かつ$\bar{q_2}$)⇒($\bar{p_1}$かつ$\bar{p_2}$)
③($\bar{p_1}$かつ$\bar{p_2}$)⇒($\bar{q_1}$かつ$\bar{q_2}$)
turnA <- function(...){
Reduce(function(a,b){a&b}, ...)
}
`%->%` <- function(a,b){!a | b}
forallは、turnAという名前に変えました。下側の演算子は、「ならば」にあたる奴ですね。
q1 <- function(l){
p1<-l[1];p2<-l[2];q1<-l[3];q2<-l[4]
(p1&p2)%->%(q1&q2)
}
f1 <- function(l){
p1<-l[1];p2<-l[2];q1<-l[3];q2<-l[4]
(!p1|!p2)%->%(!q1|!q2)
}
f2 <- function(l){
p1<-l[1];p2<-l[2];q1<-l[3];q2<-l[4]
(!q1|!q2)%->%(!p1|!p2)
}
f3 <- function(l){
p1<-l[1];p2<-l[2];q1<-l[3];q2<-l[4]
(!q1&!q2)%->%(!p1&!p2)
}
f4 <- function(l){
p1<-l[1];p2<-l[2];q1<-l[3];q2<-l[4]
(!p1&!p2)%->%(!q1&!q2)
}
ここまで、問題文を忠実に書き写しました。
li <- list()
for(i in 1:16){
li[[i]] <- c((i %/% 1) %% 2 == 0
,(i %/% 2) %% 2 == 0
,(i %/% 4) %% 2 == 0
,(i %/% 8) %% 2 == 0)
}
turnA(laply(li, q1)==laply(li, f1))
turnA(laply(li, q1)==laply(li, f2))
turnA(laply(li, q1)==laply(li, f3))
turnA(laply(li, q1)==laply(li, f4))
今回は、前回の失敗を反省して、答えをただちに出力するようにしました。
> turnA(laply(li, q1)==laply(li, f1))
[1] FALSE
> turnA(laply(li, q1)==laply(li, f2))
[1] TRUE
> turnA(laply(li, q1)==laply(li, f3))
[1] FALSE
> turnA(laply(li, q1)==laply(li, f4))
[1] FALSE
他にも他にも…その２：素因数分解
もはや問題も省略しますが、$a=756$を素因数分解する場合。
factorize <- function(x){
if(x==1){return(1)}
p<-getPrime(x,primes)
c(p, factorize(x/p))
}
getPrime <- function(x,p){
if(x==1){return(1)}
if(x %% p[1] == 0){return(p[1])}
getPrime(x,p[-1])
}
primes <- c(2,3,5,7,11)
factorize(756)
> factorize(756)
[1] 2 2 3 3 3 7 1
英語のセンスのなさがちょっと恥ずかしいですね。
あとの問題は、割とそのまま頑張っても解けるのでは、という気がしています。
primesは、素数の無限リストみたいなものを渡すのが本当はよいのかなあと思います。この問題を解くだけであれば、2,3,7でも十分です。。。