2014年6月24日火曜日

Rの論理演算子 & と &&、 | と || の違い

【1つのもの】(「&」と「|」)

 ベクトルのそれぞれの要素を評価

> c(TRUE, TRUE, FALSE, FALSE) & c(TRUE, FALSE, TRUE, FALSE)

[1]  TRUE FALSE FALSE FALSE


【2つのもの】(「&&」と「||」)

 ベクトルの1つ目の要素のみ評価

> c(TRUE, TRUE, FALSE, FALSE) && c(TRUE, FALSE, TRUE, FALSE)

[1] TRUE

前者はRの演算子や関数ではおなじみの挙動ですよね。

この動きについては、ヘルプを見れば書いてあります。呼び出し方は以下↓

> ?"&"
starting httpd help server ... done

訳は私によるアバウトなもの。

& and && indicate logical AND and | and || indicate logical OR.

「&」と「&&」は論理積(AND)を、
「|」と「||」は論理和(OR)を表します。

The shorter form performs elementwise comparisons in much the same way as arithmetic operators.

短い方の形式は、他の算術演算子と同じように、要素ごとの比較を行う。

The longer form evaluates left to right examining only the first element of each vector.


長い方の形式は、左から右へと評価を行い、両者のベクトルの最初の要素のみを比較する。

Evaluation proceeds only until the result is determined.

評価は結果が確定するまでしか行われない。

The longer form is appropriate for programming control-flow and typically preferred in if clauses.

長い方の形式は、プログラムのフロー制御に向いており、if文などに適している。


以後はけっこう蛇足な感じですので、お暇な方だけ・・・

「左から右へと評価」のところですが、こう解釈しました。例えば、

TRUE || FALSE && TRUE

という順序で並んでいた場合、まず「TRUE || FALSE」の部分が評価され、この部分がTRUEだと判明した後に、後続の「 && TRUE」が評価され、全体としてTRUEになる、と。

ただ、疑問なのは、ANDとORだけで構成される論理式において、評価順によって結果が変わるような例が思いつかないんですよね。どの順で評価しても同じ結果になるような気がして、わざわざ「左から右」と評価順を明記する必要があるのだろうかと、ちょっと気になりました。

で、もしかしたら「評価は結果が確定するまでしか行われない」と関連しているんじゃないかと思ったんですね。

左から評価していくけど、結果が確定したら、そのあとは評価しないから、評価の順序は気にしておいて、という意味かなと。

で、試してみました。

f1 <- function(){
  print("in f1")
  return(TRUE)
}
 
f2 <- function(){
  print("in f2")
  return(FALSE)
}
 
if ( f1() || f2() ){
  print("in if clause")
}


if文の条件の評価で、まずf1が呼ばれます。結果はTRUEであり、条件文の後続は||なので、f2の結果を知る必要はないですよね。なので、f2は呼びませんよ、というのがヘルプの意図するところだと。
で、実行してみると、

> if ( f1() || f2() ){
+   print("in if clause")
+ }


[1] "in f1"
[1] "in if clause"


おお、成功です。f2は呼ばれないままif文の中に入っています。

たしか、perlなんかでもこういう仕様になっていて、この仕様を使ってうまい感じでコードが書けたりしたように記憶しています。

ただ、結果がこの仕様に依存するようなコードを書くのは、自分一人専用のちょっとしたスクリプトでない限りは避けた方がいいようにも思います。

ちなみに、短い形式の方を使うと、

> if ( f1() | f2() ){
+   print("in if clause")
+ }


[1] "in f1"
[1] "in f2"
[1] "in if clause"


のように、f1の結果にかかわらず、f2は呼ばれることになります。

やっぱり、この仕様を活用するのは、ちょっと危険な匂いがしますよねえ。




2014年6月20日金曜日

Rのpch指定で使えるプロット文字の種類一覧

Rでplot関数を使った場合、デフォルトではプロット文字は○(白丸)になりますが、それ以外を使いたい場合は、pchオプションに0から25の数字を指定すればOKです。

で、どの数字を指定すればどのプロットキャラクターになるかというのは、こちらにまとめられています。

53. グラフィックスパラメータ(弐) | R-Source

でも、実際にプロットしたときの見え方や印象も確認したいですよね。cex指定との兼ね合いもあるし。

ということで、自分の環境で一気にplotしてみて、見え方を確認してみる、なんてのはどうでしょうか。

