2018年5月21日月曜日

Particle Smoothingの簡単な解説とJuliaで実装

Godsill et al. (2004)を基に、SMCを用いた状態空間モデルのスムージング(平滑化)を紹介する。モデルは隠れマルコフモデルを想定する。つまり状態変数x_tはマルコフ連鎖で、その密度はX_t|_x_t-1 = f(X_t|x_t-1)dx_tで、データy_tはx_tで条件付けると独立となり、その密度はY_t|X_t = g(Y_t|X_t)dy_tとする。なおこの論文の方法だと、モデルの密度関数が評価できないと使えない。

まず任意の時点tにおいて、Filter density p(x_t | y_1:t)を得ているとする。ここでy_1:t := (y_1, y_2,....y_t)。つまりまずSMCを用いてp(x_t | y_1:t)の近似、p_i(x_t | y_1:t) , i=1,2,...Nを得ているとする。簡単な計算よりp(x_1:T|y_1:T) = p(x_T|y_1:T)Π_t p(x_t|x_t+1:T, y_1:T) t = from 1 to T-1と分解できる。さらにモデルのマルコフ性より、
 p(x_t|x_t+1:T, y_1:T)∝p(x_t|y_1:t)f(x_t+1|x_t)となる。故に以下のアルゴリズムを得る。

まず前向きにフィルタリングを行った後、後ろ向きにスムージングを行うので、Forward Filtering-Backward Smoothingともよばれる。数学的に厳密な評価は、例えば Del Moral et al. (2009)やJasra (2013)等がある。

下にSV model、X_y=phi*X_t- + Normal(0,sig^2), Y_t=beta*exp(X_t/2)*Normal(0,1)に対する、上のParticle Smoothingの結果とJuliaのコードをのせる。パラメタはphi = 0.98, sig = sqrt(0.6), beta = 0.5とし、上のモデルにそって1000個観測データを発生させ、それを用いてN =10000個のパーティクルを用いてスムージングを行った。上にも書いた通り、まずSMCを用いてフィルタリングを行い(今回はブートストラップ・フィルタを用いた)、その後その結果を用いてスムージングをアルゴリズムに沿って行った。リサンプリングはシステマティック・リサンプリングを用いたので、合わせてそのコードものっけおておく。







 #### Code of SMC smoothing based on Godsill et al. 2004
 #### The Model is given by X_y=phi*X_t- + Normal(0,sig^2), Y_t=beta*exp(X_t/2)*Normal(0,1)

 ### input 1 data, 2 number of particles(N), 3 phi, 4 sig, 5 beta

 ### output 1 approximated filtering density, 2 approximated smoothing density

 ### I use systematic systematic resampling


 function SMC_Smoothing_SV(data, N::Integer, phi::Float64, sig::Float64, beta::Float64)

 ### initial setting

 @assert N > 0 "N must be a natural number"

 ###set sample sieze
 n = length(data)

 ###particles
 particles = zeros(n+1,N)

 ###weights
 weights = zeros(n+1,N)

 ###particles for smoothing
 sparticles = zeros(n+1,N)

 ###weights for smoothing
 sweights = zeros(n+1,N)

 ###empirical filtering density
 efilter = zeros(n)

 ###empirical smoothing density
 esmoothing = zeros(n)

 ###effective sample siez
 ESS = zeros(n)

 ###resampling criteria
 th = N / 2

 #########################Markov kernel######################################

 function Mkernel(x0, x1, phi, sig)

  return( -0.5*log(2.0*pi) - log(sig) -0.5*((x1 - phi*x0).^2) /(sig^2))

 end
 ###########################################################################


 ####################log Potential#######################################

  function Potential(x,y,beta)
      return(-0.5*log(2.0*pi) -0.5*log(exp(x)*beta^2) -0.5*y^2/(exp(x)*beta^2) )
  end
 #########################################################################

  #####initial step

  for j in 1:N

      particles[1,j] = rand(Normal(0,1))*sqrt(sig^2/(1-phi^2))
      weights[1,j] = 1/N

  end
  ###############

  #### start filterling ####

  for i in 1:n

      ### probagate particles and correct weights

      particles[i+1,:] = phi*particles[i,:] + rand(Normal(0,1),N)*sig

      for j in 1:N

          weights[i+1,j] = weights[i,j] + Potential(particles[i+1,j], data[i],beta)
      end

      ### calculate normilized weights
      for j in 1:N
          weights[i+1,j] = exp(weights[i+1,j])
      end

      sumwei = sum(weights[i+1,:])

      for j in 1:N
          weights[i+1,j] = weights[i+1,j] / sumwei
      end

      ### calculate ESS and check resampling
      sumsqwei = weights[i+1,:]' * weights[i+1,:]

      ESS[i] = 1.0 / sumsqwei

      if ESS[i] < th

          particles[i+1,:] = sysResampling(particles[i+1,:], weights[i+1,:], N)

          weights[i+1,:] = 1.0 / N

      else

          particles[i+1,:] = particles[i+1,:]

          weights[i+1,:] = weights[i+1,:]

      end

      println("Just finished filtering at time $i")
  end
   ###### end filtering #####

   ### calculate empirical filtering density
   temp = particles[2:n+1,:] .* weights[2:n+1,:]

   for i in 1:n
       efilter[i] = sum(temp[i,:])
   end

   #### start smoothing ####

   samples = copy(particles[n+1,:])

   sparticles[n+1,:] = sysResampling(samples, weights[n+1,:], N)

   for i in n:-1:1

       for j in 1:N

         sweights[i,j] = weights[i,j]*Mkernel(particles[i,j], sparticles[i+1,j], phi, sig)

       end

       ssumwei = sum(sweights[i,:])

       sweights[i,:] = sweights[i,:] ./ ssumwei

       samples = copy(particles[i,:])

       sparticles[i,:] = sysResampling(samples, sweights[i,:], N)

       println("Just finished smoothing at time $i")

   end
   #### end smoothing####

   ### calculate empirical filtering density
   temp1 = sparticles[2:n+1,:] .* sweights[2:n+1,:]

   for i in 1:n
       esmoothing[i] = sum(temp1[i,:])
   end


   return efilter, esmoothing

