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

}

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

2018年2月9日金曜日

計算機統計学周辺の研究をする人向けの講義ノートや解説論文

昨日に引き続き、表題について備忘録もかねて紹介する。当然内容は偏っている。コメント付き。

・逐次モンテカルロ
Doucet and Johansen "A Tutorial on Particle Filtering and Smoothing: Fifteen years later"  逐次モンテカルロについての一番わかりやすいチュートリアル論文。理論面はそこまで深くないけれど、この分野を基礎を知るのによい。Doucetは女好き、 Johansenとはロン毛仲間。二人ともナイスガイ。
Del Moral and Doucet "Particle methods: An introduction with application" SMCの数理的な部分のチュートリアル。ある程度この分野になれてないと多分読めない。Del Moralも女性の話を大体している。ナイスガイ。
Pollock "“Introduction to Particle Filtering” Discussion - Notes " 多分彼が学生の頃にかいたタームペーパーか何かだと思うけれど、SMCへのよい入門論文だと思う。
Kantas et al. "On Particle Methods for Parameter Estimation in State-Space Models " SMCを用いたパラメタ推定についてのサーベイ論文。
Sequential Monte Carlo Methods & Particle Filters Resources Doucetが管理しているサイト。SMCの文献が網羅的に乗せられていて便利。
Modern Computational Statistics 今イギリスにいる統計学者で最も生産的なFearnheadの講義サイト。SMCだけじゃなくてMCMCやABCの紹介も。Rのコード付き。

・MCMC
Roberts and Rosenthal "General state space Markov chains and MCMC algorithms" MCMCの理論を鍬もって開拓し続けている二人による理論面に焦点を置いたチュートリアル論文。
Geyer "Markov Chain Monte Carlo Lecture Notes" 講義ノート。測度論的表記に慣れるために良いと思う。
Sherlock "Reversible Markov chains" これも彼の講義ノート。ディリクレ形式を用いて漸近分散を評価する勉強に。
Neal "MCMC using Hamiltonian dynamics" とても有名なハミルトニアン・モンテカルロのご本人によるチュートリアル論文。
Saloff-Coste "Lectures on finite Markov chains "  state spaceが有限の時だけだけれど、色々なテクニックの解説が載っていて勉強になった。
Bierkens "Non-reversible Metropolis-Hastings" 普通の論文なのだけれど、いわゆる詳細つり合い条件を満たさないNon-reversibleなマルコフ連鎖を用いるのがこの分野の最近の流行りなので、その入門に。

・ABC
Blum et al. "A Comparative Review of Dimension Reduction Methods in Approximate Bayesian Computation" 所謂summary statisticsに焦点をおいたサーベイ論文。
Prangle "Summary Statistics in Approximate Bayesian Computation" これも↑と扱っている内容は同じ。
Fearnhead "Asymptotics of ABC" ABCの漸近的性質。この分野を鍬もって開拓しているポールによる解説。秋位にでる教科書に掲載されるらしい。
Marin et al. "Approximate Bayesian Computational methods" ABCそのものについてのサーベイろんぶん。簡潔に書かれていて入門に良いと思う。

・NPB
Orbanz "Lecture Notes on Bayesian Nonparametrics" そのまんまノンパラベイズの講義ノート。最近出たファン・デル・ファールトの教科書を読むのがめんどくさい僕みたいな人向けの、理論的側面に焦点をあてたノート。

・変分ベイズ
Blei et al "Variational Inference: A Review for Statisticians" この分野を精力的に研究している統計学者たちによるノート。どのあたりが統計学者のためなのかは分からなかったけれど、Stochastic variational inferencの話題もあって勉強になった。
Ormerod and Wand 2010 "Explaining Variational Approximations" ↑をもっと簡素化したような感じ。さっと勉強したい人むけだと思う。
Beal "Variational Algorithms for Approximate Bayesian Inference" とても有名な彼の学位論文。本当に色々なことが書かれているが、彼はアカデミアをさりファンド・マネージャーへ。

・その他
Scalable statistical inference 夏にケンブリッジであったこの分野の研究会。最先端の報告が盛りだくさんなので、今何が流行ってるのかを勉強したい人はぜひ。

