なめるな!plot

Rの勉強会Tokyo.Rで8回目の参加で初めて発表してきました。
お題はplot。誰もが利用する基本的な関数なんだけれども、意外と知られていない機能や、
業務で利用したときに躓いた点についてまとめました。

なめるな!plot

View more presentations from itoyan110.

今月末はJapan.Rという日本全国のRコミュニティメンバーが参加するイベントが
開かれるようなので、非常に楽しみです。

gsubの使い方メモ

今回のスクリプトは文字列ベクトルの先頭の大文字をとってくるという
ことである。中には空白があったりするので、そのための処理を行う必要がある。
Rでも正規表現を用いることができるということで試してみた。

strs <- c("Hello World","Programming C","LISt Programming","isUpper")
# これで大文字でない箇所を空文字に置き換える
strs.abbr <- gsub("[^A-Z]","",strs) 
strs.abbr
# 先頭の大文字だけを抽出する
strs.split <- strsplit(strs," ");
head.concat <- function(x){
  head.chars <- substring(x,1,1) # 単語の先頭の文字群を得る
  head.upper <- gsub("[^A-Z]","",head.chars) # 大文字だけにする 
  paste(head.upper,collapse="") # 大文字のベクトルを結合
}
strs.abbr.head <- lapply(strs.split,head.concat)
strs.abbr.head <- unlist(strs.abbr.head)
strs.abbr.head

実行結果は以下のようになる。2番目の例では、途中の大文字は消去されている
ことが確認できる。

> strs.abbr
[1] "HW"   "PC"   "LISP" "U" 
> strs.abbr.head
[1] "HW" "PC" "LP" ""  

グループごとに複数の統計量を返す方法の検討

以前紹介したaggregate関数はグループごとにスカラー値を返す
関数をオプション(FUN)に指定するが、分位点を求めるquantile関数のように、
複数の値を返す方法はないかと尋ねられたので、以下のコードを作成した。

set.seed(1234)
let <- sample(letters[1:2],100,rep=T) # lettersはa-zが入っている
letN <- paste(let,c(rep(1,50),rep(2,50)),sep="")
dat <- iris[1:100,1] # アヤメデータから借用

d <- data.frame(let,letN,dat) # データの作成

dsp <- split(d,list(paste(d$let,d$letN,sep="_")))
res <- lapply(dsp,function(x){return(quantile(x$dat))})

dはダミーデータである。
dspはdをsplitしたものであるが、list(paste(d$let,d$letN,sep="_"))
ごとに区切るようにしている。この表現は何をキーにしてリストごとに分けるか
ということを指定しているのであるが、list(d$let,d$letN)では、
直積でリストを作成してしまうので、それを避けるために特殊な切り方にしている。

例えば1列目にはa,bという要素があり、2列目にはa1,a2,b1,b2という要素があるので、
リストは1,2列目の要素が
(a,a1),(a,a2),(a,b1),(a,b2),(b,a1),(b,a2),(b,b1),(b,b2)
となるように分けようとする。

しかしながら実際に1列目がa,2列目がb1という要素が無くても、
空のデータフレームからなるリストを作成しようとするので、それを嫌がって、
空のデータフレームが生成されるのを防いでいるのである。

最後の行にあるlapplyでdspに適用する関数は、無名関数と呼ばれるものを利用している。
通常の関数はmyfunc <- function(x) { ... } という形式が一般的であるが、
myfuncという関数名に相当する部分がなく、関数の定義のみで
記述されていることから無名関数と呼ばれている。
機能を何か所も使いまわすことがなく、かつ短い記述で済む場合に用いると
便利だと思われる機能である。

カーネル平滑化のメモ

今回はカーネル平滑化による曲線フィッティングについてである。カーネル平滑化というのは、簡単に言うと、x軸の各値で、周囲の点にガウス分布の形をした重みをかけて、y軸の値を推定するノンパラメトリックな平滑化手法である。

重みというのは、例えばx=3.2の点を見ており、周囲にx=1,3,4のデータがあれば、データに対する重みは x=3のデータ > x=4のデータ > x=1のデータの順にするということである。ガウス分布の裾を大きくすると、x軸のどの値でもデータ全体を見るようになるので、次第にx軸に平行な平滑化の直線が得られることが予想される。KernSmoothというライブラリのksmoothでカーネル平滑化は実行できるので、example(ksmooth)に手を加えて、それを確かめてみよう。

例に用いられているデータcarsは、車の速度と車が減速してから停止するまでに移動した距離の関係を表したものである。

library(KernSmooth)
require(graphics)

colset <- 2:6
bandset <- c(2,5,10,20,100)

# for文を書かないで1要素ずつ処理
myfunc <- function(bc,xx,yy,type="normal"){
 pp <- ksmooth(xx,yy,type,bandwidth=bc[1])
 lines(pp,,col=bc[2],lwd=2)
 return(pp)
}