end





function sysResampling(particle, weight, N)


  x = zeros(N)

  c = zeros(N)

  u = zeros(N)

  c = cumsum(weight)

  u[1] = rand(Uniform(0,1)) / N

  i = 1

  for j in 1:N

      u[j] = u[1] + (j-1) / N

      while u[j]>c[i]
        i = i+1
       end

       x[j] = particle[i]
    end

    return x

end



2018年5月15日火曜日

修士留学について

最近ちょいちょい修士留学についての相談をなぜか受けるので、それについて適当に書いてみる。学部から留学している人とかは知らない。イギリスの大学院については実際に見ていること、アメリカについては話をきいた感じなので、その辺の区別をつけてくれるとハッピー。あと分野によってもばらけるので、その辺も考慮してほしいなぁ。所謂STEM分野だとこんな感じ位かな位の話で。

まずイギリスとアメリカの大学院の大きな違いは、イギリスでは修士課程と博士課程は基本的に分かれていて、かつ多くのイギリス人学生は修士号を取らないで学部卒で博士課程に進む。一方で留学生が博士課程に留学する時は修士号をすでに自分の国で取得していることが普通だと思う。博士課程自体は3年半から4年程で終わる。アメリカと違い、いわゆるオーバードクターは基本的に許されていないと思う。

アメリカの場合は基本的には修士課程と博士課程がくっついている5年制コースとなっている。またプレリムだとかそういう名称の進級試験があって、それに合格すれば3年に進級でき、出来なければクビとなる。まあなんでアメリカの大学で修士号ってのは残念賞扱いをする人も多い。最近は修士課程コースも出来ているようだけれど、修士号だけは博士課程をドロップアウト、クビになった証となるケースも多いので。基本的には5年制だと思うけれど、最近は長期化傾向にあるようで6年とか7年いる人もけっこういるようだ。

アメリカでもイギリスでも同じことが言えると思うけれど、日本と比較して修士号の価値は低いと思う。これは博士号を取得してアカデミアだけではなく、企業で働くのがキャリアプランとして普通のことなので、院=修士までという日本とはちょっと事情が違う。というかSTEM分野で博士号を取得したら企業に行く人の方が多いと感じる。アカデミアで働くより断然給与がいいし、分野によっては企業の方が強い自前の研究所を運営していることもある。

学費については、かなり高いと思う。1年制のフルタイムだと、例えば僕が所属する大学の学部だと年間で370万円程する。これに生活費なりを足すと、一年間で700万円程はかかると思う。アメリカだと例えばNYCにあるコロンビア大学だと、同じ分野で年間500万円ほどする。アメリカだと修士課程は2年間が普通だと思うので、1000万円以上のお金がかかる。また博士課程と違って修士留学を支援する奨学金などは皆無に等しいので、ほぼ自費留学となる。

