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" 著者の一人に貰った一冊。何かと必要になるガウス過程の基礎やそれを用いた統一的な機械学習へのアプローチが分かりやすく解説されててよかった。

2018年2月1日木曜日

SMCサンプラー

#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>
/*
 *
 */


/* The initial distribution for each co-ordinate*/
double eta_1( void );

/* The potential function */
double potential( double x, double y );

/* The Markov kernel samples */
double M_kernel_sample( double x );


/* Retrieving the samples at the beginning of the trajectory */
void retrieve_back_samples( double * back_samples, double ** samples, int ** ancestory, int N, int dim, double * final_G );

/* A complete step over the first time instance t=0 */
void smc_propagate_0( double ** samples, int ** ancestory, int N, int dim, double ** final_G );

/* One SMC step along space */
void smc_propagate ( double ** samples, int ** ancestory, int N, int dim, double * final_G, double data_point );

/* The function over which to produce estimates */

double f( double x );


double ** data_y;


gsl_rng * r;


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

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

  int dim;

  printf("\nGive the dimension:\n"); scanf("%d", &dim);

  int n;

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

  int i, j;

  /* load the data */

  FILE * read_data;

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

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

  for ( i = 0; i < dim ; i++ ) data_y[i] = malloc ( n * sizeof(double) );

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

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

  }

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

  int N;

  printf("\nGive the number of particles:\n"); scanf("%d", &N);

  double ** samples = malloc( dim * sizeof(double *));

  int ** ancestory = malloc( dim * sizeof(int *) );

  int curr_coordinate;

  int i_space; int i_time;

  FILE * save_samples;

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

  int complete_samples;

  printf("\nPress (1) for complete samples, press (O) for averages :\n"); scanf("%d", &complete_samples);

  //// We execute time t=0 first;
  ////

  /* I follow the Del Moral "notation", thus onced finished all dim's at step 1
     I have NOT resampled according to the final weight, I have only weighted */

  double * final_G;

  smc_propagate_0( samples, ancestory, N, dim, &final_G );

  ////

  //// I return samples over the first co-ordinate;
  //// Recall, we need one more resampling, and need to retrieve the path backwards

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

  retrieve_back_samples( back_samples, samples, ancestory, N, dim, final_G );

  if ( complete_samples ) {

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

    fprintf( save_samples, "\n"); }

  else {

    double estimate = 0.0;

    for ( i = 0; i < N ; i++ ) estimate = estimate + f(back_samples[i]);

    estimate = estimate/N;

    fprintf( save_samples, "%15.9lf \n", estimate);  }


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

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

  //// We then proceed to the other times;

  for ( i_time = 1 ; i_time < n ; i_time++ )  {

    for ( i_space = 0 ; i_space < dim ; i_space++ )  smc_propagate ( samples, ancestory, N, dim, final_G, data_y[i_space][i_time] );

    retrieve_back_samples( back_samples, samples, ancestory, N, dim, final_G );

    if ( complete_samples ) {

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

      fprintf( save_samples, "\n"); }

    else {

      double estimate = 0.0;

      for ( i = 0; i < N ; i++ ) estimate = estimate + f(back_samples[i]);

      estimate = estimate/N;

      fprintf( save_samples, "%9.5lf \n", estimate);  }

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

  }

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


  fclose(save_samples);


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


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

    free( samples[i] ); free( ancestory[i] ); free ( data_y[i] );    }

  free(samples); free(ancestory); free( data_y );


  return(0);

}



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


void smc_propagate_0( double ** samples, int ** ancestory, int N, int dim, double ** final_G ) {

  int i_samples = 0; int i_space = 0; int k;

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

  int * curr_ancestor = malloc( N * sizeof(int) );

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


  ////////// We fill in the 0-th co-ordinate


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

    curr_samples[i_samples] = eta_1();

    curr_G[i_samples] = potential(curr_samples[i_samples], data_y[0][0]);

    curr_ancestor[i_samples] = i_samples;                    }


  samples[0] = curr_samples;  ancestory[0] = curr_ancestor;


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


  ////////// We fill in the remaining co-ordinates

  double * new_samples;

  int * new_ancestors;


  for ( i_space = 1 ; i_space < dim ; i_space++ ) {

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

    new_ancestors = malloc( N * sizeof(int) );

    /* Resample and Propagate */

    gsl_ran_discrete_t * normalise = gsl_ran_discrete_preproc (N, curr_G);

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

      size_t index =  gsl_ran_discrete (r, normalise);

      new_samples[i_samples] = eta_1();

      new_ancestors[i_samples] = index;

      curr_G[i_samples] = potential(new_samples[i_samples], data_y[i_space][0] );  }

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

    /* We must now shift the matrices one position down */

    for ( k = i_space ; k > 0 ; k-- ) { samples[k] = samples[k-1];  ancestory[k] = ancestory[k-1]; }

    samples[0] = new_samples; ancestory[0] = new_ancestors;

    gsl_ran_discrete_free (normalise);

  }

  (*final_G) = curr_G;

}


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