いろいろなpch指定とcex指定
sizes <- seq(0.5, 2.0, 0.25)
plot( 0,0, type="n", xlim=c(0,25), ylim=c(0.5,2.0), xlab="", ylab="")
for(i in sizes){
  par(new=T)
  plot( 0:25, rep(i,26), pch=0:25, cex=i,
        ylim=c(0.5, 2.0), xlab="pch", ylab="cex" )
}


なんだか目の検査みたい。あと、pch=11のダビデの星(六芒星)はどう見ても不恰好で、一生使うことはないでしょう。

曲線をプロットした感じは、こんな風になります↓


いろいろなpch指定で曲線をプロット
# プロットする曲線の定義
f <- function(x){ 2 * (sin( x - 0.5*pi ) + 1 ) }
 
x <- seq(0, 2*pi,length=30)
plot( 0,0, type="n", xlim=c(0,2*pi),ylim=c(0,30), xlab="", ylab="")
for(i in 0:25){
  par(new=T)
  plot( x, f(x)+i, pch=i, ylim=c(0,30), ylab="pch" )
}


ちなみに、数字でなく直接文字をダブルクォートで指定してもOK。(例: pch="あ")

むりやり変なサンプルを作ってみる↓


文字(日本語)で曲線をプロット
S <- "文字列を使って曲線をプロットしていくようなこともできますよ"
n <- nchar(str)
x <- seq(0, 2*pi, length=n)
plot( 0,0, type="n", xlab="", ylab="", xlim=c(0,2*pi), ylim=c(-1,1) )
for(i in 1:n){
  par(new=T)
  plot ( x[i], sin(x[i]), pch=substr(str,i,i), xlim=c(0,2*pi), ylim=c(-1,1) )
}


文字をプロットするんだったらtext関数を使った方がすっきり書けるのに、というつっこみはご勘弁を。pchのテーマで始めたもんですから・・・


【おまけ】
pchはPlotting CHaracterの略、cexはCharacter EXxpansionの略だったりします。ほら、覚えやすくなったでしょ。




2014年6月16日月曜日

Rを使って曲線で囲まれた領域に色をつける

曲線のグラフを塗りつぶすやり方を前に書いたんですが↓

Rのsegments関数でグラフを塗りつぶす - Rプログラミングの小ネタ

Rのpolygon関数でグラフを塗りつぶす方法を理解する - Rプログラミングの小ネタ

上記は曲線とx軸に囲まれた領域に着色するというものでしたが、ほとんど同じやり方で曲線と曲線に囲まれた領域の色塗りもできますよね、という話。


curve(sin(x), 0, 4*pi)        # サイン曲線を描画
curve(cos(x), 0, 4*pi, add=T) # コサイン曲線を描画
 
n <- 1000                    # x軸上に取る点の数
x <- seq(pi, 3*pi, length=n) # プロットするx座標を求める
y1 <- sin(x)              # 曲線上のy座標を求める 曲線1
y2 <- cos(x)              # 曲線上のy座標を求める 曲線2
 
# 曲線1上→曲線2上の順で、頂点を指定する
polygon( c(x,rev(x)), c(y1,rev(y2)), col="gray")


上記は曲線上にたくさんの点を置いていき、それらを頂点とする多角形をpolygon関数で描画するというやり方。頂点の列がクロスしても、そのように塗りつぶしてくれるので、どちらの曲線が上にくるとか下にくるとか気にしなくてもOKです。

↓もう一つのやり方


# polygon関数のところ以外は前述の例と同じコードを実行後
 
# 曲線1上の点と曲線2上の点を結ぶ線分で領域を埋め尽くす
segments( x,y1, x,y2, col="gray")


setmentsは、曲線1の上にとった点と曲線2の上にとった点を線分で結ぶ関数なのですが、点を多数とると、この線分で領域が埋め尽くされて、まるで塗りつぶしたようになるというやり方。線分は上から下に引いても、下から上に引いても同じようになるので、こちらも曲線の上下関係を気にする必要はありません。

前者のやり方だと点は100個も取れば滑らかに見えるようになりますが、後者のやり方だと100個では足りないですね。(上記例では1000個にしました)




2014年6月11日水曜日

Rの箱ひげ図(boxplot)で外れ値を表示させない

例えば17人のテストの点数のデータがあって、↓こんなだったとしましょう。

