Statistical Programming

有名な確率分布

確率分布とは確率変数と確率変数に対応したその変数の起こる確率の分布。具体例で考えるとわかりやすい。サイコロの目が確率変数で、それぞれの目の出る確率(もちろん1/6)。それらをx軸が確率変数、y軸がその変数の確率といった形で分布化したものが確率分布。確率密度関数といっても差し支えないと思う。この場合、x軸は1〜6の離散値。y軸はその変数の確率なので全て1/6。よって横一直線の関数ができる。
<一様分布>上のサイコロの例のように常に確率密度関数が常に一定の値をとる分布を一様分布U(a,b)という。この分布では以下のことが言える。

f(x) = k [a<=x<=b]    ||   0   [x<a || x>b]  ,  µ = (a+b) / 2 , σ^2 = ( (b-a)^2 ) /12    

※kは一定であり確率の総和が1であることから求める。

<二項分布>
一回の試行である事象Aの起こる確率がpである時にこの試行をn回繰り返す。事象Aがk回起こる確率分布B(n,p)は
    nCk * p^k * ( (1-p)^(n-k) )
※nCkは組み合わせを表現。変数ではない
また平均、分散は
    µ = np
    σ^2 = np(1-p)

正規分布
みんな大好き正規分布は頻度派統計ほどじゃないがベイズでも大きな役割を果たす。正規分布確率密度関数N(µ,σ^2)は
f(x) = 1/(√(2π)-1) * e^(- (x-µ)^2 / 2σ-2 )
eはネイピア数でe=2.71828....
<ベータ分布>
事前分布、事後分布として利用されることの多い分布だがベイズ統計に置いてはその分布の意味自体ではなく分布の形だけが重要なことも多いらしい。ベータ分布の確率密度関数Be(p,q)は
f (x) = k * x^(p-1) * (1-x)^(q-1)
※kは定数、0<x<1,0<p,0<q
µ = p / (p+q)
σ^2 = pq / ( (p+q)^2 * (p+q+1) )

kは確率の総和が1になることから求められる(規格化の条件)。

ベイズ統計入門

普通の統計では統計モデルを仮定しそれによってデータをを説明する。そのモデルを規定する重要な要素がそのデータの母数(平均µ,分散σなど)。ベイズ統計ではその母数とベイズの定理を組み合わせる。

ベイズの定理

P(H|D)=( P(D|H)*P(H) ) / P(D)

このHを母数θと捉えていく。つまりP(θ|D)はDの時にデータの母数がθである確率となる。さらに、連続値を扱うためにベイズの定理を確率ではなく分布、確率密度関数として解釈する必要が出てくる。

解釈の変更
事前確率P(θ)=>事前分布π(θ)
尤度P(D|θ)=>尤度関数f(D|θ)
事後確率P(θ|D)=>事後分布π(θ|D)

また、分母P(D)はデータを得る確率でありそれはデータを得られた後、定数になるのでベイズの定理を簡単な式に少し変更することが出来る。事後分布は分布なんだから右辺とイコールじゃなくて比例でも分布の関係性に大きく影響はでないし求めたいのはその分布のy軸の値(たぶん確率とはまた違うが確率等をもとに出来上がったもの)自体ではなくて(それも大事だが)どのθが最も可能性が高いのかについて。だから分母はある種無視していい。と思う。たぶん。

ベイズ統計公式
事後分布π(θ|D) ∝ ( 尤度f(D|θ) * π(θ) )

このθはベイズの定理のHと違って複数の離散値も連続値も表すことが可能。∝は比例を意味しとる。
古典的な頻度派との違いは主に母数の扱いとその焦点。

頻度派<=>ベイズ
母数:定数<=>確率変数
焦点;定数の値<=>確率変数の分布

要は、頻度派では特定の母集団がいてその母数は定数(ただしその母数の本当の値は神のみぞ知る)。僕らの手元にはそのサンプルデータがあり行うことはそのサンプルデータがどれくらいのその母集団と近いか、どれくらいその母集団にとって一般的なサンプルなのかを見極めることによって母集団を推定するとかとか。翻ってベイズ統計では手元にデータがありそれをもとに母集団の母数の分布を求める。母集団はいろいろな可能性が考えられるけど現在のデータからするとこういう母集団の分布がありうるよねーみたいな。何が大きく違うっていうと特定の母集団があると決めつけるか、母集団もいろいろな可能性があるよねーとゆるりと考えるか。頻度派ではそのいろいろな可能性の部分をサンプルデータに任せてて、ベイズでは母集団がその任を負う。