次に講義のクオリティについて。これについては非常に低いと言わざるを得ないと思う。理由はまずイギリスもアメリカも修士コースというのはドル箱であって、なるべくたくさん学生をあつめて、金をかき集めようとしている。なので自然と先生と生徒の比率がとんでもないことになり、かつ学力の差が激しすぎるので、レベルも自然と低めになる。また先生も基本的に忙しいので修士コースの教育に熱心の人なんてめったに見かけないし、チューターとよばれる教育義務がある博士課程の学生が担当する講義も結構ある。個人的には日本の修士課程の方がしっかりしていると感じる。前述したとおり、博士課程に進学を望むようなイギリス人は修士課程に基本的には来ないので、そういうのもレベルの低下に拍車をかけていると思う。

だから学部で勉強してきた人ほど、イギリスの修士課程は物足りないと感じるんじゃないだろうか。とくにイギリスの修士課程は1年で終わるくせに、スッカスカの日程なので、講義があるのはコースにもよるけれど実質4か月位だと思う。

具体的にどういう生徒が多いかというと、金持ちの中国人の子供で埋め尽くされている。学生の8割位は中国系じゃないのかな。なんでこうなるかというと、イギリスの場合はEU圏か否かで学費が変り、EU圏外からの留学だと、圏内からの学生より3倍ほど高い学費を取る。なので大学としてはそういう学生に来て欲しいわけだ。

アメリカのことはしらないけれど、イギリスで修士をとった人の主なその後のキャリアは帰国だ。先にも述べたけれど、まず修士号の価値がそんなにたかくないので、高度な専門職にあんまりつけない。勿論修士号を取得して金融街とかで例えばクオンツとかで働いている人もいるけれど、国籍がイギリスであることが条件だと思う。それと多分だけれど、労働ビザももらえないと思う。

ならどういう人にお勧めできるかというと二つあって、一つは日本でいい仕事に就きたい人だと思う。キャリアアップ?とかそういうのだと思う。

二つ目はこっちの大学で修士号を取得して、博士にアプライする人とかかな。こっちの大学で学位を取得することで応募できる奨学金とかもあるし、推薦書も書いてもらえるかもしれないので。周りの友達にも、日本人ではないけれどそういう人は数名いる。まあどちらにせよ相当の支出があるので、コストをペイできるかは微妙だなーと思う。

なので1からバリバリ勉強したい人は、問答無用でアメリカの博士課程に行くのがいいと思う。もう基礎的な勉強はいいやって人は、イギリスの博士課程に行けばいいと思う。

まあでも取り合えず留学して、異国の地で色々と経験することは人生経験としては個人的には非常に有意義なことだと思う。こっちの修士課程に留学して能力がバリバリあがったり、こっちで仕事できるかといわれると、厳しいのかなと思う。

あと話はそれるけれど、日本よりいわゆる人文系の分野で修士号なり博士号をとっても、アカデミア以外での就職が本当に厳しい、需要がない。 

2018年4月18日水曜日

米国統計学大学院Phd受験の方法

・はじめに
受かったところ→イェール、ミネソタ、UCLA、シカゴ(ごねた)等
落ちたところ→コロラド州立大学、UCバークレー(二回も不合格通知を送ってきた悪の総本山)等

・GPA編
4点のスケールで3.5位が最低限。これ以下は素直にPhdは諦めましょう。来世があります。

・出願とそれまでのスケジュール
適当に勉強して遊びましょう。気が付けば成績がいいはずです。出願日のデッドラインを確認して出願しましょう。

・出願後の過ごし方
彼女がいないと留学は辛いです。彼女を作るか、いるなら大切にしましょう。

・英語(TOEFL)対策
気合でToefl ibtで100点以上取りましょう。統計学専攻は80点台とかが最低点の所もあるけれど、100点あった方がいい。90点が本当に最低ライン。出来なかったらPhdは諦めましょう。来世があります。

・GRE対策
Generalの算数のセクションは気合で満点とりましょう。GRE mathも満点近いスコアをとりましょう。 出来なかったらPhdは諦めましょう。来世があります。

・SoP  対策
そっれぽいのを気合でかきましょう。

 ・国内奨学金応募
沢山出して気合で何処か受かりましょう。授業料も生活費も5年位出してくれるとこが望ましいです。気合があれば何処か受かると思います。MA留学だと殆どの団体はお金をくれないです。何処かに受かっているとPhdの合格確率が間違いなくあがります。

・推薦状依頼
日頃から色々な先生と仲良くして、ここぞとばかりにお世話になりましょう。いい推薦書がもらえない人はPhdは諦めましょう。来世があります。

・出願先の決め方
適当

追記
・コネ
いい成績をとり先生と親しくなり推薦書を書いてもらえる間柄になるのも、出願先の大学の先生とコンタクトをとって押しかけるのも、すべてコネ。コネを作るのが一番有効で一番難しい。

2018年3月17日土曜日

Particle MCMCの簡単な解説とJuliaによる実装例