x <- c( 40, 45, 50, 55, 60, 62, 63, 64, 65, 66, 67, 68, 70, 75, 80, 85, 90)


そのまま、boxplotを使うと↓こんな感じになります。

boxplotの表示結果(外れ値扱い 有り)

boxplot(x)


で、
  「この小っちゃい丸は何だ?」
  「外れ値です」
  「外すなよ」

みたいなこともあると思います。

外れ値と言っても、プロットされているわけだから情報が落ちているわけじゃないし、構わないような気もしますが、まあ、特別扱いしたくないって時もあるでしょう。

で、ヒゲをどこまで伸ばすかというのは、boxplot関数のrangeオプションで指定します。デフォルトだと、箱の長さの1.5倍までとなっております。

このrangeをめちゃくちゃ大きくしておけば、事実上、ヒゲは最大のデータ、最小のデータまで伸びることになりますが、そんなことしなくても、特別な値 0 を指定することで、最大・最小まで伸ばすことができます。

boxplotの表示結果(外れ値扱い 無し)

boxplot(x, range=0)


rangeの指定によって、ヒゲの長さがどう変化するか見てみましょう。

boxplotで外れ値扱いの範囲を変化させたもの(外れ値の表示 有り)
par(mfrow=c(1,4))
stripchart(x, vertical=T, pch=1, main="散布図")
boxplot(x, range=1  , main="range=1")
boxplot(x, range=1.5, main="range=1.5(デフォルト)")
boxplot(x, range=2.0, main="range=2.0")


一番左のはデータをそのままプロットした散布図です。

左から二番目はrange=1となっているので、箱の長さの1倍分までヒゲを伸ばします。25パーセンタイルが60点、75パーセンタイルが70点なので、箱のサイズは10。つまり上方には、80(=70+10)までヒゲが伸びることになります。

いつでも80まで伸びるというわけではなく、80以下で最も大きいデータまで伸びます。今回はちょうど80点というデータがあるので、そこまで伸びているというわけです。■

[2014年6月12日追記]「外れ値として扱うが表示はしない」という風にしたい場合は、 outline=F というオプションを使います。


boxplotで外れ値扱いの範囲を変化させたもの(外れ値の表示 無し)
x <- c( 40, 45, 50, 55, 60, 62, 63, 64, 65, 66, 67, 68, 70, 75, 80, 85, 90)
par(mfrow=c(1,4))
stripchart(x, vertical=T, pch=1, ylim=c(40,90), main="散布図")
boxplot(x, range=1  , outline=F, ylim=c(40,90), main="range=1\noutline=F")
boxplot(x, range=1.5, outline=F, ylim=c(40,90), main="range=1.5(デフォルト)\noutline=F")
boxplot(x, range=2.0, outline=F, ylim=c(40,90), main="range=2.0\noutline=F")





2014年6月8日日曜日

Rで複数のベクトルからランダムに1つを選択する

読者からの質問に答えるコーナーを作ろうと思いましたが、そんなものは届かないので、Yahoo知恵袋に寄せられた質問に勝手に答えることにしてみました。

統計ソフトRにおいて『ベクトル(1,3),(2,6),(3,10)の中からランダムに一つとり、そ... - Yahoo!知恵袋

統計ソフトRにおいて『ベクトル(1,3),(2,6),(3,10)の中からランダムに一つとり、そのベクトルをxとする』ということをしたいのですが、可能でしょうか?

すでにベストアンサーが選ばれており、↓こんなやり方が提示されておりました。

a = matrix(c(1,3,2,6,3,10),2,3)
x = a[,floor(runif (1,min=1,max=4))]


このスクリプトは、質問者の要望どおりに動きます。

では、私は何をすれば・・・?

(心の声: ケチをつければいいんじゃよ)

そうします。

回答者は、選択肢となる複数のベクトルが列ベクトルになっているような行列を作って、その中からランダムに列を選択するという方法をとっています。

でも、もっと単純に、複数のベクトルから成るリストを作って、その要素から選択するという方が、理解しやすいように思います。

また、runifという連続型の確率分布の関数を使って、その結果をfloor関数で離散値(整数)に変えています。最初から離散型のsample関数を使えばいいのでは?

あと、代入に「 = 」(イコール)を使っていますが、推奨されている 「 <- 」 を使ったほうがいいかなあって。これは、「左辺と右辺は等しい」のような意味ではなく、「左辺に右辺を代入する」って意味だぞってことを明確にするためです。