void smc_propagate ( double ** samples, int ** ancestory, int N, int dim, double * final_G, double data_point ) {

  /* Resample and Propagate*/

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

  int * new_ancestors = malloc( N * sizeof(int) );

  int i_samples; int i_back; int k;

  gsl_ran_discrete_t * normalise = gsl_ran_discrete_preproc (N, final_G);

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

    size_t index =  gsl_ran_discrete (r, normalise);

    /* To propagate, I now need to retrieve the ancestor of the index-particle d-1 steps earlier. */

    int back_index = index; double back_sample;

    for ( i_back = 0; i_back < dim-1 ; i_back++ ) back_index = ancestory[i_back][back_index];

    back_sample = samples[dim-1][back_index];

    new_samples[i_samples] = M_kernel_sample(back_sample);

    new_ancestors[i_samples] = index;

    final_G[i_samples] = potential(new_samples[i_samples], data_point);

  }

  gsl_ran_discrete_free (normalise);

  /* We must now shift the matrices one position down */

  free( samples[dim-1] ); free( ancestory[dim-1]);

  for ( k = dim-1 ; k > 0 ; k-- ) { samples[k] = samples[k-1];  ancestory[k] = ancestory[k-1]; }

  samples[0] = new_samples; ancestory[0] = new_ancestors;

}


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


void retrieve_back_samples( double * back_samples, double ** samples, int ** ancestory, int N, int dim, double * final_G ) {

  gsl_ran_discrete_t * normalise = gsl_ran_discrete_preproc (N, final_G);

  int i_back; int i;

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

    size_t index =  gsl_ran_discrete (r, normalise);

    /* I now need to retrieve the ancestor of the index-particle d-1 steps earlier. */

    int back_index = index;

    for ( i_back = 0; i_back < dim-1 ; i_back++ ) back_index = ancestory[i_back][back_index];

    back_samples[i] = samples[dim-1][back_index];

  }

}


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


double eta_1( void ) {

  return( gsl_ran_gaussian(r,1.0) );  }


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


double potential( double x, double y ) {


//  return(exp(-x*x/2));

  return ( exp( -(y-sin(x))*(y-sin(x))/(2.0*0.4*0.4) ) );

}


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


double M_kernel_sample( double x ) {

  /* I have chosen parameter values so that the stationary law is N(0,1) */
  double y = 0.5*x + sqrt(3.0)/2.0*gsl_ran_gaussian(r, 1.0);

  return(y);  }


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


double f( double x ) {  return(x);  }

2017年12月3日日曜日

イギリスとアメリカのPh.dの出願やシステムの違い等

この手の話はもう結構情報が共有されていると思うけれど、アメリカとイギリスの大学両方に出願した人間は珍しいかもしれないので、誰かの役に立てばいいなという感じで書いてみる。僕は統計学部にしか出願していないので、他の学部については詳しくは知らないけれど、まあどこもそこまで変わらないと思う。進学先はイギリスなので、アメリカの実際の博士課程での生活は知らない、体験してないので。大学名はいちいち書かないけれど、出願結果の打率は7割程度だった。三割ちょっとしか打てないイチローがその時は小さく見えました。

1 イギリスの博士課程の出願について。求められた書類は、成績表と、IELTSのスコアと、推薦書位だった。後で書くGREのスコアは任意だし、そこまで重要視されていないと思う。IELTSは6.5から7.5位が求められているレンジで、各セクション毎に最低点が課せられている。ただ英語のスコアが足りてなくても、条件付き合格をくれて、期間内にスコアが出れば正式な合格をくれるみたいな制度もあって、英語に関しては結構融通が効く。所謂GPAの明確な足切りラインみたいなのも無かった。まあでも暗黙の了解としてみんな出身大学ではトップ層だと思う。