Andreiu et al. (2010) "Particle Markov chain Monte Carlo methods"を簡単に説明する.SMCを用いた隠れマルコフモデルの(ベイズ的な)推定で一番良く使われる手法で,名前の通りSMCとMCMCを組み合わせて事後分布からのサンプリングを行う.

xをlatent変数(マルコフ連鎖),yをデータ,cをパラメタとする.つまり隠れマルコフモデル(状態空間もでる)を考える.この時ポステリア,p(c|y)∝p(c)∫p(x|c)p(y|x,c)dxを推定したいとする.ここでp(c)はプライア.∫p(x|c)p(y|x,c)dx = p(y|c)なので,これは周辺尤度を求める必要があることを意味している.一般にlatent変数をintegrate outすることは難しいので,普通のMCMCだと(x,c)空間で探査することになり,xとcに相関がある場合など,うまくいかない,ミックスする速度が遅い等の問題が発生する.

ここでzというauxiliary(確率)変数を導入する.モンテカルロ推定値と思って構わない.z>0と仮定し,さらにzはp(y|c)の不偏推定量だとする.つまり,∫zp(z|c)dz = p(y|c)だとする.この時zp(z|c)p(c)をターゲットとしてMCMCを行えば,このターゲットはmarginalに∫zp(z|c)p(c)dz = p(c)p(y|c)となり,のぞみの分布(に比例したもの)が手に入る.このカーネルが満たすべきいくつかの条件を満たしていること等は簡単に示せる.また TVの意味で事後分布にあるレートで収束することも保障されている.

SMCの重要な性質として,ウェイトの平均がnormalizing constantの不偏推定量になることがある.つまり先ほどのZとして SMCで近似した周辺尤度を用いれば良いことがわかる.パラメタ更新のプロポーザル分布として対象な方法,例えばランダムウォークを用いたとすると,このときMH確率は min{1, A}, A:=(p(c*)Z*/p(c)/Z)となる,ここでZはSMCで近似した周辺尤度,*はプロポーザルを意味する.これをParticle independent Metropolis-Hastingsと呼ぶ.文脈によってはAndreiuの別の研究を踏まえで,Pseudo-Marginal MCMCなどとも呼ぶ.

アルゴリズムを簡単にまとめると,1:初期値のパラメタを所与としてSMCを回して尤度の不偏推定値を得る→2:この尤度を用いてMCMC(今回はMH法)でパラメータの更新を行う→3:更新されたパラメタのもとでSMCを回し尤度の不偏推定値を得る→4:この尤度を用いてMCMCでパラメータの更新を行う→5:3と4を好きなだけ繰り返す,となる.論文では他にギブス・サンプラー型のアルゴリズムも紹介されている.これらの総称をParticle MCMCという.

経験的にこの手法はパラメタ空間が低次元だとよく動く.ただ勿論数学的に示されるべきことは示されているけれど,高次元だと挙動が怪しくなってくる.あと当然バッチ推定なのでコストもかかる.ただコストがかかるといってもたかがしれてるのでCとかJulia等で実装すれば別に気になる程でもないと思う.寝る前に回しておくとかすればいい.

今回も簡単なSVモデルを例とし,真値として(phi, sigmx, beta) = (0.8, 0.36, 1.0)とセットしてn=1000個の観測値を生成した.SMCパートはブートストラップ・フィルタ+適応的リサンプリングを,MCMCパートは単純なランダム・ウォークMCMCを用いた.またパーティクルの数Nを500,MCMCの回数Mを10000とした.1000番目までの観測値をバーンイン期間とし,取り除いた.

何にも工夫しないランダムウォークMHを使ってるので,結果はあんまり良くないと思う.一応バーンイン後のパラメタのパスの平均はそれぞれ,0.78, 0.34, 1,14とまあ悪くはないけれど,とくにβは色々と怪しい.PMCMCのチューニングパラメタをどう選ぶべきかみたいな論文もあるけれど,興味があまりないので読んでいない.オリジナルの論文では5万回位MCMCパートを繰り返していたりするので,アイディアは素晴らしいけれどなかなか重い手法だと思う.ただ状態空間モデルのパラメタ推定は非常に大変なので,まあこんなものかなとも思う.

以下結果(ヒストグラム,スカッター,パラメタのパス,自己相関)とコード.なお以下のコードを使うのには前に紹介したブートストラップ・フィルタのコードが必要.






# PIMH (Andreiu et al.2010) for a simple SV model
# I use random walk MH
# x_t=Φ x_{t-1}+N(0,σ_x^2)
# y_t= β*exp(x_n/2)*N(0,1)
#
# You need the function Bootstrap_Filter to run SMC
#
# Input
#       data: observation
#       N: Number of particles
#       M: Number of MCMC steps
#
#
#  Output: path of phi, sigx beta and number of updates