plot(cars$speed,cars$dist);
targetY <- 50
abline(h=targetY)
bc <- cbind(bandset,colset)
coords <- apply(bc,1,myfunc,cars$speed,cars$dist)

legend(x=min(cars$speed),y=max(cars$dist),
 legend=paste("bd",bandset,sep=""),pch=19,col=colset)

以上の内容がフィッティングである。次に、このフィッティングカーブを持ちいて、距離から速度を推定することを行う。

ksmooth関数はx軸とy軸の値をセットにして返すので、lines関数で曲線に近い直線を描くことができる。その精度は、オプションで、点数を大きい値に指定することによって実現できる。

y軸の値に最も近いx軸のインデックスを返す関数を以下に示す。実行するときには、上のコードの下に付け加えればよい。

# y軸の値がvalのときのx軸の値を返す
# 初めて超えた点を近似値として返す
retEstIntersection <- function(vec,val){
 tarIdx <- which(vec$y>val)
 if( length(tarIdx) < 1 ){
  return(0)
 }
 idx <- min(which(vec$y>val))
 return(idx)
 # return(which.min(abs(vec$y-val))) # 一番近い点を求める場合
}

# targetYに最も近いデータの行番号が返ってくる
# これをksmoothの返り値のリストに放り込めば，交点のx軸の近似値が返せる
ret <- sapply(coords,retEstIntersection,targetY,simplify=T)
# 上の行はこんな書き方でも代用できる
# ret <- lapply(coords,retEstIntersection,targetY,simplify=T)
# result <- unlist(ret)

以上のコードを実行すると、下記のものが得られる。

1)y軸の値がtargetYを初めて超えたx軸のインデックスを返した結果

> ret
[1] 64 67 70 80  0

ここで0は平滑化で得られた曲線とy=targetYと交わらなかったことを意味する。

2)プロット

フィッティングの精度を見てみると、bandwidthが小さいとデータに引っ張られやすく、大きいとx軸に平行な直線に近いものになり、精度が悪くなる。今回の例では、bandwidth=5もしくは10が一番きれいにフィッティングできているように見える。

applyに渡す関数を自前で書いてみる

Rの高速化のポイントにあげられるものとしてforを使わないということがある。
ベクトル単位で一度に処理せよということである。
今回はベクトル単位で処理する際によく用いられるapplyの挙動について紹介してみる。

まずは肩慣らしに、普通のapplyを使ってみよう。

> mat <- matrix(1:12,nrow=4)
> mat
     [,1] [,2] [,3]
[1,]    1    5    9
[2,]    2    6   10
[3,]    3    7   11
[4,]    4    8   12
> mat[2,3]<- NA
> mat
     [,1] [,2] [,3]
[1,]    1    5    9
[2,]    2    6   NA
[3,]    3    7   11
[4,]    4    8   12
> apply(mat,1,mean)
[1]  5 NA  7  8
> apply(mat,1,mean,na.rm=T)
[1] 5 4 7 8

関数にオプション(na.rm=T)を与えるときには、関数の後ろにオプションを書けばよい。
次に自前の関数を用意してapplyがどのような挙動になるか見てみる。

今回は、あるベクトルと行列に対し、完全に内容が一致する行のインデックスを
取り出す操作を行うことにしよう。

mat <-matrix(1:12,nrow=4)
mat[2,3] <- NA
pattern <- c(3,7,11) # 調べようとする内容
myfunc <- function(m,pat){
 print(c(is.matrix(m),is.vector(m)))
 print(m)
 return( which(all(m==pat)) )
}
apply(mat,1,myfunc,pattern)

自作の関数(myfunc)では、最初の引数mにapplyで一番最初の引数、
2番目の引数patには、FUNの後にある引数が対応している。
上記のスクリプトを実行すると、以下の結果が得られる。

> apply(mat,1,myfunc,pattern)
[1] FALSE  TRUE # print(c(is.matrix(m),is.vector(m)))の内容
[1] 1 5 9 # print(m)の内容
[1] FALSE  TRUE
[1]  2  6 NA
[1] FALSE  TRUE
[1]  3  7 11
[1] FALSE  TRUE
[1]  4  8 12
[[1]] # applyの返り値、integer(0)は長さ0の整数型ベクトル
integer(0)

[[2]]
integer(0)

[[3]] # マッチしたので、数値が返ってくる。しかし、3ではなく1になる
[1] 1

[[4]]
integer(0)

1が返ってきたという上記の内容から、実は1行ずつ判定しているということがわかる。
しかも関数の中で、最初の引数であった行列はベクトルにdropしているのである。
だから、3番目が一致しているということを判別したければ、さらに上記のコードに

拡張を行う必要がある。

res <- apply(mat,1,myfunc,pattern) # 上のコードの最終行
ret <- lapply(res,function(x) length(x)>0) # 対応する行は長さが1以上
idx <- which(unlist(ret))