イギリスの博士課程の大きな特徴は、日本と同じで修士と博士が完全に分かれているところにあるとおもう。経済学部だとアメリカ式の5年一貫の博士課程コースもあるみたいだけれど、これは例外。大学にもよるけれど、大体3年半から4年で博士課程が終わる。日本と違ってオーバードクターが延々出来るわけではなく、ある程度の期間で取得しないとキックアウトされる。これはオーバードクターしている学生が多いと、国からもらえる補助金が減らされるためだからだと聞いたけれど、定かではない。

またイギリスだと学部が3年、修士が1年で終わるので、日本の学生が学士を取るのと同じ期間で修士号を取れる。博士課程の出願に修士号は必ずしも必要とされないけれど、外国人は基本的にどっかの国で修士号を取得している。なので企業に就職したことが無く、学部から博士課程にそのまま進学してくるイギリス人の学生は22歳前後とかなので、兎に角若い。また修士課程も日本ほどありがたがられていないというか、ビジネス色が強いので、レベルはかなり低いと感じる。これは色々な国から来てる留学生も同じことを皆言ってるので、世界的に見てもやってる内容のレベルは低いと思う。生の人参を講義中にかじり始める傾奇者な学生もいた。

出願期間は割と長い期間あって、定員が埋まり次第その年の募集は終了する。合否の決断もまず習いたい先生にコンタクトをとって、その先生が書類などを見てインタビューをして、インタビューの結果、指導してもいいと思ったら学部レベルに話をあげる感じだった。なので入試委員会があるとかではなく、合否が先生毎に委ねられている感じだと思う。はっきり言ってこっちの先生とコネがあれば比較的簡単に合格できると思う。コネを作るのが難しいのだけれど。僕はコネつくりの結果、ゴマをすりすぎて指紋が溶けてなくなりました。

一番のネックは、大学の奨学金制度が全然充実していない所だと思う。学部は基本的に生活費や学費を出してくれない。EU圏の学生向けのは結構あるけれど、それ以外の学生向けのは殆どないので、基本的には自前で金策する必要がある。学部によっては研究室の先生が資金を持っている場合は、それから出してくれることもあるみたい。ただイギリスで修士号を取得すると、イギリスで応募できる奨学金の幅が色がるので、そういう理由でこっちで修士号を取る人もいる。Cash Rules Everything Around Me.

分野にもよるけれど、博士号を取得して企業に就職するキャリアが普通なことと、博士取得にかかる年数が比較的短いので、修士号の価値は相対的に低い。こっちで修士号だけ取得して、こっちの有名企業とかに就職できるとかは基本的には夢物語だと思う。学費も高いので、こっちで修士号だけ取得して得られるペイが、コストを上回るかと言われると、正直微妙だと感じる。履歴書を飾っていけ。


2 アメリカの博士課程の出願について。アメリカの大学への出願は兎に角大変だった。求められたのは成績表、Toefl ibtのスコア、推薦書、SOP、GRE generalとMathのスコアだった。まずイギリスの大学と違って、各スコアに明確な足切りラインがあることが多い。Toeflなら90点前後、GRE mathなら上位5%、GPAなら3.5/4.0とかそんな感じで明確に書かれていることが多い。大学によってはToeflの各セクションに最低点が設定されていることもあった。ちなみに僕はイギリスにきて4歳児が英語ペラペラで心が折れそうでした。

また日本やイギリスと違って、修士課程というものが基本的になく、5年一貫の博士コースしか博士号を取得する道はない。修士号をどこの大学で取っていようが関係なく、問答無用で5年一貫コースに入る。たぶん2年目に進級試験みたいなのを受けさせられ、その出来が悪いとキックアウトされると思うどこの学部も。最初の2年は必修みたいなものなので多分パスは出来ないのだろうけれど、あとの3年は論文さえかければ1年でも2年でもいいのだと思う。稀に超天才が3年とかで博士号を取得するケースを耳にする。以上の背景から、修士号=ドロップアウトなので、基本的に修士号だけ取得は残念賞扱いと話を聞いた。最近は修士課程も色々な学部で出来ているようだけれど、まあビジネスだと思う。MBA?言うまでもないですね。