次に具体例を。コインの表の出る確率をθとする。コインを3回投げて3回とも表だった。θの事前分布がベータ分布Be(2,2)の時にθの事後分布は?そしてその事後分布から推定するθの値は?

これを考える前に前提知識としてベイズ統計のアプローチについてもう少し知る必要がある。アプローチは大別して2つある。

1. 自然な共役分布;多少無理があっても事前分布と事後分布の分布のタイプをそろえることで、公式を用いて簡単に計算出来るモデルを作る
2. MCMC法;多少計算に正確さが欠けても複雑なモデルをそのまま扱う

今回は共役分布を利用する。となるとベータ分布の公式を把握する必要がある。ベータ分布Be(p,q)は当たり前ながら確率密度関数

ベータ分布の公式
f(x)=k*x^(p-1)*(1-x)^(q-1)kは定数,p>0&&q>0
k= 1 / Be(p,q)
µ(平均値)=p/(p+q)
σ^2(分散)=pq/( (p+q+1)*(p+q)^2 )
M(最頻値)=(p-1)/(p+q+2)

この公式を用いて解いて行く。まず尤度f(D|θ)は表の出る確率がθであることから

尤度f(D|θ)=θ^3
事前分布π(θ)はBe(2,2)が与えられているので
π(θ)=k^(2-1)*(1-θ)^(2-1)
    =kθ(1-θ)
    ∝θ(1-θ)
事後分布∝尤度*事前分布なので
    =(θ^3)(1-θ)
    =θ^(5-1)*(1-θ)^(2-1)
    =Be(5,2)

よって事後分布はベータ分布Be(5,2)であることがわかった。
上記のことから2つの事象しか現れないベルヌーイ試行において以下のことが言える。

ベルヌーイ事象の事後分布(事象A,Bしか存在しない)
事象をn回行い事象Aがr回でたとする。事象Aが起こる確率をθ、事前分布をBe(p,q)とすると

事後分布=Be(p+r,q+n-r)

最後にベータ分布の平均、分散、最頻値を求める。

平均=5/(5+2)=5/7
    =0.74
分散= 5*2 / ( (5+2)^2 * (5+2+1) )
    = 5/196
    =0.026
最頻値=(5-1)/(5+2-2)
   =0.8

よって平均を用いた場合0.74,最頻値の場合0.8の確率で表が出るコインだと考えられる。

意中の相手への告白がうまくいくかどうか過去の事例を解析し成否を占う。with Scala

告白したい!でもうまくいくか不安…だから掲示板にでも相談しよう…こんなケースってよくありがちだけどそれらのデータを蓄積して解析してるシステムってあんまないなーって思う。だからちょっと作ってみよう!と思い立ってのこと。とは言っても迷惑メールフィルターとほぼ一緒の作りだし、使っているのはベイズの基本的な知識とScalaだけ。あとは形態素解析(文章を意味のある最小限の単語に分割する)のライブラリとしてkuromojiにはお世話になってる。具体的にどうやって解析して行くかというと、相談内容の単語の出現頻度をベイズ的に処理して過去のデータからいかにその投稿内容やコメントが成功もしくは失敗に近いかを計る。成功か失敗、より近い方を予測として投稿者に教えてあげる。

1. 投稿者が相談内容を投稿;「告白うまくいくでしょうか?」

2. 相談内容の文字列を意味の通る最小限の単語のリストに分割。

3. 過去のデータからそれぞれの単語の出現頻度を求める

4. ベイズの定理を用いて、その投稿内容がどれほど成功もしくは失敗に近いか求める

5. より確率の高い方を予測として投稿者に知らせる;"成功します"、"失敗します"

ではそもそもベイズの定理とは何か。

P(H|D)=( P(D|H)P(H) ) / P(D) 

これだけ。上の関係式は

Dの時にHが起きる確率=

( Hの時にDが起きる確率*Hが起きる確率 ) / Dが起きる確率 

を意味しとる。このP(D|H)を尤度、P(H)を事前確率、P(H|D)を事後確率という。なぜこの関係式が得られるかは基本的な確率公式を組み合わせたらすぐわかるからまあ割愛するとして、次にベイズの更新について。ベイズの更新とは同じ事象を繰り返し行う場合に前回の事後確率は今回の事前確率として扱うことができるということ。

