Statistical Programming

ピアソンの相関係数を実装 in Scala

レコメンデーション等をするためには各アイテムの相関度などを調べる手法があり、そこで必要なのがユークリッド距離やピアソンの相関係数など。今回実装してみたのはピアソンの方。

f:id:kyu9q:20140218202227p:plain

 

2つの変数X,Y(X,Yに属する値のリストみたいなもん。実装ではStream[Double]という形を取っている)がありそれらの相関度をしらべたいときに上の公式を用いて相関係数を求めることが出来る。各偏差(Xi-平均のところ)を求めるところが共通しているからそこはメソッドとしてまとめれそう。それと分母分子両方に1/(n-1)が入っててこれは相関係数だけを求めるなら省くことが出来るけど標準偏差とかも正確に出したいから省略しないでおく。下のは省略した形。

f:id:kyu9q:20140218203345p:plain

 

とりあえず適当にデータフレームのクラスを作ってみる。

 

class Dataset(x:Stream[Double],y:Stream[Double]){} 

 後々のことを考えるとこのクラスはどうかと思うけどとりあえず2変量のケースにのみ考える。xは例えばあるクラスにおける生徒達の中間試験の結果、yは期末試験の結果。次に各変数の平均値を求める。

 def ave(data:Stream[Double]):Double

=data.reduce( (a,b)=>a+b)/data.length

 平均を使って各変数のそれぞれの値の偏差を求める。分母で各偏差の2乗も求めなければならないのでさらに各要素の偏差の2乗を求めるメソッドもつくる。

 

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))) 

 

各偏差の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))) 

//各偏差の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)

//2つのデータセットのピアソン相関関係

 

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 

 ちなみにもしX,Yの要素数が等しくない場合zipメソッドのおかげで要素数の短い方に合わせて長いほうの余剰をカットする。