また印象的だったのは兎に角"多様性"を入試の段階で重視することで、具体的には出願時に、人種、性別、信仰する宗教、親の年収、自分が思む自分の性別、軍歴の有無、性の対象となる性別等等、僕からすればどうでもいい事を兎に角聞かれた。特に州立大学だと、その州の出身者を優先的にとならなくてはいけないとか色々あるので、はっきりいってフェアな競争では全くない。例えばその年の定員が20人だったら、10人はアメリカ人で、半分は必ず女性にして、アジア系は3人までで、ヒスパニック系は何人で等など、自身のバックグラウンドに沿った枠の中で競争することになる。僕みたいな平凡な学生にとってはこれは結構死活問題だった。ワンちゃん改宗もありかもしれない。

またイギリスと違って、合格の裁量権が各先生にあるわけではなく、入試委員会が基本的にその年の学生の合否を決めている。まあ凄い先生に推薦書で押してもらえれば別に関係ないとは思うけれど、イギリスよりはあまりコネは効かないのかなと思うし、委員会の構成によって合否の結果も変わりそうなので、運に左右される部分も大きいと思う。己の運命力を信じて、願っていけ。

お金は基本的に合格=学費免除+生活費支給なので、イギリス程金策をする必要性は無いけれど、逆に言えば自前の奨学金を持って出願すれば学部が払うコストが減るので、合格する確率は高くなると思う。実際に最初は不合格だったけれど、日本での奨学金がその後もらえることになって、そのことを伝えたら合格になったケースもしっている。合否も交渉の余地が多分あって、僕も最初は修士までのコースしか合格しなかったけれど、他の大学で色々合格したことを伝えたら、結局博士コースで合格をくれた大学があった。行かなかったけれど、そんなことは僕の知ったことではない。



最後に話をまとめると、イギリスだろうがアメリカだろうが博士課程に合格したければ、いい成績をとって、各種試験でいいスコアをとって、いい推薦書を書いて貰う、これだけだと思う。とくに運命力と大事なのかな。願っていけ。

2017年11月19日日曜日

逐次モンテカルロ法(パーティクル・フィルター)について

僕は別に逐次モンテカルロ法(SMC)の熱心の信者というわけではないけれど,少し知識があるので覚書を書いてみる.簡単化のためにデータyが与えられた時の状態変数xのフィルタリング問題を例とするが,適用範囲はこの限りではない.

1 カルマンフィルターと比較して,SMCは非線形かつ非ガウスの一般的な隠れマルコフモデルを扱える.

2 MCMCとSMCとの違いについて.ターゲットが事後分布のときに,SMCは事後分布と共に周辺尤度も同時に計算できる.またこれらは各点で得られる.推定も名前の通り逐次的に行うので,バッチ推定ではなくいわゆるオンライン推定を行っている.計算コストもSMCの方が一般的に少ない.最終時点での分布のみに興味がある場合は各点での結果を保存する必要はない.バーンインの確認とか,ステップサイズはどうするのだとか,採択率をどの程度が最適だとかとか,そういったものも必要ない.あと分散のバウンド等もMCMCより比較的容易でかつシャープなものが得られていると思う.

3 SMCの基本的なアルゴリズムはインポータンスサンプリングとリサンプリングを逐次的に行うだけ.この組み合わせで色々できる.サンプリングを行ってリサンプリングを行うアルゴリズムより,リサンプリングを行ってサンプリングを行うアルゴリズムの方が推定精度はよい.数学的にはファイマン・カック方程式のある意味での近似として捉えられる.この解釈によってとにかく色々便利な数学的ツールが統一的にもたらされる.ちょっと具体的には,隠れマルコフモデルのフィルタリングの場合,ファイマン・カック・パス測度が事後分布にあたり,ポテンシャル関数がxを所与としたyの密度となる.