2018年2月8日木曜日

計算機統計学周辺の研究をする人向けの教科書

この分野の論文を読んだり、書いたりするときに、自分が読んでよかったり持っていて役に立っている教科書をリストアップする。コメント付き。当然話題は偏っている。CかC++等でプログラミングできることも必須。Rはやっぱり遅い!お絵かきにはとても良い。Juliaもいいらしいけれど、いかんせんセクシーな方のジュリアが検索で引っかかるのが難点か。

・測度論・確率論
Capinski and Koop "Measure, Integral and Probability" 邦訳があれば解答とかあるのでそっちの方がいい。これくらいの測度論の知識がないと日常生活すら送るのが困難なのでは。
Ash "Probability and Measure Theory" 関数解析の話もあっていい。
Meyn and Tweedie "Markov Chains and Stochastic Stability" マルコフ連鎖の知識がないとこの分野で生きてけない。とても分かりにくい一冊だけれど、これしかない。ディリクレ形式の話がないのが不満。
Hall and Hyde "Martingale Limit Theory and Its Application" マルチンゲールでぶん殴るための一冊。便利。ノーベル・マルチンゲール賞をあげたい。

・関数解析・確率解析
山田 "工学のための関数解析" MCMCがらみで必要になってくるスペクトルの話が余りない以外は十分。
Kreyszig "Introductory functional analysis with applications" 分厚いけれどわかりやすい。
Bobrowski "Functional Analysis for Probability and Stochastic Processes" いわゆるマルコフ・オペレータのお勉強に使った。わかりやすい。確率論とのつながりを意識した関数解析の教科書。
Oksendal "Stochastic Differential Equations" 結構必要になるSDEの知識の会得。
Iacus "Simulation and Inference for Stochastic Differential Equations" 実際にどうシミュレートするか。↑には載ってない小ネタがのってたり、参考文献さがしにもいいし、結構役に立つ一冊。
Kloeden et al "Numerical solution of SDE through computer experiments" SDEの数値解析のバイブル。
Rao "Statistical Inference for Fractional Diffusion Processes" 記憶を持つSDEの理論的な話が載ってる。


・統計学
Cappe et al. "Inference in hidden markov models" 一般的な隠れマルコフモデル(状態空間モデル)の数理的なことが書かれている。分かりやすくてとてもいい。
Douc et al. "Nonlinear time series" ↑と同じ内容をちょっと新しい話題までカバーしたかんじ。↑より簡潔に書かれているのでちょっと分かりにくい。
Bernardo and Smith "Bayesian Theory" ベイズ統計学で一冊と言ったらこれだと思う。意思決定論的なアプローチで最初の方は退屈なので読んでない。

・計算機統計学
Liu "Monte Carlo Strategies in Scientific Computing" 色々なモンテカルロ法の話が載っている。けど分かりやすく書こうとし過ぎて逆になんか分かりにくい。
Robert and Casella "Monte Carlo Statistical Methods" この分野のバイブルと言えばバイブルかもしれないけれど、そう呼ぶにはちょっと頼りない感じ。けどモンテカルロ法に興味がある人はマストバイ。
Efron and Hastie "Computer Age Statistical Inference: Algorithms, Evidence, and Data Science" 普通の尤度解析の話からSVMやらランダム・フォレストやらニューラルネットワークの話まで、広く薄くカバーしていて結構いい。
Sisson et al. "Handbook of Approximate Bayesian Computation" 一部で微妙に流行っているABCの唯一の教科書。まだ出版されてないけれど、arxiv等で一部掲載される論文が見れる。
Del Moral "Feynman-Kac Formulae: Genealogical and Interacting Particle Systems with Applications" 色々なモンテカルロ法の数学的基礎が統一的に書かれている。バイブル。
Sarkka "Bayesian Filtering and Smoothing" スムージングやフィルタリングの基礎が書かれている。カルマンフィルタや逐次モンテカルロについても。読むのはめちゃめちゃ簡単。
Rasmussen and Williams"Gaussian Processes for Machine Learning" 著者の一人に貰った一冊。何かと必要になるガウス過程の基礎やそれを用いた統一的な機械学習へのアプローチが分かりやすく解説されててよかった。