function PIMH(data, N::Integer, M::Integer)

    @assert N > 0 "N must be a natural number"
    @assert M > 0 "M must be a natural number"
    n = length(data)
    par1 = zeros(M)
    par2 = zeros(M)
    par3 = zeros(M)
    likl = zeros(M)
    acpp = zeros(3)
    cur_state = zeros(4)

    ########################## Markov Kernel###################################

    function Mkernel(x,phi,sig)


        return(phi * x + rand(Normal(0, sig)))

    end
    ###########################################################################

 ########################Potential##########################################

    function Potential(x,y,beta)

     return(1.0/(sqrt(2.0*pi*beta^2*exp(x)))*exp(-(y^2)/(2.0*beta^2*exp(x))))

    end

 ###########################################################################




    ####initial step####

    #### initial values for parameters
    cur_state[1]= rand(Normal(0,1))
    phi = exp(cur_state[1]) / (1+exp(cur_state[1]))

    cur_state[2] = rand(Normal(0,1))
    sigma = exp(cur_state[2])*sqrt(1-cur_state[2]^2)

    cur_state[3] = rand(Normal(0,1))

 ####given initials, run SMC####
   output = Bootstrap_Filter(data, N, phi, sigma, cur_state[3])

   cur_state[4] = output[2,] # get log-likelihood

  ###end the first SMC step###


 ###start PMCMC ####

 for k in 1:M ### iterations for MCMC step


     ####propose candidates, run SMC, and run MCMC for par1

      pro_par1 =   rand(Normal(cur_state[1],5))

      phi = exp(pro_par1) / (1+exp(pro_par1))

      sigma = exp(cur_state[2])*sqrt(1-phi^2)

      result = Bootstrap_Filter(data, N,  phi, sigma, cur_state[3])

      prop_likl = result[2,]

      newprior1 = sum(logpdf(Normal(0,1),[pro_par1, cur_state[2], cur_state[3]]))

      oldprior1 = sum(logpdf(Normal(0,1),[cur_state[1], cur_state[2], cur_state[3]]))

      acp = exp((prop_likl - cur_state[4] + newprior1 - oldprior1))

      u = rand(Uniform(0,1))

      if u < acp
          cur_state[1] = pro_par1

          cur_state[4] = prop_likl

          acpp[1] = acpp[1] + 1

      end

      ####propose candidates, run SMC, and run MCMC for par2

      pro_par2 =  rand(Normal(cur_state[2] ,5))

      phi = exp(cur_state[1]) / (1+exp(cur_state[1]))

      sigma = exp(pro_par2)*sqrt(1-phi^2)

      result = Bootstrap_Filter(data, N,  phi, sigma, cur_state[3])

      prop_likl = result[2,]

      newprior2 = sum(logpdf(Normal(0,1),[cur_state[1], pro_par2, cur_state[3]]))

      oldprior2 = sum(logpdf(Normal(0,1),[cur_state[1], cur_state[2], cur_state[3]]))

      acp = exp((prop_likl - cur_state[4] + newprior2 - oldprior2))

      u = rand(Uniform(0,1))

      if u < acp
          cur_state[2] = pro_par2

          cur_state[4] = prop_likl

          acpp[2] = acpp[2] + 1
      end



     ####propose candidates, run SMC, and run MCMC for par3

      pro_par3 =  rand(Normal(cur_state[3],5))

      phi = exp(cur_state[1]) / (1+exp(cur_state[1]))

      sigma = exp(cur_state[2])*sqrt(1-phi^2)

      result = Bootstrap_Filter(data, N,  phi, sigma, pro_par3)

      prop_likl = result[2,]

      newprior3 = sum(logpdf(Normal(0,1),[cur_state[1], cur_state[2], pro_par3]))

      oldprior3 = sum(logpdf(Normal(0,1),[cur_state[1], cur_state[2], cur_state[3]]))

      acp = exp((prop_likl - cur_state[4] + newprior3 - oldprior3))

      u = rand(Uniform(0,1))

      if u < acp

          cur_state[3] = pro_par3

          cur_state[4] = prop_likl

          acpp[3] = acpp[3] + 1
       end

      par1[k] = cur_state[1]
      par2[k] = cur_state[2]
      par3[k] = cur_state[3]
      likl[k] = cur_state[4]


     #####end PMCMC
     println("Just finished time $k")
   end

    Phi = exp.(par1)./(1+exp.(par1))

    Sig = exp.(par2).*(1-Phi.^2)

    Beta = par3

 ##end function
     return Phi, Sig, Beta, acpp
end


2018年3月13日火曜日