で、私のお勧めするスクリプトは下記です↓

a <- list( c(1,3), c(2,6), c(3,10) )
x <- unlist( sample(a, 1) ) # リストaの要素から1つ選ぶ



リストの要素をsample関数で取り出すと、部分リストとして値が返されるみたいなので、unlist関数でベクトルにしてやらなきゃいけないところが、ちょっと一手間ですが、それ以外はごく自然でシンプルだと思います。


【蛇足】
ベストアンサーの回答者の上から目線がちょっと気になりました。
subspace2013さん ,R言語の乱数ね。

(中略)

さて、これが期待されていた実行結果カナ? 何が解らなかったのだろう?





2014年6月5日木曜日

Rのcontour(等高線)を使って陰関数のグラフを描く

Rで、y = f(x) という形の関数のグラフを描きたければ、curve関数を使えばいいですよね。

curve関数で描いた3次曲線のグラフ

curve( x^3 - x, from=-1.5, to=1.5, asp=1)

で、簡単にいきそうにないのは、f(x,y) = 0 という形。陰関数ってやつですね。x^2 + y^2 = 1とかも一例ですね。この例みたいな円の方程式くらいだったら、yについて解いちゃうとか、x=cos(t),y=sin(t)みたいに媒介変数を使う方法もあるんですが、この方法が任意の陰関数について使えるというわけではない。

で、どうやるかというと、答えは↓ここにあります。

陰関数で表された曲線の描画 - 裏 RjpWiki


z = f(x,y) みたいに考えて、xy平面の各格子点に対してzの値を計算する。そして、contour関数を使って等高線を引く。その際、levels=0とオプション指定すれば、f(x,y) = 0のところに線が引かれるという、ちょっとトリッキーなやり方ですが、これで任意の陰関数のグラフが描けるので、とてもうれしいです。

この手法を理解するために、あえて簡単なサンプルをあげてみます。


contour関数で描いた円のグラフ
n <- 10 # メッシュ分割数
x <- seq( -1, 1, length = n )
y <- seq( -1, 1, length = n )
z <- matrix(0, n, n)
for (i in 1:n) {
  for (j in 1:n) {
    z[i, j] <- x[i] ^ 2 + y[j] ^ 2 # ここに陰関数を書く
  }
}
contour(x, y, z, drawlabels = F, levels=1, asp=1)

円がちょっと多角形っぽく見えるのはnが小さいからです。大きくすればなめらかな曲線になります。

この手法、contour関数を使ったことないとピンとこないですよね。

計算されたzの値を見てみると分かりやすいと思います。

> round(z,2)
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,] 2.00 1.60 1.31 1.11 1.01 1.01 1.11 1.31 1.60  2.00
 [2,] 1.60 1.21 0.91 0.72 0.62 0.62 0.72 0.91 1.21  1.60
 [3,] 1.31 0.91 0.62 0.42 0.32 0.32 0.42 0.62 0.91  1.31
 [4,] 1.11 0.72 0.42 0.22 0.12 0.12 0.22 0.42 0.72  1.11
 [5,] 1.01 0.62 0.32 0.12 0.02 0.02 0.12 0.32 0.62  1.01
 [6,] 1.01 0.62 0.32 0.12 0.02 0.02 0.12 0.32 0.62  1.01
 [7,] 1.11 0.72 0.42 0.22 0.12 0.12 0.22 0.42 0.72  1.11
 [8,] 1.31 0.91 0.62 0.42 0.32 0.32 0.42 0.62 0.91  1.31
 [9,] 1.60 1.21 0.91 0.72 0.62 0.62 0.72 0.91 1.21  1.60
[10,] 2.00 1.60 1.31 1.11 1.01 1.01 1.11 1.31 1.60  2.00

↓xy平面上にzの値をプロットしてみました

zの値をxy平面にプロット

plot( x, y, xlim=c(-1,1), ylim=c(-1,1), type="n", asp=1)
for (i in 1:10) {
  for (j in 1:10) {
    text( (i-5.5)/4.5, (j-5.5)/4.5, sprintf("%1.2f", z[i,j]))
  }
}

上記のz値の値に対して、z = 1 のところに等高線を引くと下記のようになり、

zの値のプロットと、z = 1 の等高線(contour)