4 よく見かける誤解として,ブートストラップ・フィルターをSMCそのものとして紹介しているのを見かけるけれど,ブートストラップ・フィルターはSMCの一種ってだけであって,インポータンス・デンシティとして状態変数xが従う密度(マルコフ・カーネル)を用いたのがブートストラップ・フィルターとなる.この時ウェイトがxを所与としたyの密度となる.ウェイトは単なるラドン・ニコディウム微分なので,測度に関する絶対連続の仮定を満たせばインポータンス・デンシティとして別に何を使っても良い.

5 リサンプリングを行わない逐次インポータンスサンプリングは,サンプルサイズが増加していく毎に推定値の分散が指数的に発散していく.いわゆる劣マルチンゲール列になることが示すことができる.なのでサンプルサイズが大きい時に得られた事後分布の近似精度はかなりよくないことが多い.リサンプリングを行うとこの問題は消える.

5 しかしリサンプリングを行うといわゆる余計なモンテカルロ・エラーを加えるので基本的に推定精度は悪化する.なので各点で毎回リサンプリングをするのではなく,ある閾値(effective sample sizeを用いる)を越えたら行うようにするのが一般的である.またリサンプリングの結果として高い尤度を持つ同じ粒子がなんども選ばれて,事後分布の近似精度が悪くなる,Path Degeneracyと呼ばれる問題が不可避的に発生する.なのでパス全体に依存するようなSMCはサンプルサイズnが大きい場合に使い物にならなくなるので,そのようなアルゴリズムはなるべく避ける必要がある,粒子の数Nを各点で無限に増やせば数学的には必ず動くアルゴリズムだけれど,実際にはこのPath Degeneracyという問題はとにかく厄介.なので結局は原因は違えども,逐次インポータンスサンプリングと同じような問題にある程度は苦しむことがある.

6 リサンプリングを行うことによってシステムをリセットするような効果が生まれる.まず初期分布の影響は指数的に消える,また各点nで事後分布と周辺尤度の推定誤差をサンプルサイズとは独立に非漸近的にバウンドすることがきでる.このバウンドは基本的に次元が大きくなると緩くなる.これらをforgetting propertiesと呼び,リサンプリングを行う強い動機となっている.リサンプリングは長期での推定精度を保証するために,短期的に推定精度はある程度犠牲にする,いわば支払うコストみたいなものである.

7 SMCで近似された周辺尤度や事後分布は不偏推定量かつ一致推定量で,漸近正規性もみたすが,nに対する漸近理論とNに対する漸近理論を分けて考える必要性がある.基本的にはNに対しての漸近理論が多い.

8 感覚的には,逐次インポータンスサンプリングでは確実に起こる現象を,リサンプリングを加えることによって,SMCではアルゴリズムの工夫次第では何とかできるようになった位なもんだと思う.

9 SMCでパラメタの推定を行うのは難しい.SMCが周辺尤度の普遍推定量を与えることを利用してMCMCを組み合わせて行ったり,SMCをネストして用いて推定したり,Robbins–Monro型のアルゴリズムを利用した方法等が考案されている.

10 高次元での振る舞いは未だよくわかってないことが多い.ブートストラップ・フィルターの場合はある程度のレートで各点での粒子を増やせばある意味で安定することは示されているけれど,まあ自明な結果だと思う.

11 MCMCよりは容易に並列計算が行えるので,計算速度はある程度はやい.例えば LibBiというソフトでは,前述したSMCを用いた状態空間モデルのパラメタ推定を結構高速に行ってくれる. 

12 SMCに限らず,モンテカルロ法は積分を近似する手法であって,別にフィルタリングやスムージング等だけに応用されるものではない.例えばSMCを用いてナヴィエ・ストークス方程式を解く,みたいな応用研究も存在する.

2017年11月18日土曜日

どうでもいい予測の意味の違い

(統計的)機械学習で使われるモデルはブラックボックスなので、予測は上手くできるかもしれないけれど、なぜそうなったかは説明できないし、例えば将来何かルールが変わったときに、人々がどう行動を変えるのかを説明できないという批判を主に経済学のとある分野に携わっている人がしているのをよく見かける。結論から言うと、who caresじゃないのかな。

程度に差はあるけれど、まずアルゴリズム自体にはある程度数学的な基礎付けがある。MCMCやLASSOみたいに上手くいく理由がきちんと数学的に説明できるものもあるし、古典的な変分ベイズみたいにKL最小化が結局何を意味してるんだとか、他にもディープホニャララみたいに実証的に上手く働くのは分かるけれど、何故なのかはイマイチわからないとかそういう部分は多分にある。まあでもそれらの数理的な側面はこれから研究されて明らかになっていけばいい重要なトピックだと個人的には思うし、それに異を唱える人もそういないと思う。