Juliaで逐次モンテカルロ

Juliaで書いたSV modelへのSMCのコード。シンプルなブートストラップフィルタ+適応的リサンプリングを採用。使うためにはDistributionsパッケージが必要。今回は(phi, sig, beta) = (0.91, 1.0, 0.5)とした。サンプルサイズとパーティクルの数はそれぞれ1000とした。大体の環境で1秒もかからず計算終わるはず。以下結果とコード。







# Bootstrap Filter (Gorodn et al.1993) for a simple SV model
# I use adaptive resampling
# x_t=phi x_{t-1}+N(0,sigma^2)
# y_t= beta*exp(x_n/2)*N(0,1)
#
# Input data: observation
#       N: Number of particles
#       phi: paramater
#       sig: paramete r
#       beta: parameter
#
#  Output: Effective sample size, log-likleihood, posterior mean and variance.


function Bootstrap_Filter(data, N::Integer, phi::Real, sig::Real, beta::Real)

    @assert N > 0 "N must be a natural number"
    @assert phi < 1
    @assert phi > -1
    n = length(data)

    samples = zeros((n+1),N)
    weights = zeros((n+1),N)
    th = N / 2
    pm = zeros(n,1)
    pv = zeros(n,1)
    logl = 0.0
    ESS = zeros(n)



    ########################## Markov Kernel###################################

    function Mkernel(x,phi,sig)


        return  phi * x + rand(Normal(0, sig))

    end
    ###########################################################################

    ########################Potential##########################################

    function Potential(x,y,beta)

        return (1.0/(2.0*pi*beta^2*exp(x)))*exp(-(y^2)/(2.0*beta^2*exp(x)))

    end

    ##########################################################################

    #~~initial step
    for j in 1:N

        samples[1,j] = rand(Normal(0,sig/(1-phi)))

        weights[1,j] = 1.0 / N

    end
    ###

    logl = log(sum(weights[1,:]) / N)

    ###main loop

    for i in 1:n

        for j in 1:N

            samples[i+1,j] = Mkernel(samples[i,j], phi, sig)

            weights[i+1,j] = log(weights[i,j]) + log(Potential(samples[i+1,j], data[i], beta))

        end




            ###normalize weights###

            maxwei = maximum(weights[i+1,:])

            for j in 1:N
                weights[i+1,j] = exp(weights[i+1,j]-maxwei)
            end

            ### calculate log-lilelihood

            logl = logl + log(sum(weights[i+1,:])/N)

            sumwei = sum(weights[i+1,:])

            for j in 1:N
                weights[i+1,j] = weights[i+1,j] / sumwei
            end

            ####

            ### calculate posterior mean and variance

            pm[i] = weights[i+1,:]'*samples[i+1,:]


            pv[i] = weights[i+1,:]'*samples[i+1,:].^2  - pm[i]^2


            ### check resampling

            sumsqwei = weights[i+1,:]' * weights[i+1,:]

            ESS[i] = 1.0 / sumsqwei

            if ESS[i] < th

                copy_samples = copy(samples[i+1,:])

                index = wsample(1:N, weights[i+1,:], N)

                samples[i+1,:] = copy_samples[index]

                weights[i+1,:] = 1.0 / N

            else

                samples[i+1,:] = samples[i+1,:]

                weights[i+1,:] = weights[i+1,:]

            end





             ### end resampling





          println("Just finished time $i")




  end


     return ESS,logl, pm, pv


end

2018年3月2日金曜日

エビデンスに基づいた政策運営

エビデンスに基づいて政策運営をしましょう,みたいなフレーズをよく目にする.例えば東大も,政策評価研究教育センターを立ち上げた.僕は政治に興味が持てないし,ポリコレだ不平等だ貧困問題だどうだを真剣に話している人達は苦手というか,友達になれないというか,そういう人達を茶化してしまう人なので,まあ正直どうでもいいのだけれど,思うことをつらつら書いてみたい.ここでいうエビデンスとは統計的手法にもとづいてなんらかの行為の識別された因果を一応指す.

1 まずエビデンスとやらに基づいて政策運営をして,そうでない政策運営よりよくなったエビデンスがあるのかなと思う.よくなるの基準が曖昧なので,まあなんらかの意味ででいいのだけれど.こういうことを言うと延々と遡って同じことがいえるので不毛といえば不毛かもしれないけれど,なんでそう統計的手法なり基づいて何かすれば裁量的な運営よりよくなると自信をもっていえるのかがわからない.