例えば、2つの壷A,BがあるとしてAには白い玉が2個と赤い玉が1個、Bには白い玉が1個と赤い玉が2個入っている。A,Bいずれかから玉を2つ取り出してみた。この時点ではA,Bどちらから取り出したかわからないが取り出した玉は赤い玉1個と白い玉1個だった。さて、自分はどちらの壷から玉を取り出していたのかそれぞれの確率を求めたいとする。
※P(D)はすべてのパターン(ここでは壷がAであるパターンとBであるパターン)の尤度P(D|H)の和と同値である。つまり今回は
P(D)=(P(D|A)*P(A)+(P(D|B)*P(B)

壷から玉を一つ取り出してそれが赤い場合

P(A|赤)=

( (P(赤|A)*P(A) ) / ( (P(赤|A)*P(A)+(P(赤|B)*P(B) )

=( (1/3)*(1/2) ) / ( (1/3)*(1/2)+(2/3)*(1/2) ) 

=0.3333333

P(B|赤)=

( (P(赤|B)*P(B) ) / ( (P(赤|A)*P(A)+(P(赤|B)*P(B) )

=( (2/3)*(1/2) ) / ( (1/3)*(1/2)+(2/3)*(1/2) ) 

=0.6666667

つまり、赤い玉一つが出て来た段階で壷Aから取り出した確率は3分の1、Bは3分の2ということがわかった。同じことを2つ目の玉にも行う。ただし先ほどは事前確率P(A),P(B)両方とも2分の1としたが現時点でP(A)=0.3333333,P(B)=0.6666667ということが判明しているのでそれらの値を事前確率として活用する。これがベイズの更新。

壷から玉を一つ取り出してそれが赤い場合

P(A|白)=

( (P(白|A)*P(A) ) / ( (P(白|A)*P(A)+(P(白|B)*P(B) )

=( (2/3)*0.3333333 ) / ( (2/3)*0.3333333+(1/3)*0.6666667 ) 

=0.5

P(B|白)=

( (P(赤|B)*P(B) ) / ( (P(赤|A)*P(A)+(P(赤|B)*P(B) )

=( (1/3)*0.6666667 ) / ( (2/3)*0.3333333+(1/3)*0.6666667 ) 

=0.5

結果、半々の確率でした〜ということがわかる。ここまででだいたいベイズの基本はおっけ。ベイズの定理の良いところは人間の直感に近いプロセスを踏んでいるところだと思う。

基本的にはこれと同じこと単語レベルで行う。先ほどはDは玉の色、赤か白だったけどこれを今回は単語として捉える。

P(成功|単語が得られたよ)= ( P(単語が得られたよ|成功)P(成功) ) / P(単語が得られたよ)

P(失敗|単語が得られたよ)= ( P(単語が得られたよ|失敗)P(成功) ) / P(単語が得られたよ)

先ほどのベイズのやり方と今回のやり方の唯一違うところは分母を計算する必要がないということ。先程示した通り異なるパターンの間では分子が異なり分母は等しい。故にどちらの確率の方が高いかを比べる上で分母の計算自体は必要なくなる。分子だけを見ればいい。

P(A|D):P(B|D)=P(D|A)*P(A):P(D|B)*P(B)

2つを比べて確率の高い方が真ということになる。この2つの確率を比べるわけだけど全体の単語数同じになるわけでそうするとP(単語が得られたよ)が成功時と失敗時、同じ結果になる。だから今回はそこは省いてしまう。というお話、

P(成功|単語が得られたよ)=  P(単語が得られたよ|成功)P(成功)

 

P(失敗|単語が得られたよ)=  P(単語が得られたよ|失敗)P(成功) 

まずはそのために形態素解析のライブラリであるkuromojiを活用して文章を単語に分割する必要がある

kuro.scala
def mkTokenS(content:String):Stream[Token]
=Tokenizer.builder().build().tokenize(content).toStream.filter( (x:Token)=>x.getBaseForm() != null)
    //StringからToken作成

以前形態素解析の記事を書いたので詳しいことは割愛。
先ほどのベイズのやり方と今回のやり方の唯一違うところは分母を計算する必要がないということ。先程示した通り異なるパターンの間でくなる。分子だけを見ればいい。
P(A|D):P(B|D)=P(D|A)*P(A):P(D|B)*P(B)
2つを比べて確率の高い方が真ということになる。

NaiveFilter.scala

package textmining
import org.atilika.kuromoji._
import scala.collection.mutable.ListBuffer

object Tes extends App{
  println(NaiveFilter.data)
  println(NaiveFilter.Sdata)
  println(NaiveFilter.Fdata)
  println(NaiveFilter.stringSdata)
  println(NaiveFilter.stringFdata)
  NaiveFilter.tokenSdata.map(x=>x.getBaseForm)
  NaiveFilter.tokenFdata.map(x=>x.getBaseForm)
  println(NaiveFilter.occurS)
  println(NaiveFilter.occurF)
  println(NaiveFilter.dicF.map(x=>x.getBaseForm))
  println(NaiveFilter.proS)
  println(NaiveFilter.proF)
  println(NaiveFilter.Saft)
  println(NaiveFilter.Faft)
  NaiveFilter.compare("いっくおー")
}
object NaiveFilter{ 

  def addtest(post:String)={

  }

  var data:ListBuffer[(Boolean,String)]
          =ListBuffer( (true,"成功すると思うよ"),(true,"絶対成功するって!!"),(false,"きついと思うなあ失敗するんじゃない?")
              ,(true,"君可愛いんだし言っちゃえよ"),(false,"失敗するに一票"))
          //事前の値は学習データ。dataのみmutable

  var Sdata=data.filter(x=>x._1==true)  //成功data
  var Fdata=data.filter(x=>x._1==false) //失敗data

  val stringSdata=Sdata.map(x=>x._2).mkString   //成功例の文字列データ
  val stringFdata=Fdata.map(x=>x._2).mkString   //失敗例の文字列データ

  val tokenSdata=kuro.mkTokenS(stringSdata)     //成功例のトークンリスト
  val tokenFdata=kuro.mkTokenS(stringFdata)     //失敗例のトークンリスト

  val numTokenS:Double=tokenSdata.length
  val numTokenF:Double=tokenFdata.length

  val dicS=tokenSdata.toSet //成功例の重複しないトークンセット
  val dicF=tokenFdata.toSet //失敗例の重複しないトークンセット

  val occurS=dicS.toList.map(x=>tokenSdata.count(y=>y.getBaseForm==x.getBaseForm))  //dicSの全ての要素の出現回数を算出
  val occurF=dicF.toList.map(x=>tokenFdata.count(y=>y.getBaseForm==x.getBaseForm))  //dicF 

  val proS=occurS.map(x=>x.toDouble/numTokenS) //成功の時の各トークン(重複しない)の出現回数/成功時の全トークン数==>成功時各尤度
  val proF=occurF.map(x=>x/numTokenF) //失敗の時の各トークン(重複しない)の出現回数/失敗時の全トークン数==>尤度

  val Saft=proS.reduce( (a,b)=>a*b)  //成功時の事後確率=成功時の全てのトークン(重複しない)の尤度を掛け合わせたもの
  val Faft=proF.reduce( (a,b)=>a*b)  //失敗時

  def compare(post:String)={
    val preS=Saft   //学習データに置ける事後確率が事前確率へと変わる。これは成功時の事前確率
    val preF=Faft   //こちらは失敗時
    (true,post)+=:Sdata
    (false,post)+=:Fdata
    println(preS)
    println(Saft)
  }
}

さて、ここまででいくつか問題点がある。まず上のコードではvalを多用しているがvalは再利用性がない。

    var a=3                   //> a  : Int = 3
    def b=a                   //> b: => Int
    val c=a                   //> c  : Int = 3
    a=5
    a                         //> res0: Int = 5
    b                         //> res1: Int = 5
    c                         //> res2: Int = 3

このように変数の中の変数の値を更新したとしてもvalの場合その変数に反映されない。ここの例ではbもcもaの値を返すが、それらを定義した後にaの値を変更したとするとbの値はaと連動して変更されているがcの値は最初のaの値がそのままになっている。つまりvalを使うと似たようなコードを書かないと行けなくなる。というわけでvalを使っているところをdefに書き換えて変数じゃなくて関数にしてしまう。それに加えて、成功事例と失敗事例でほぼ同じコードを書いているのだからこれも部品として再利用出来る形にしたい。あと、事前確率の計算が入ってないの。というかこれじゃあ各尤度を求めるだけになっとる。んで下のが改善版。

NaiveFilter.scala
package textmining
import org.atilika.kuromoji._
import scala.collection.mutable.ListBuffer

object Tes extends App{
  NaiveFilter.examine("可愛い失敗")
}
object NaiveFilter{ 

  var data:ListBuffer[(Boolean,String)]
=ListBuffer( (true,"可愛成功"),(true,"可愛成功"),(false,"失敗"),(false,"失敗"))
//元の値達は学習データ。data,Sdata,Fdataはmutable。falseの学習データには失敗しか入ってないけどtrueの方に可愛と成功が入ってる。これによって失敗という単語の重みは成功や可愛という単語よりも大きくなる。故に可愛いと失敗が1語づつのStringを調査するときは失敗に傾く。

  var Sdata=data.filter(x=>x._1==true)  
//成功data。mutable
  var Fdata=data.filter(x=>x._1==false) 
//失敗data。mutable

  def preS=(Sdata.length-1).toDouble/data.length    //examineでSdataには最初に調べる項目が追加されるがdata自体には関数の最後に追加される。それによって成功事例が多く見積もられるのを防ぐためにSdataに-1をしてる
  def preF=1-preS
  //成功、失敗それぞれの事前確率。examineの中で使う。ただし追加したpostが影響しないようにvalを使おうと思ったけどdataとSdataとの関係が崩れてしまうからやめた

  def getString(list:ListBuffer[(Boolean,String)])
=list.map(x=>x._2).mkString
  def stringSdata=getString(Sdata)  //成功例の文字列データ
  def stringFdata=getString(Fdata)  //失敗例の文字列データ

  def getToken(stdata:String)=kuro.mkTokenS(stdata)
  def tokenSdata=getToken(stringSdata)      
//成功例のトークンリスト
  def tokenFdata=getToken(stringFdata)      
//失敗例のトークンリスト

  def getnum(tokens:Stream[Token])=tokens.length
  def numTokenS:Double=getnum(tokenSdata)
  def numTokenF:Double=getnum(tokenFdata)


  def norep(tokens:Stream[Token])=tokens.toSet
  def dicS=norep(tokenSdata)    
//成功例の重複しないトークンセット
  def dicF=norep(tokenFdata) 
//失敗例の重複しないトークンセット

  def occur(dataset:Set[Token],tokens:Stream[Token])
:List[Int]=dataset.toList.map(x=>tokens.count(y=>y.getBaseForm==x.getBaseForm)) 
  def occurS=occur(dicS,tokenSdata)  
//dicSの全ての要素の出現回数を算出
  def occurF=occur(dicF,tokenFdata)  
//dicF 

  def prob(oclist:List[Int],numToken:Double)
:List[Double]=oclist.map(x=>x.toDouble/numToken)
  def proS=prob(occurS,numTokenS) 
//成功の時の各トークン(重複しない)の出現回数/成功時の全トークン数==>成功時各尤度
  def proF=prob(occurF,numTokenF) 
//失敗の時の各トークン(重複しない)の出現回数/失敗時の全トークン数==>尤度
  //ここまでが尤度の算出

  def zipDicPro(dic:Set[Token],prolist:List[Double])
=dic.map(x=>x.getBaseForm).zip(prolist).toMap
  def zipS=zipDicPro(dicS,proS)
  def zipF=zipDicPro(dicF,proF)

  //単語リストと各単語の尤度をzipしてMap化。

  def aft(prolist:List[Double]):Double
=prolist.reduce( (a,b)=>a*b)

  def examine(post:String)={
    val Snow=(true,post)
    Snow+=:Sdata
    val Fnow=(false,post)
    Fnow+=:Fdata
    //両方にadd
    println("成功事例 : "+Sdata)
    println("失敗事例 : "+Fdata)

    val dic=getToken(post).toSet.toList  
//postの重複のなトークン集合
    println("重複しない単語リスト : "+dic.map(x=>x.getBaseForm))
    val allyudoS=dic.map(x=>zipS(x.getBaseForm))
    println("成功時各尤度 : "+allyudoS)
    val allyudoF=dic.map(x=>zipF(x.getBaseForm))
    println("失敗時各尤度 : "+allyudoF)
    println("成功事前確率 : "+preS)
    println("失敗事前確率 : "+preF)
    val Saft=aft(allyudoS)*preS
    val Faft=aft(allyudoF)*preF

    println("成功事後確率 : "+Saft+" <=> 失敗事後確率 : "+Faft)
    if(Saft>Faft){
            println("成功するでしょう") 
            Fdata-=Fnow
            Snow+=:data
    }
    else{
            println("失敗するでしょう")
            Sdata-=Snow
            Fnow+=:data
    } 

  }

学習データには失敗しか入ってないけどtrueの方に可愛と成功が入ってる。これによって失敗という単語の重みは成功や可愛という単語よりも大きくなる。故に可愛いと失敗が1語づつのStringを調査するときは失敗に傾く。また、examineでSdataには最初に調べる項目が追加されるがdata自体には関数の最後に追加される。それによって成功事例が多く見積もられるのを防ぐためにSdata.lengthに-1をしてる。わかりやすいように短くてわかりやすい単語を学習データにもexamineのパラメータにも使ったけど別にふつーーの口語体の文章を入れたも全然おっけ。今回詰まったとこは尤度を求めるところ。examineのパラメータにあって学習データにない単語が出てきた時にその尤度は当然ながら成功事例では全て0、失敗事例では全て1になってしまった。しかもそれらを最終的に掛け合わせるのだから結果がすごいことに…それを防ぐために若干の誤差が生じるけどexamineのパラメータの値を成功、失敗両方のデータに加えて必ずすべての単語の出現頻度を1以上にしている。この実装は効率が著しく悪いのと、examineのパラメータが長くなると巨大なリストを作成しなくてはならず計算が一定を超えるとできなくなるのでファイルへの読み書きやリストを一定の長さで分割したりくっつけたりする必要がある。それと最も大きな問題は文字数が多くなるにつれて確率がやたらと小さな数字になってしまい桁数溢れを起こしてしまうということ。この問題に対してはすべての尤度と事前確率に対して対数Logを取ることで解決している。

a*b*c*.......*n

Log(a*b)=Log(a)+Log(b)

Log(a*b*c......*n)=Log(a)+Log(b)+.......+Log(n) 

 本来の尤度の積は一番上のやつ。だがこれは積なのでどんどん桁が小さくなってしまう。故に積の計算をLogの足し算へと変換する。大小関係は変わらないので問題なし。対数と取った場合の最終的なコードはGitHubにおいてる。

StatsBox/src/main/scala/classifier at master · kyu999/StatsBox · GitHub

言語と彼女

何度見ても言い得て妙だと思う。プログラミンング言語を彼女に例えた文言。

PHP は、あなたが高校時代のある夏、不器用ながらも付き合った初めての彼女です。もっと真剣な関係を築こうとしてはいけません。この子は複雑な問題を抱えています。

PerlPHPのお姉さんです。あなたからしたらちょっと年上かもしれませんが、90年代は人気者として鳴らしていました。Larry Wall と長い間付き合っているので、今は理想も低くなったし、かなりのブッサイクになりました。彼は「お前らが何を言おうと知らん、俺は彼女を愛している!」と言います。そんな人は彼だけです。

Rubyスクリプティングファミリーの中でもイケてる子です。初めて彼女を見たときは、その美しさに息を呑んだことでしょう。それに、一緒にいて楽しかった。当時はちょっとのろまで天然という印象がありましたが、ここ数年でだいぶ大人びてきました。

PythonRuby良識あるお姉さんです。彼女はエレガントで、洗練されています。もしかしたら完璧すぎるほどに。ほとんどの男性は、「Pythonが好きじゃないヤツとかいるの?」と思うでしょう。Pythonが好き、そうでしょうねぇ。エッジーでロマンティックなRubyの「つまんない版」として見ているだけですよ。

Java は成功したキャリアウーマンです。彼女と働いたことがある一部の人は、彼女の地位は実力ではなく、中間管理職に好印象を持たせる要領の良さによるものだと感じています。あなたは聡明な彼女と一緒に身を固めたいと感じるかもしれません。ただ、「そこはそうじゃないわよ、ああもう、あなたっていつも間違ったType Interfaceを使うしセミコロンも抜けてるじゃないの」と口うるさく言われる日々に備えてください。

C++Javaの従姉妹です。あらゆる側面でJavaと似ていますが、一番の違いは、もっと純粋な時代に育っていて、Protectionを使いたがらないところです。ここでいうProtectionは、Automatic Memory Managementのことですよ。何だと思ったんですか?

C はC++のお母さんです。白ヒゲの古いハッカーの前で彼女の名前を言ってみたら、間違いなく目をキラキラさせながら回想するでしょう。

Objective C もC家の一員です。あの変な宗派にだいぶ前に入信していて、外の世界の人と付き合おうとしません。

Haskell, Clojure, Scheme や彼女たちの仲間は、大学生の頃のとある夏にこの上なく楽しい時間を一緒に過ごしたであろう、流行に敏感なアート系のかしこい女の子たちです。本気で試された、と初めて感じさせてくれた女性たちでした。もちろん、彼女たちと真剣な関係になることはありえません(自分が一番よく分かっていますよね)。なのに、今でも考えてしまいます。「もしあるとしたら?」

C# については、彼女の家族の評判を理由に避けるかもしれません。でも彼らによると、ここ数年の間にだいぶまともになったそうです。一度関わったら、もう家族の一員なんだ。分かるかい?データベースがいるって?彼女の兄弟のMSSQLが何とかしてくれるよ。住む場所がない?ハッ、彼女の父さんがAzure街にマンションを買ってくれるってさ。なになに、フレンドリーすぎる親戚たちのことで思い悩んでいるんだって?いいや、逃さないよ。君はもう家族の一員なんだ。聞いてるか?

Javascript - ねえ、彼女はたしかPHPに会うよりもずっと前にファーストキスした子だよね?今はどうしているのかなぁ。ここ数年の間にずいぶんとキャリアが軌道に乗ったみたいだ。昔のよしみで、いつか近況報告をしあえたらいいんだけど…(jQuery で全身ドレスアップした彼女を見かける)うわぁ…すごくあかぬけて綺麗になったなぁ… 

 jQueryでドレスップってのがまた面白い表現。

2群の検定 with R

今まで「Scala Days」というタイトルでScalaについて主に書いてたんだけど最近python使わなくちゃいけなくなったし統計に関することもメモりたいからタイトル名を「Statistical Programming」へと変更。というわけでとりあえず統計の復習を少しばかり。

2つのグループに関するデータがあり、グループ間に有為な差が見受けられるか否かを調べたい。そういったケースのプロセスについて。まずその2つのグループが独立した異なるグループ(独立2群)かそれとも単一のグループだが異なる時期に実験を行ったもの(関連2群)なのか判断する必要がある。

独立2群の例;A組の生徒は3ヶ月間補習を受ける。B組の生徒は補習を受けない。3ヶ月後に両方の生徒は試験を受けその試験結果に有為な差がないかどうか調べる

関連2群の例;A組の生徒はまず試験を受ける。その後3ヶ月間特別な補習を受ける。補習が終了したらA組の生徒は再び試験を受け、最初の試験結果と補習後の試験結果に有意な差が見受けられるかを調べる

関連2群の場合、

(補習後の試験結果)ー(補習前の試験結果)

をした後に1標本t検定を行う。

> data=read.csv("試験.csv",header=TRUE)
> data
    om  am cs
1   94  97 96
2   81  94 99
3   85  89 86
4   81  81 77
5   79  85 80
6   96  99 99
7   67  67 79
8   96 100 99
9   70  77 81
10  82  84 90
11  82  92 97
12  80  98 95
13  89  94 98
14  75  88 93
15  91  94 95
16  57  86 82
17  95  92 99
18  90 100 99
19  86  99 98
20  89  85 89
21  73  92 96
22  85  73 68
23  97  97 97
24  85  98 87
25  80  84 88
26  99  99 97
27  87  93 89
28  91  90 94
29  88  98 98
30 100 100 99
31  87  76 96
32  96  96 94
33  84  96 94
34  78 100 97
35  98  98 97
36  48  55 59
37  58  83 84
38  91 100 94
39  89  90 95
40  67  76 83
41  69  86 98

> origin=data$om
> result=data$cs
1標本t検定
> differ=result-origin
> t.test(differ)

    One Sample t-test

data:  differ
t = 5.4028, df = 40, p-value = 3.256e-06
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
  4.885244 10.724512
sample estimates:
mean of x 
 7.804878 

p-valueは0.05より明らかに小さいから補習の前後に有為な差があると考えられる。

翻って独立した2群である場合、まずコルモゴロフ・スミルノフ検定
を使って正規分布かどうかを調べてみる。帰無仮説は「標本分布が正規分布に従う」なのでp-valueが十分に大きかったら正規性が保証される。

コルモゴロフ・スミルノフ検定
>ks.test(origin,"pnorm",mean=mean(origin),sd=sd(origin))

    One-sample Kolmogorov-Smirnov test

data:  origin
D = 0.124, p-value = 0.5545
alternative hypothesis: two-sided

 警告メッセージ: 
In ks.test(origin, "pnorm", mean = mean(origin), sd = sd(origin)) :
   コルモゴロフ・スミノフ検定において、タイは現れるべきではありません 

ちなみに上記の例のようにデータに同値の値が存在している場合p-valueが正しく算出されないことがあるため別のパッケージをダウンロードしたほうがいい。
もしここで正規性に疑いがあると判断出来る場合、t検定ではなくWilcoxonの順位和検定(Mann-WhitneyのU検定)によって2群の相違を調べる。この手法はt検定とは異なりノンパラメトリック。帰無仮説は「両群が同じ母集団から抽出された」。

Wilcoxonの順位和検定
> wilcox.test(origin,result)

    Wilcoxon rank sum test with continuity correction

data:  origin and result
W = 479, p-value = 0.0008022
alternative hypothesis: true location shift is not equal to 0

 警告メッセージ: 
In wilcox.test.default(origin, result) :
   タイがあるため、正確な p 値を計算することができません 

p-valueは十分小さく仮説は棄却できる。よって2群の間に有為な差が見受けられる。

次に、2群の分散が極端に異なる場合を避けるためにF検定を行う(分散分析で登場するF検定とは計算方法が異なる)。帰無仮説は「二群の母分散は等しい」であり、p-valueが十分大きかったら等分散性が保証される。

等分散性を調べるF検定
> var.test(origin,result)

    F test to compare two variances

data:  origin and result
F = 1.7627, num df = 40, denom df = 40, p-value = 0.07678
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 0.9399903 3.3053492
sample estimates:
ratio of variances 
          1.762667 

0.07678と十分大きいため仮説は採択される。

※t検定において等分散性はさほど影響を及ぼさないことがわかっているらしいっす

最後に正規性が保証されている場合、Welchのt検定を行う。

> t.test(origin,result)

    Welch Two Sample t-test

data:  origin and result
t = -3.3155, df = 74.335, p-value = 0.001417
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -12.495062  -3.114694
sample estimates:
mean of x mean of y 
 83.29268  91.09756 
alternative hypothesisとは帰無仮説が棄却された時に得られる対立仮説のこと。p-valueが十分小さい時にこの対立仮説を採用する。ノンパラメトリックな検定は使い勝手がいいけれどその分適切に用いられたパラメトリックな検定より正確性は落ちるらしいからできることならパラメトリックな検定を使用出来る時は使用しよ

素数とビット計算 with python

素数を求める関数を作ること自体は簡単だけど速く素数を見つけ出す計算はなるほどなと思わされた。そんなの思いつかんって。こちらが普通の素数を求める関数。特に変わったとこはなし。

prime
def is_prime(n):
    for c in range(2,n):
        if n%c ==0:
            return False
    return True

次に早く計算が出来るように改良されたアルゴリズム。さっきは引数のnの一つ手前までループしてたけど今回は 1+√n の一つ手前まで。つまり√n まで。ここだけトリッキーなところだけど理論自体はそんなに難しくない。まず、nが14だとすると、1,2,7,14でもって14を割り切ることができる。この時、2が14を割り切るとわかってしまえば7も14を割り切るってわかる。つまり対になっているものが存在するのだからほぼ半分くらい計算を省略できるんでない?っつーことか。ここで疑問に思うのがなんでn/2じゃなくて√nなんだろうということ。この対になっているのを求めるにはこの場合14/2=7をすればいい。これらによって2で割り切れるとわかると同時に7でも割り切れることがわかる。つまり、 n%a==0 が成り立つ時 n%(n/a)==0 もまた成り立つ。
ここで横軸a,縦軸y、そしてf=a,f=n/aの2つの関数(数学的な)を思い描いてみると、f=aは右肩上がりの直線、f=n/aは右肩下がりの曲線になっている。先ほどのアルゴリズムはこのf=aの直線のみを活用してnまでループしていたけど今回はこの2つの関数の交点にあたるところまでループを行えばいい。たぶん。つまり、n=14の時f=aの関数で2を調べることと、f=n/aで7を調べることが同義であるということか。

a=n/a
a*a=n
a=√n

よって√nまでのループでおk。視覚に訴えないとこれきついな…

fast_prime
def fast_is_prime(n):

    upper = int(n**0.5)

    prime=True

    for c in range(2,upper+1):

        if n%c ==0:

            prime = False

            break

    return prime

first-n_prime()関数は単純に引数nまでの素数をリストアップするだけ。

first_n_prime
def first_n_prime(n):
    ans=[ ]
    candidate=2
    while len(ans) < n:
        #find a new prime
        if fast_is_prime(candidate):
            ans.append(candidate)
        candidate+=1
    return ans

次はビットの計算。Cでは単純にビット演算子を使ったけど今回はまた別な方法。digitを求めるのに%を使っているのが憎いね。自分で実装してみてちょっと困ったのがStringで使えるメソッドとListで使えるメソッドがあんまり共通してなかったこと。StringでfindメソッドなところをListでは object in list みたいな書き方をしてるし。findでは探し物が見つからなかった時に-1を返すのにlistのindexメソッドの方ではエラーになるし。Scalaみたくコレクションをまたいで同じメソッドで処理ってのがあんまりできないんだな。

bit
 def int2bin(x):
    L=[ ]
    while x>0:
        digit=str(x%2)
        L.append(digit)
        x/=2
    L.reverse()
    return "".join(L)
another_bit
def int2binNew(x):
    ans=""
    while(x>0):
        digit=str(x%2)
        ans=digit+ans
        x/=2
    return ans