とすると、idxには3が入って、めでたく必要な情報が得られたことになる。
今回はマトリックスという列の属性がすべて共通なものであったため、
容易に比較ができたが、データフレーム（列ごとに属性が違う）の場合は、
比較対象の属性の中で一番表現範囲の広い型にしてしまってから比較を
行わなければならない。
大雑把ではあるが、文字列型>実数型>整数型>論理型の順で表現できる
範囲が広くなると考えてよい。 入力がどうなるか判断がつかない場合は、
とりあえず文字列型に変換すればよい。 たとえばmatや調べたいパターンが
データフレームであった場合、以下のような修正が考えられるであろう。

myfunc <- function(m,pat){
 mm <- as.character(unlist(m))
 patpat <- as.character(unlist(pat))
 return( which(all(mm==patpat)) )
}

追記：上ではapplyの理解のためにmyfuncの返り値にwhichを付けているが、
実はつけない方が楽にかける。以下のように書けば、データ長で、
行番号の有無を判定する作業が無くて済む。

myfunc <- function(m,pat){
 return( all(m==pat) )
}
ret <- apply(mat,1,myfunc,pattern)
result <- which(ret)

twitteRを使ってみた

RからtwitterとやりとりするtwitteRというパッケージを使ってみました。
Twitter clients for Rを基に操作した忘備録です。

sess <- initSession("ユーザ名","パスワード")

返ってきたsessを用いて、自分のアカウントに関連した操作する。
誰でも見れない部分については、sessを引数にして操作する。

publicTweets <- publicTimeline()

twitter全体で最新のつぶやきを取得。デフォルトでは20とってくる。

korindo_tweets <- userTimeline("korindo")
korindo_tweets
[[1]]
[1] "korindo: 近所の居酒屋でランチなう。「面白そうなのでそろえちゃった☆」とノンアルコールビールが揃っていた。先にビールの種類を増やそうぜ"
(以下略)

あるユーザのタイムラインを取得する。この場合はkorindoさんのタイムラインを
取得している。相変わらずの酒豪ぶりですね。

friendsTweets <- friendsTimeline(sess)

sessに該当する人（自分自身以外にはないでしょうが）の友人の呟きを取ってくる。
当然セッションが必要で、最初に使用したinitSession()を利用する必要が出てくる。

tokyor <- searchTwitter("#TokyoR",num=10)

ある文字列に関するつぶやきを検索。この場合は#TokyoRというハッシュタグを持つ
つぶやきを探してくる。

bgates <- getUser("billgates")

あるユーザを具体的に見たいときは、getUser()を用いる。返り値はユーザ名。
実際に上手く動作するかどうかは、ユーザがプロファイルをロックしているか
どうかに依存している。

his_friends <- userFriends("billgates")

ユーザの友人を取ってくる。

his_followers <- userFollowers("billgates")

ユーザのフォロワーを取ってくる。

つぶやきや、ユーザの友人をとってくることができるので、テキストマイニングの分野や
嗜好分析、レコメンドに使えそうなパッケージですね!

macの設定

外出用にMacBook Proを買ってみた。

UNIXな開発環境をvmware以外で実現してみたかった。
以下個人的なメモ。たまに更新するかも。

3/26

• 右クリックはcontrol+マウスクリック
  
• デフォルトのターミナルはXX-YY-no-macbookと出ているが、これはシステム環境設定の共有で名前を変更すればよい。
  
• gccなどの開発環境が欲しかったので、Xportを入れてみたらexitの挙動が変わった。原因不明
  
• パーティションの変更は「ディスクユーティリティ」からHDのパーティションを選択し、ドライブ名の右下の部分をいじる。
  
• C言語に欠かせないバックスラッシュはoption(alt)+\(￥マーク)でいける
  

ダブルガウシアンの作成

ダブルガウシアンというのは、その名の通り、ガウス分布を二つ重ねたものです。
ピークが2つ登場することに特徴があります。
ダブルガウシアンから値を取得することをRでコーディングしてみました。

ちなみに、単純にN(0,1)+N(1,2)のようにするのはダメ。
再生性よりN(0+1,1+2)=N(1,3)になり、ピークが1箇所になってしまいます。

gauss1 <- rnorm(200,-100,10);# N(-100,10)から200点ランダムに取得
gauss2 <- rnorm(200,200,10);# N(200,10)から200点ランダムに取得

reproductiveGauss <- gauss1 + gauss2; # ガウス分布の再生性よりN(100,20)
# こちらが合成した2つのガウス分布の確率密度関数に相当する
# 100点のみを取得しているが、数がある程度あれば、問題ない
doubleGauss <- sample(c(gauss1,gauss2),100,rep=FALSE);
par(mfrow=c(2,1));
hist(reproductiveGauss,breaks=20);
hist(doubleGauss,breaks=20);

その結果は下図のようになります。結果の違いがすぐに見て取れますね。