scala IDE for Eclipseバージョン上がりましたね。
ピアソンの相関係数を実装 in Scala
レコメンデーション等をするためには各アイテムの相関度などを調べる手法があり、そこで必要なのがユークリッド距離やピアソンの相関係数など。今回実装してみたのはピアソンの方。
２つの変数X,Y（X,Yに属する値のリストみたいなもん。実装ではStream[Double]という形を取っている)がありそれらの相関度をしらべたいときに上の公式を用いて相関係数を求めることが出来る。各偏差（Xi-平均のところ)を求めるところが共通しているからそこはメソッドとしてまとめれそう。それと分母分子両方に1/(n-1)が入っててこれは相関係数だけを求めるなら省くことが出来るけど標準偏差とかも正確に出したいから省略しないでおく。下のは省略した形。
とりあえず適当にデータフレームのクラスを作ってみる。
class Dataset(x:Stream[Double],y:Stream[Double]){}
後々のことを考えるとこのクラスはどうかと思うけどとりあえず２変量のケースにのみ考える。xは例えばあるクラスにおける生徒達の中間試験の結果、yは期末試験の結果。次に各変数の平均値を求める。
def ave(data:Stream[Double]):Double
=data.reduce( (a,b)=>a+b)/data.length
平均を使って各変数のそれぞれの値の偏差を求める。分母で各偏差の２乗も求めなければならないのでさらに各要素の偏差の２乗を求めるメソッドもつくる。
def dev(data:Stream[Double],average:Double)
:Stream[Double]=data.map(x=>x-average)
def devto2(devdata:Stream[Double])
:Stream[Double]=devdata.map(x=>(math.pow(x,2)))
各偏差の２乗の合計の平方根×要素数-1は標準偏差。標準偏差を求めるメソッド。reduceは何かと重宝する。
def devto2(devdata:Stream[Double])
:Stream[Double]=devdata.map(x=>(math.pow(x,2)))
ここまでで分母が完成。次に共分散を求めるために先程作ったXの偏差たちとYの偏差たちをzipしちゃう。
def zipdev(dev1:Stream[Double],dev2:Stream[Double])
:Stream[(Double,Double)]=dev1.zip(dev2)
これ効率的にどうなんだろ。でもまあ後々のこと考えたらこっちのほうが好み。効率優先しすると後々何してんのかわからなくなることあるし。とりあえず可読性の高い書き方してからそれを徐々に効率よくしていく、みたいな。あとは各要素を掛け合わしてやつの合計を求めれば共分散がでる。ピアソンの相関係数は分子/分母で完了。
def cov(zipdata:Stream[(Double,Double)])
:Double=zipdata.map(x=>x._1*x._2)
.reduce( (a,b)=>a+b)/(zipdata.length-1)
def pear(covari:Double,variX:Double,variY:Double)
:Double=covari/(variX*variY)
これが全体。
class Dataset(x:Stream[Double],y:Stream[Double]){
def ave(data:Stream[Double]):Double=data.reduce( (a,b)=>a+b)/data.length
//平均
def dev(data:Stream[Double],average:Double)
:Stream[Double]=data.map(x=>x-average)
//各偏差
def devto2(devdata:Stream[Double])
:Stream[Double]=devdata.map(x=>(math.pow(x,2)))
//各偏差の２乗
def sd(dv2:Stream[Double])=math.sqrt(dv2.reduce( (a,b)=>a+b)/(dv2.length-1))
//√分散＝標準偏差
def zipdev(dev1:Stream[Double],dev2:Stream[Double])
:Stream[(Double,Double)]=dev1.zip(dev2)
//2つのデータセットの偏差をTuple化
def cov(zipdata:Stream[(Double,Double)]):Double
=zipdata.map(x=>x._1*x._2)
.reduce( (a,b)=>a+b)/(zipdata.length-1)
//2つのデータセット共分散(=対応するXとYの偏差の積の平均）
def pear(covari:Double,variX:Double,variY:Double)
:Double=covari/(variX*variY)
//２つのデータセットのピアソン相関関係
val xav=ave(x)
val yav=ave(y)
val xdev1=dev(x,xav)
val ydev1=dev(y,yav)
val xdev2=devto2(xdev1)
val ydev2=devto2(ydev1)
valsdX=sd(xdev2)
valsdY=sd(ydev2)
valzipXY=zipdev(xdev1,ydev1)
valcovXY=cov(zipXY)
valpearson=pear(covXY,sdX,sdY)// peason=covariance / (SD of X * SD of Y)
def summary={
println(" ")
valxsum=Stream("","X; "+x,"mean -> "+xav,"deviation -> "+xdev1,"deviation^2 -> "+xdev2,"standard deviation -> "+sdX)
valysum=Stream("","Y; "+y,"mean -> "+yav,"deviation -> "+ydev1,"deviation^2 -> "+ydev2,"standard deviation -> "+sdY)
valsumxy=Stream("","covariance -> "+covXY,"peason's correlation -> "+pearson)
Stream(xsum,ysum,sumxy).foreach(x=>x.foreach(println))
}
}
Rみたくsummaryメソッドつくってみた。こんな感じで実行してみる。
object UnitTes extends App{
val x=Stream(0.5,3,5,3,3,6,7,8,5,2)
val y=Stream(5.0,3,5,3,3,6,7,8,5,2)
val data=new Dataset(x,y)
data.summary
}
結果。
X; Stream(0.5, 3.0, 5.0, 3.0, 3.0, 6.0, 7.0, 8.0, 5.0, 2.0)
mean -> 4.25
deviation -> Stream(-3.75, -1.25, 0.75, -1.25, -1.25, 1.75, 2.75, 3.75, 0.75, -2.25)
deviation^2 -> Stream(14.0625, 1.5625, 0.5625, 1.5625, 1.5625, 3.0625, 7.5625, 14.0625, 0.5625, 5.0625)
standard deviation -> 2.3481671339342287
Y; Stream(5.0, 3.0, 5.0, 3.0, 3.0, 6.0, 7.0, 8.0, 5.0, 2.0)
mean -> 4.7
deviation -> Stream(0.2999999999999998, -1.7000000000000002, 0.2999999999999998, -1.7000000000000002, -1.7000000000000002, 1.2999999999999998, 2.3, 3.3, 0.2999999999999998, -2.7)
deviation^2 -> Stream(0.0899999999999999, 2.8900000000000006, 0.0899999999999999, 2.8900000000000006, 2.8900000000000006, 1.6899999999999995, 5.289999999999999, 10.889999999999999, 0.0899999999999999, 7.290000000000001)
standard deviation -> 1.9465068427541912
covariance -> 3.638888888888889
peason's correlation -> 0.7961297534564015
時系列解析入門v2 DetrendingとDifferencing
前回の記事では時系列解析の概要、定常性、そして自己相関について触れたから今回はその続きのDetrendingから。そもそもDetrendingとは、時系列データが持っているTrend(傾向)を取り出すことでより定常なデータを入手すること。今回利用するのは回帰分析を利用した単純なアプローチで、時系列データを回帰分析し回帰直線を求め、その後時系列データから回帰直線を差し引いたデータ(残差のデータ)を得る。Trendが除去されている分残差のデータはより定常に近づく。
tsdata=ts(data)
plot(tsdata)
reg=lm(tsdata~time(tsdata))
abline(reg)
acf(auto correlaton function)関数で定常性を比べてみる。
residata=resid(reg)
acf(residata)
１つ目が元々の自己相関。2つ目がDetrending後の自己相関。
明らかに前者より後者の方が定常に近くなっている。ただ、未だ非定常なのでDifferencingによってさらにTrendを除去する。具体的には、すべての値Xtに対してとなりの値Xt-1との差を求めている。
ΔX=Xt-Xt-1
例えばX=[4,2,4,2,6,7,3]だとするとΔX=[-2,2,-2,4,1,-4]となる。もちろんそれぞれ２つの値の差を求めるわけだからΔX全体の要素数はXの要素数の-1となる。この例で言えばlength(X)=length(ΔX)+1となるわけだ。
difdata=diff(residata)
acf(difdata)
この自己相関を見る限りどうやらこの時系列データはLag=1まで自己相関が見られる（＝昨日の値と今日の値に相関があるが一昨日の値と今日の値には直接的な相関があるとは言えない）。の解釈でいいのかな。
Juliaという名の欲張り言語
最近Juliaというプログラミング言語を知りめっさ同感した。僕らは確かに統計処理、機械学習、データマイニング等がしたいのだけれどそれと同時にスピードも欲しいしmatlabみたく手軽に書きたい。ほんと仰る通りでどの言語も一長一短なんだよ。
なぜ僕らはJuliaを作ったか（原文）
2012年2月14日(火) | Viral Shah, Jeff Bezanson, Stefan Karpinski, Alan Edelman端的に言えば、僕らは欲張りだからだ。
僕らはMatlabのパワーユーザーだ。LispハッカーやPython使いやRuby使いもいる。髭が生える前からMathematicaを使っていたのもいるし、未だに髭が生えてない仲間もいる。常識的な人にはオススメしないくらい多くのグラフをR言語で描いてきた。そしてC言語は僕らのユートピアだ。
いま挙げた言語は大好きだ。どれも素晴らしいしパワフルだけど、科学計算、機械学習、データ・マイニング、大規模な線形代数演算、分散・並行コンピューティング、といった僕らがやるようなものにはどれも一長一短で、仕事に完璧にはまる機能もあれば何とも使い物にならないものもある。どれもトレード・オフなんだ。
僕らは欲張りだ。これじゃ十分じゃない。
僕らが欲しい言語はこんな感じだ。まず、ゆるいライセンスのオープンソースで、Cの速度とRubyの動的さが欲しい。Lispのような真のマクロが使える同図象性のある言語で、Matlabのように分かりやすい数学の記述をしたい。Pythonのように汎用的に使いたいし、Rの統計処理、Perlの文字列処理、Matlabの線形代数計算も要る。シェルのように簡単にいくつかのパーツをつなぎ合わせたい。チョー簡単に習えて、超上級ハッカーも満足する言語。インタラクティブに使えて、かつコンパイルできる言語が欲しい。
（そういえば、C言語の実行速度が必要だってのは言ったっけ？）
こんなにもワガママを言った上だけど、Hadoopみたいな大規模分散コンピューティングもやりたい。もちろん、JavaとXMLで何キロバイトも常套句を書きたくないし、数千台のマシンに分散した何ギガバイトものログファイルを読んでデバッグするなんて論外だ。幾層にも重なった複雑さを押しつけられるようなことなく、純粋なパワーが欲しい。単純なスカラーのループを書いたら、一台のCPUのレジスターだけをブン回す機械語のコードが生成されて欲しい。A*Bと書くだけで千の計算をそれぞれ千のマシンに分散して実行して、巨大な行列の積をポンと計算してもらいたい。
型だって必要ないなら指定したくない。もしポリモーフィックな関数が必要な時には、ジェネリックプログラミングを使ってアルゴリズムを一度だけ書いて、あとは全ての型に使いたい。引数の型とかから自動的にメソッドを選択してくれる多重ディスパッチがあって、共通の機能がまったく違った型にも提供できるようにして欲しい。これだけのパワーがありながらも、言語としてシンプルでクリーンなものがいい。
これって、多くを望みすぎてるとは思わないよね？
僕らがごまかしようのないほど欲張りなのは分かってるけど、それでもぜんぶ欲しいんだ。二年半ほど前、この欲にまみれた言語を作り始めた。まだ完成してないけど、そろそろ1.0のリリースの時期だ。僕らが作った言語の名前はJulia。すでに僕らの無礼な要求に9割方は応えてくれてるけど、ちゃんとした形になるためには僕ら以外の要求も聞かないといけない。だから、君がもし欲張りで理不尽でわがままなプログラマなら、ちょいとこいつを試してもらいたいんだ。
これであとはJVM言語だったら完璧だな。
時系列解析入門
大雑把に言うと、時系列解析とは過去のデータから未来を予測する分野。時系列データは通常のデータと異なり各値が強い相関を持つケース（非定常性）や”誤差が累積してゆく（ランダムウォークなど）”ケースが多く存在しており回帰分析など一般的に良く知られている手法とは相性が悪いことが多い。故に時系列解析特有のアプローチが必要。Rを使った時系列解析の大まかな流れとして
- Identification : Data preparation, then Model Selection
- Estimation and testing
- Application
具体的な手順の１つのアプローチとして以下のような流れがある。
- データをRで読み込む
- ts( )メソッドでdataを時系列データ化する
- 時系列データをplot( )して可視化
- acf( ) メソッドによって定常性を調べる
- もし非定常なら回帰分析を利用したDetrendingする
- それでも非定常ならdiff( )メソッドでDifferencingを行う
- acfと残差を確認する
- モデル選択；定常＝＞AR,MAモデル、非定常=>ARIMA,SARIMAモデル
- 仮説、モデルの妥当性を検討（もし芳しくないようならモデル選択まで戻る）
- モデルを基に予測
定常性(Stationary)とはなにか。その前に時系列データは確率過程として扱うということを理解する必要がある。時系列データにおいては時間軸のそれぞれの値を確率変数、データ自身をその実現値として捉える。サイコロを例に挙げると考えやすい。サイコロを振って出た目の分だけお金をもらえるとする。そのとき各出目が確率変数、もらえるお金が実現値。サイコロの例では確率変数が６つと限られているが、それを拡張しその確率変数の列が果てしなくつづいていると考える。そのときその確率変数の列を確率過程という。定常性に話を戻すと時系列解析に置ける定常は定常な確率過程のことを指し示している。
一言でいうと
時間が経っても全体で見ればその時系列は変わらないという性質のこと。
厳密に言うと、
期待値と自己共分散が時間tに依存せず自己相関が２点の距離(Lag) hにのみ依存する場合、それは定常な確率過程である（弱定常。強定常は分布が完全に一致し計量経済学においては稀であるらしい）。
数学的に言うと、
例えばある八百屋さんがお客さんの数を毎日数えているとする。今日までの来客数の平均と５日前までの来客数の平均が定常過程では同じになる。もしそうでなければお客さんの平均は時間tに依存していることになりこの場合非定常になるわけだ。自己共分散も同じように考えて、一貫しているなら定常、そうでないなら非定常。
acfはauto correlation functionの略で、自己相関を表している。この自己相関を
簡単に言うと、
一昨日のデータと今日のデータがどれほど相関があるか（もちろん一昨日と今日だけではない）
別の言い方をすれば
異なる時点間でその過程が何かしらの関係性を持っているかどうか
数学的に表現すると
この自己相関はもちろんどれほど関数を横にずらすかによっって変わってくる。上のacfメソッドはずらす距離(Lag)ごとの相関度を計算してくれる。Lag==0は関数をずらしていないことを意味するので必ず相関度１になる。もしLagが大きくなるにつれて徐々に相関度が下がって行き最終的に誤差の範囲（Rでは点線で示される）に収まるようならそれは非定常。翻ってLag==1から急激に相関度が下がり範囲内に収まるならそれは定常的なデータ。これでデータの定常性を調べることができる。自己相関は見ての通り自己共分散を基に計算される。しかし自己共分散はその単位の大きさが値の大きさに影響してしまうのに対して自己相関は基準化されたもの。Detrendingはまた次回〜
これが非定常な自己相関。次のがより定常的な自己相関。
関数とメソッド in Scala
関数とメソッド。一見同じように見えるが、メソッドはオブジェクトやクラスに結びついているもの。関数はそれらに結びついていない単独な第一級オブジェクト（少なくともScala利用者の立場的には）であり値として扱うことが出来る。僕みたいなめんどくさがりやには結構便利な機能で例えばあるメソッドの処理時間を調べたい場合。
def examine(post:String):Unit={
//何も返さない適当な処理
}
def time(f:String=>Unit,post:String)={
val start=System.currentTimeMillis
f(post)
val end=System.currentTimeMillis
println("it takes "+(end-start))
}
今回上のexamineの処理時間を調べたいとする。timeメソッドはStringを受け取りUnitを返す（つまり何も返さない）関数とStringを引数にとるメソッド。そしてこれら２つのメソッドはオブジェクト、NF2に属しているとする。２つの定義は上でしたので後は実行するだけ。
val t=NF2.time _
valf=NF2.examine _
//オブジェクトのメソッドは関数ではないので関数化するには メソッド名+ _
t(f,"失敗成功")
t(f,"彼は私に気があるのでしょうか？")
楽で良いー。