2 これはちょいちょいってるけれど,同じことが2度と起こらない,例えば政策みたいなもの効果の因果関係が仮に識別できたとして,それがそこまで意味があることなのかなと思う.例えば10年前にある所得以下の人に給付金をばらまく政策をしたとして,その効果の因果関係がわかったとする.けれどこの分析結果をもって,今同じ政策をしたら10年前よ同じような結果が得られるとは言えない.これは時系列データへの解析と同じ考えで,ある時点での事象はその時一回きりしか観測できないので,何かしらデータを生成している過程に仮定をおかないと,統計解析が意味をもたない.時系列解析の場合はエルゴード性だとかを使って正当化するのだけれど,因果推定の文脈ではこれをどう正当化しているのかが分からない.フィルタリングの話で言えば,因果推定はスムージングとしては機能するけれど,予測等には使えないように思う.ある政策Aは時点nでC位影響があった,なかったとか,できてその程度の議論じゃないだろうか.

3 2のようなことを言ってみたところ,専門的に研究されている人から,だから経済構造を仮定して分析を行うんだと言われ.まあ確かにそうだようなぁと思ったのだけれど,構造を仮定して分析すればすればで別な問題もつきまとう.同じことを何度も言ってバカみたいだけれど,なんでその経済構造の仮定が正しいよいえるのかとか色々.あとまあ基本的にceteris paribusで分析するので,短期的な話ならある程度説得力もあるのだけれど,長期的な影響の話になると正直とても胡散臭い,

4 この手の主張を強くして,キャンペーンをはる人達の中に,とても都合の良い人だなぁと思ってしまう方が何人もいる.というか己の価値観なり社会性議論ありきで,何かしている人が多いと感じてしまう.日頃は因果と相関の違いをくちうるさく述べているのに,自分の主義主張に合う研究だったらコホート分析とかそういった類のものまでエビデンスとして紹介してしまうような人たちだ.学者として提言するならば,己の主義主張はひとまず置いておいてほしいと思う.何をエビデンスとするのかは本人の自由なので好きにすればいいと思うのだけれど,きちんとエビデンスの定義なりを共通しておかないと,エビデンスに基づいた政策運営なりは非常に不毛だと思う.本当にはげているのなら申し訳ないのだけれど,せめて議論くらいには毛を生やしてほしい.

5 例えば優生学的なものも,データを使えばある程度は正しさを示せると思う.ただデータ的に正しいからといって,こういったものを正しいと言ってしまうことの是非はまた別な問題として扱われるべきだと思う.倫理的な話.

6 エビデンスって言葉がバカっぽい.

2018年2月10日土曜日

Bootstrap Filter for a simple SV model (Gordon et al. 1993)




/*
*bootstrap filter (Gordon et al. 1993)  with adaptive resampling
*x_t = rho*x_t-1 + sigma*N(0,1)
*y_t=beta*exp(x_t/2)*N(0,1)
*/

/*outputs = (particles, normilized weighsts, posterior)*/

#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <math.h>
#include <time.h>

gsl_rng * r;

double * data_y;

void initial_step (double * x0, double * wei);

void bootstrap_filter (double *x0, double *x1, double *wei0, double *wei1, double datum, double rho, double phi, double beta);

void Mkernel (double *x0, double *x1, double rho, double phi);

double potential (double *x, double y, double beta);

void resample( double *samples, double *weights, int N);

///