plot( x, y, xlim=c(-1,1), ylim=c(-1,1), type="n", asp=1)
for (i in 1:10) {
  for (j in 1:10) {
    text( (i-5.5)/4.5, (j-5.5)/4.5, sprintf("%1.2f", z[i,j]))
  }
}
contour(x, y, z, drawlabels = F, levels=1, add=T)

これが、この手法の仕組みとなります。

levelsオプションにベクトルで複数の値を指定すると、より等高線っぽくなります。drawlabels(線が示す値)も表示してみましょう↓

z = 0.1, 0.5, 1.0, 1.5 の等高線(contour)

contour(x, y, z, drawlabels = T, levels=c(0.1, 0.5, 1, 1.5), asp=1)

この手法の仕組み、理解していただけましたでしょうか。これで、いろんな陰関数のグラフを描きまくってください。




2014年6月1日日曜日

Rで正規分布のグラフを書いてビリギャルの偏差値40アップを視覚化

そういえば、偏差値を気にしなくてはいけなかった頃(中高生とか)、私は偏差値の意味を理解してなかったなあ、なんて。

今回の模試は易しい問題が多くて、全受験者の平均点が前回よりも○点高くなった、なんてことはよくある話なので、テストの点数自体は何の指標にもならないですよね。

でも、平均と標準偏差を基準に合わせるように調整すれば、複数回の模試を比較できるじゃん、って発想ですよね。で、偏差値というのは、平均が50、標準偏差が10になるように調整した点数であると。

さて、↓この「学年ビリのギャル」(注:写真は本人ではありません)は、偏差値を40上げたらしいけど、それってどれくらいすごいのか。

学年ビリのギャルが1年で偏差値を40上げて慶應大学に現役合格した話
学年ビリのギャルが1年で偏差値を40上げて慶應大学に現役合格した話

同じ「偏差値40アップ」でも、20→60と、30→70では具合が違います。調べてみると、どうやら約30→約70みたいな感じだったようです。(60じゃ慶應入れませんものね)

偏差値という指標が有効だという発想は、テストの点数が正規分布に近い分布をしているという前提に基づいていますよね。つまり、偏差値に換算したときに、↓こんな感じに分布しているであろうという。

curve( dnorm(x, mean=50, sd=10), from=0, to=100 )



偏差値、30と70のところに線を書いてみます↓


curve( dnorm(x, mean=50, sd=10), from=0, to=100 )
abline( v=30, col="red" )
abline( v=70, col="blue" )


最初、赤線のあたりにいたのが、真ん中の山のあたりにいる有象無象どもを追い抜いて、青線のあたりまで来たという感じです。

横軸が、y = 0 のところにないので、正規分布曲線の下側の領域の面積が分かりにくいですね。ということで、polygon関数でグラフを塗りつぶす方法を理解する - Rプログラミングの小ネタ の方法にしたがって、色をつけてみました↓


n  <- 100
 
xs <- seq(0, 30, length=n)
ys <- dnorm(xs, mean=50, sd=10)
polygon( c(xs,rev(xs)), c(rep(0,n),rev(ys)), col="red")
 
xs <- seq(70, 100, length=n)
ys <- dnorm(xs, mean=50, sd=10)
polygon( c(xs,rev(xs)), c(rep(0,n),rev(ys)), col="blue")


偏差値10以下とか偏差値90以上というのが、ほとんどいない(でも、もちろんゼロでもない)ということが分かりますね。

数値的に知りたい場合は、pnorm関数を使うと分かります。pnormは下側確率を出してくれるので、下の書式では偏差値30以下の割合が分かります。

> pnorm( 30, mean=50, sd=10 )
[1] 0.02275013

つまり、偏差値30以下の人は、全体の2.28%くらいいるということが分かります。

> pnorm( 70, mean=50, sd=10 )
[1] 0.9772499

上記は、偏差値70以下の人が、全体の97.7%くらいを占めるという意味です。

70以上の人の割合が知りたければ、100%から引けばいいので、

> 1 - pnorm( 70, mean=50, sd=10 )
[1] 0.02275013

正規分布は左右対称なので、偏差値30以下の人と、偏差値70以上の人の割合は同じになります。

さて、先述のビリギャルがどのくらい順位をあげたのかというのは、具体的な数字を当てはめると分かりやすいでしょうか。

例えば、全体が10000人だったとして、さきほどのpnormで算出した値を当てはめてみると、最初 9772位だったのが、受験時には228位になった、みたいな感じになります。