なのでここでいうブラックボックスとはデータ生成過程そのもについての理論がないと仮定して話をする。

まず一部の経済学者が、何故を自分たちは説明できると思っているようだけれど、それは外からみればデータ生成過程が何らかの経済学の理論の結果として得られるものだと仮定してるからだと映る。その理論が正しい保障がそもそもないし、そういう理論の妥当性の検証をするのがエコノメの役割だったと思うのだけれど、何故か正しいことは所与として解析を行う構造推定だとかなんだとかいうのが流行っているらしい。経済学のコミュニティ内でそういう意味付けが求められるというのは分かるけれど、データ科学のコミュニティではそんなに求められてはいないと思う。コミュニティに近いエコノメ理論家ですら彼らが何をやってるのかよく分からない人が多いんじゃない。

勿論データ科学でもデータ解析をするさいに使う手法を決めるのだから、暗にある程度はデータ生成過程に仮定を置く。iidだとか隠れマルコフだとか用いるプライアーだとか。そのうえでどのモデルが良くフィットするのだとか、チューニングパラメターをどうすればいいのだとか、どのアルゴリズムが上手く働くだとか、そういった議論は客観的、実証的な側面から文脈をかえて活発に行われているし、一つの大きな研究トピックだと思う。

構造推定の裏にある理論の妥当性を信じない人達から、何でそのモデルが妥当だと言えるのか、何で様々なモデルをフィットさせてみないのかと聞かれたら、多分この分野の人達は外の人を納得させるような議論は出来ないと思う。アウトプットが上手く今あるデータを説明出来てるから正しいんだと言うならば、本末転倒だし。これは別に批判しているわけではなく、殆どの人はそんなことは気にしてないのではというだけである。

将来何かしらのルールが変わったら、その時にまた推定、学習させて予測すればいいだけで、何で現在において起こるかどうかも分からないことを想定した、反実仮想実験みたいなことをする必要があるのか、大体の人には分からないと思う。第一、例えば将棋で二歩が出来るようになったとしたら、棋士ですら将来自分がどう戦略を変える分からないと思う。フェアな例えではないかもしれないけれど、将来サッカーでキーパー以外も手を使えるようになったらどう戦略がかわるかなんて、考える意味すら希薄だし、誰が分かるのだろう。

想定外のことがおきたら、プログラムも予測できなくなるけれど、人もできなくなると思う。なのでそれが俺達は出来ると言える意味が良く分からないし、だからそれが強みだと言われても、違うコミュニティでは受け入れられ辛いと思う。過去のことについて、あの時AではなくBをしていたらCではなくDになっていただろうとか、その程度のことは強い仮定の下で言えると思うけれど。

僕が知っている限り、今のデータ科学で求められてるのは、なるべく少ない計算コストで精度よくデータなりモデルなりを予測・近似する手法で、特に実際にデータ解析をしている人程そういう思考が強い。MCMCでさえ正しいのは分かるけれど高次元だと計算が遅いから、と言われて現場では毛嫌いされている。無限っていつだ?と真顔で聞かれる。なのでサンプルサイズやパラメータの次元が1000とかその程度で、一回の推定に一週間だとか時間がかかって、かつ予測・近似精度が格別にいいわけでもない手法なんて、はっきりいって誰も気にしない。勿論使える状況もあるだろうけれど、いわゆるビッグデータの解析には全く向かない。良く分からない理論を所与として得られるデータの意味付けなんかも、やはり誰も気にしていないと思う。


最後に話はそれるけれど、アルファ碁は動学構造モデルの計量経済学における最高の成功例であるという主張をみた。同じ基盤技術を使ってるというだけで、多分野での研究成果を、自分が関わってる分野の研究の成果みたいに表現して、本当に下劣だと感じた。ディープマインドとかで同じ発言を是非してみてほしい。

自分たちの研究分野や研究成果を誇りに思うのは重要だと思うけれど、それが他のコミュニティでも当然受け入れられるだろうと思うのは自惚れだと感じる。