int main (void){

/********** Initializing the Random Number Generator *********/

const gsl_rng_type * T;
gsl_rng_env_setup();
T = gsl_rng_default;
r = gsl_rng_alloc (T);
printf("give the seed: ");
unsigned long int seed; scanf("%lu", &seed);
gsl_rng_set (r,seed);

/*************************************************************/

/* n includes the first step */
int n;
printf("\nGive the time-steps:\n"); scanf("%d", &n);

/* number of particles*/
int N;
printf("\nput number of particles: \n"); scanf("%d", &N);


/*rho*/
double rho;
printf("\nrho?: \n"); scanf("%lf", &rho);

/*phi*/
double sigma;
printf("\nsigma?: \n"); scanf("%lf", &sigma);

/*beta*/
double beta;
printf("\nbeta?: \n"); scanf("%lf", &beta);


int i,j;
//////////////////////////////////////////////////////////////////

double ** samples = malloc ((n+1) * sizeof (double*));

for (i = 0; i < (n+1); i++) samples[i] = malloc ( N * sizeof(double));

double *posterior = malloc (n * sizeof (double));

/* load the data */

  FILE * read_data;

  read_data = fopen("data.txt","r");

  data_y = malloc (n * sizeof(double));

  for ( i = 0; i < n ; i++ ) fscanf(read_data ,"%lf", data_y+i);

////////////////////////////////////////////////////////////////

FILE * save_samples = fopen("samples.txt","w");

FILE * save_weights = fopen("weights.txt","w");

FILE * save_posterior = fopen("posterior.txt","w");

double ** weights = malloc( (n+1) * sizeof(double *));

for ( i = 0; i < (n+1) ; i++ ) weights[i] = malloc( N * sizeof(double) );

/////////////////////////////////////////////////////////////////////////


/*initial step*/

for (i = 0; i < N; i++) initial_step (&(samples[0][i]), &(weights[0][i]));

for ( i = 0; i < N ; i++ ) fprintf( save_weights, "%15.9lf ", weights[0][i]);

for ( i = 0; i < N ; i++ ) fprintf( save_samples, "%15.9lf ", samples[0][i]);

fprintf( save_weights, "\n");

fprintf( save_samples, "\n");

//////////////////////////////////////////////////////////////////////////

/*main loop*/

for (i = 0; i < n; i++){

  for(j = 0; j < N; j++){

bootstrap_filter (&(samples[i][j]), &(samples[i+1][j]), &(weights[i][j]), &(weights[i+1][j]), data_y[i], rho, sigma, beta);

 }

  //Normalizing and Resampling Check////////////////////

  double sumwei = 0.0;

  double sumsqwei = 0.0;

  double empiricaldist = 0.0;

  for (j = 0; j < N; j++) weights[i+1][j] = exp(weights[i+1][j]);

  for (j = 0; j < N; j++) sumwei += weights[i+1][j];

  for (j = 0; j < N; j++) weights[i+1][j] = weights[i+1][j] / sumwei;

  for (j = 0; j < N; j++) sumsqwei += (weights[i+1][j]*weights[i+1][j]);

  for (j = 0; j < N; j++) empiricaldist += samples[i+1][j]*weights[i+1][j];

// calculate posterior before resampling
  posterior[i] += posterior[i] + empiricaldist;

 fprintf( save_posterior, "%15.9lf ", posterior[i]);

 fprintf( save_posterior, "\n");

 double ESS = 1.0 / sumsqwei;

  if (ESS < N / 2) {

  resample (samples[i+1], weights[i+1], N);

  //reset weighsts
  for (j = 0; j < N; j++) weights [i+1][j] = 1.0 / N;

  }

/////////////////////////////////////////////////////////////////

  for (j = 0; j < N; j++){

    fprintf( save_weights, "%15.9lf ", weights[i+1][j]);

    fprintf( save_samples, "%15.9lf ", samples[i+1][j]);

  }


  ///////////////////////////////////////////////////////

  printf("Just finished time t=%d\n", i+1);


  fprintf( save_weights, "\n");
  fprintf( save_weights, "\n");

}

////////// Finished with Execution

fclose(save_samples);

for (i = 0 ; i < (n+1); i++){

  free(weights[i]);

  free(samples[i]);

}

free(weights), free(samples), free(posterior);


return(0);
}

////////////////////////////initial_step////////////////////////////////////////

void initial_step (double *x0, double *wei){

  *wei = 1.0;
  *x0 = gsl_ran_gaussian(r,1.0);
}


////////////////////////////////////////////////////////////////////////////////

////////////////////////////bootstrap_filter//////////////////////////////////////////////////////////////////////////////////////////////////////////

void bootstrap_filter (double *x0, double *x1, double *wei0, double *wei1, double datum, double rho, double phi, double beta){



/// probagate and mutate///////////////////

 Mkernel (x0, x1, rho, phi);

*wei1 = *wei0 + log(potential (x1, datum, beta));

//////////////////////////////////////////////



}

/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

///////////////////////Mkernel/////////////////////////////////////////////////

void Mkernel (double *x0, double *x1, double rho, double sigma){

*x1 = *x0*rho +sigma*gsl_ran_gaussian(r,1.0);

}

///////////////////////////////////////////////////////////////////////////////


///////////////////////potential///////////////////////////////////////////////

double potential (double *x, double y, double beta){

return ( exp( -y*y/( 2.0*beta*beta*exp(*x) ) ) );

}


///////////////////////////////////////////////////////////////////////////////


/////////////////////// resample//////////////////////////////////////////////

void resample( double *samples, double *weights, int N) {

  int i;

  double *copy_samples = malloc( N * sizeof(double) );

  for ( i = 0; i < N ; i++ ) copy_samples[i] = samples[i];

  gsl_ran_discrete_t * gsl_weights = gsl_ran_discrete_preproc(N, weights);

  for ( i = 0; i < N ; i++ )  {

    size_t index = gsl_ran_discrete(r, gsl_weights);

    samples[i] = copy_samples[index];

  }

  gsl_ran_discrete_free(gsl_weights);

  free(copy_samples);

}

///////////////////////////////////////////////////////////////////////////////