バナナでもわかる話

開設当初は計量経済学・統計学が専門の大学院生でした。今はデータを扱うお仕事をしています。統計学・経済学・投資理論・マーケティング等々に関する勉強・解説ブログ。ときどき趣味も。極力数式は使わずイメージで説明出来るよう心掛けていますが、時々暴走します。

階層ベイズ生存解析を用いたwebサイトの訪問者分析に関するStanでの実装

今回も少し難しめの話。
前回生存解析についてまとめたので、生存解析を使っているマーケティングモデルでも勉強してみようかと探してみたところ、次の論文があったので、論文モデルを整理しながらStanで実装してみようと思います。

前回の生存解析の記事
www.bananarian.net

今回まとめる論文
山口景子 『頻度の時間変化を考慮した階層ベイズモデルによるウェブサイト訪問行動の分析』
https://www.jstage.jst.go.jp/article/marketingscience/22/1/22_220102/_pdf/-char/ja

Stanについてよく知らない人はこの本がオススメ


目次

スポンサーリンク



webマーケティング研究の現状

webマーケティング研究

これは論文の内容ではありませんが、マーケティングサイエンス分野におけるwebマーケティング研究の現状について概観しておこうと思います。

まず、私が調べた限りですが、webマーケティングに関するマーケティング研究の関心は*1最近非常に高まっています。

しかし、その一方でwebマーケティング分析において欠かせないユーザー(or 顧客 or 消費者)のアクセスログデータは基本的に、各企業が管理しており、マーケティング研究者の手元にはありません。その為研究者がこれらの研究顧客分析を行うにあたっては、企業と共同研究を行う必要があり、ハードルが高く、報告も少ないのが現状です。

アクセスログの経営学的意義

経営学においては、消費者の購買行動(購買前、購買後含む)や購買心理を集計的な観点、個別的な観点から分析し、その結果を企業経営に役立つ形でインプリケーションとして提示することが求められます。

しかし購買前行動については、実際に使用できるデータ(消費者が商品購買前に何をしていたかに関するデータ)がほとんどなく、現象として観察することは出来ても計量が難しく、実際に役立てるのが難しいと言われてきました。

ただ今はというと、webサービスやECサイトの登場と共に、膨大なユーザーアクセスログが蓄積され、購買前の行動を十分に観察出来るようになりました。そのため、今経営学分野ではアクセスログを使った研究をしたいというニーズが高まりつつあるわけです。

one to one マーケティング

また、細かなユーザー情報を収拾できるようになったために、ウェブマーケティング分野ではone to oneマーケティングと言って、「消費者一人一人のニーズや状況を汲み取ったマーケティングコミュニケーションを行おう」という状況に移行しつつあります。

ウェブサイトの訪問頻度のモデル化

PV数

ウェブサイトのPV数はそのサイトの広告価値に影響します。そこでPV数を高めたり、今後の潜在的なリスクを見つけるためにどういった要因がPV数に影響を及ぼしているか分析出来たら嬉しいということになります。

さてここで、山口さんはとりあえずPV数を「訪問頻度」と「訪問一回当たりのPV数」に分解し、今回はそのうちの「訪問頻度」について考えてみたというわけですね。

モデル

尤度

 T日間において、ユーザー i=1,\cdots,Nがそれぞれ J_i回サイトに訪れたという状況を考えます。ユーザーiがj回目にサイトに訪れた日数を t_{(i,j)}と置くと、各回の訪問間隔 y_{(i,j)}=t_{(i,j+1)}-t_{(i,j)}はパラメータ \lambda_{ij}の指数分布に従うと仮定します。

ここで、 j=J_iの場合の y_{(i,j)}が本来は観測されないことに注意してください。観察期間 Tとの差を考え、指数分布の切断を考えます。


まず、 j≠J_iの場合の密度関数
 f(y_{(i,j)};\lambda_{ij})=\lambda_{ij}exp(-\lambda_{ij} y_{(i,j)})

 j=J_iの場合は次のように切断を考えます。

 y_{(i,j)}=T-t_{(i,j)}
この時、観測されていない変数 y_{(i,j)}^*=t_{(i,j+1)}-t_{(i,j)}は指数分布に従うため、
 f(y_{(i,j)}^*;\lambda_{ij})=\lambda_{ij}exp(-\lambda_{ij} y_{(i,j)}^*)

ここで、少なくとも y_{(i,j)}^*≧y_{(i,j)}であることはわかっているので、

 Prob(y_{(i,j)}^*≧y_{(i,j)})=1-Prob(y_{(i,j)}^*<y_{(i,j)})=1-F(y_{(i,j)}) = exp(-\lambda_{ij} y_{(i,j)})



更に、顧客の情報を織り込みたいので、次のようにパラメータを規定します。
1を含む顧客情報に関する(K+1)次元ベクトル x_{ij}を用いて、 \lambda_{ij}を次のように考えます。

 \lambda_{ij} = exp(\beta_i x_{ij})

ただし \beta_iは(K+1)次元パラメータベクトル。

事前分布と階層パラメータ

更に事前分布として次のように設定する。

 \beta_i ~ MultiNormal(\delta,\Sigma_{\beta})

 \delta~  MultiNormal(適宜設定)
 \Sigma_{\beta} ~  InvWishart(適宜設定)

共変量の選択

この論文では、目的上経時情報が知りたいので、そのような共変量を選択し、セミパラメトリックなcox比例ハザードモデルを行い、妥当な共変量を判断した後に、再度実行を行っている。

cox比例ハザードモデルについてはこちらをどうぞ
www.bananarian.net

Stanで実装してみた

R で擬似アクセスログを作ってみた - あらびき日記
access_logs.R · GitHub

アクセスログやユーザー情報のデータなんて当然私は持ち合わせていないので、上記サイトさんを参考にアクセスログデータもどきを作らせていただきました。

アクセスログ作成プロセスはこの記事では省略しますので、上記サイトさんを参考にしてください。

model="

functions{
	real TruncatedExp_log(real Y,real lam){
		real pr;
		pr = exp(-lam*Y);
		return(log(pr));
	}
}

data{
	//個人ID数
	int N;
	//特性数
	int CNUM;
	//各IDに対応する訪問回数
	int T[N];
	//T[N]の和
	int SUMT;


	//個人の特性
	matrix[SUMT,CNUM] CX;
	//Truncate情報
	int<lower=0,upper=1> Tru[SUMT];
	//DiffTime=Y
	vector<lower=0>[SUMT] Y;
}

parameters{
	//係数パラメータ
	vector[N] beta0;
	matrix[N,CNUM] beta;
	
	//事前分布のパラメータ
	vector[CNUM+1] BETA;
	cov_matrix[CNUM+1] SIGMA;
}

transformed parameters{
	vector<lower=0>[N] lambda0;
	vector<lower=0>[SUMT] lambda;
	vector[CNUM+1] Ibeta[N];

	for(i in 1:N){
		lambda0[i] = exp(beta0[i]);
		Ibeta[i][1]=beta0[i];
		for(c in 1:CNUM){
			Ibeta[i][c+1]=beta[i,c];
		}
		if(i==1){
			for(t in 1:T[i]){
				lambda[t]=lambda0[i]*exp(dot_product(beta[i],CX[t]));
			}
		}
		else{
			for(t in 1:T[i]){
				lambda[T[i-1]+t]=lambda0[i]*exp(dot_product(beta[i],CX[T[i-1]+t]));
			}
		}
	}
}

model{
	//尤度関数
	for(i in 1:SUMT){
		if(Tru[i]==0){
			Y[i] ~exponential(lambda[i]);
		}
		else{
			Y[i] ~TruncatedExp(lambda[i]);
		}
	}			
	//事前分布
		for(i in 1:N){
			Ibeta[i] ~ multi_normal(BETA,SIGMA);
		}

	//階層事前分布
	//適宜調整(仮に下のような感じにしてみた。下のよりももっと無情報に近づけるのが良い)
	BETA ~ multi_normal(rep_vector(0,CNUM+1),diag_matrix(rep_vector(10,CNUM+1)));
	inverse(SIGMA) ~ wishart(10,diag_matrix(rep_vector(10,CNUM+1)));
}

"

library(rstan)
fit=stan(model_code=model,data=data,seed=1234,par=c("Ibeta"))

雑感

かなり解釈しやすいモデルで、使える場面が多々ありそうだな~と思いました。
論文にもあるように、各人のパラメータを使ってセグメンテーションを行うと直感的にもわかりやすいかもしれないです。

*1:理由は後述

消費者の情報探索行動の目的を反映した購買モデルをStanで実装する

今回も中々高度な話をやっていきます。

消費者は何か物を買うときに「情報探索行動」を行います。
しかし、情報探索行動は特に物を買う気の無い消費者も行います。
例えば、特に買いたいものはないけど、流行をチェックするために服屋をぶらつくだとか、本屋の好きなジャンルのコーナーを見るのを日課にしているといったことがそれです。消費者がどちらの状態で情報探索を行っているのかによって、店側はその消費者に取るべきアプローチが変わってくるはずです。

目次

スポンサーリンク



参考文献

今回参考にした論文は次の通り
・Bloch et al .1986 Consumer search: An extended framework
・Newman and Stelin .1972 Prepurchase information seeking for new cars and major household appliances
・Doni 2017 情報探索の目的を考慮した購買決定モデル

主にDoni 2017の内容の要約とStanでの実装になります。
試しで使用するデータセットは持ち合わせてないので、とりあえず型だけ合っている架空のデータ(乱数)を使用して実行できるかどうか試しています。


ちなみにStanについてよく知らないという方は、この本がおすすめです。

情報探索行動

まず、情報探索行動について整理しておきます。情報探索行動は主に「自身の記憶」に対してアプローチする内的情報探索と、「外部」に対してアプローチする外的情報探索に分けられます。今回まとめるDoni 2017では外的情報探索のことを情報探索と呼び、モデル化を行っています。

購買前探索

ニーズを自覚している消費者(何かしら購買目的のある消費者)が行う外的情報探索を特に「購買前探索」と呼びます。

購買前探索は、購買に伴う不確定要素を削減し、集めた情報をもとに購買選択肢の評価や決定を行うために行います。

例えば、〇〇日に東京観光へ行くから、ホテルを予約したい。
予算の範囲内で極力快適なホテルに泊まりたい。
だから色々な情報サイトや口コミサイト、公式の情報を集め、不確定要素を削減し、どのホテルを予約するか決めようといったものがそれです。

進行的探索

ニーズを自覚していない消費者が行う外的情報探索を特に「進行的探索」と呼びます。

これは、製品知識の獲得や、その他欲求の解消が目的で、衝動買いなどが発生したり、今後の購買のための情報集めであることもありますが、基本的には購買目的で行っている行動ではありません。

モデル

情報探索の目的(購買前探索 or 進行的探索)が不明であるというのを隠れマルコフモデルに落とし込むことで定式化を行います。また、ある時点の消費者の購買魅力度(効用)の確定項はそれ以前の具体的な情報探索行動をもって定式化し、購買魅力度のパラメータに対して階層構造を与え、この魅力度のパラメータは消費者の属性が影響していると考えます。
以下詳しく説明していきます。

データ

~与えられている情報~
消費者 i=1,2,\cdots,I
消費者が店(又はECサイト)に訪れた機会 t=1,2,\cdots,T_i

消費者iが各機会tにおいて購買をしたか否か
 I_{it} =1⇒商品を購入
 I_{it} = 0⇒商品未購入

各消費者iの属性行列 Z_i
ex) 性別などなど

各消費者iの各機会tにおける情報探索行動に関する行列 X_{it}
ex)アクセスログなどなど


~与えられていない未知の情報~
消費者iが機会t時点でどのような目的 c_{it}で情報探索行動を行っていたか

消費者の購買行動のモデル

まず、消費者iは購買魅力度 U_{it}^{c_{it}}を持つと仮定します。
ここで、消費者が商品の購買を決める( I_{it} =1)とは次のような状況であると考えます。

消費者iと機会tに依存するある閾値 \tau_{it}が存在し、
 U_{it}^{c_{it}} >\tau_{it}^{c_{it}} ⇒I_{it} =1
 U_{it}^{c_{it}} ≦\tau_{it}^{c_{it}} ⇒I_{it} =0

更に、 X_{it}のうち、購買前探索に関係しそうな変数を X_{it}^{(1)}、進行的探索に関係しそうな変数を X_{it}^{(2)}と名付けます。この変数分類は、過去の購買前探索・進行的探索に関する学術論文を参考に決定します。

例えば次のような情報です。
・購買前探索を行っている消費者は同一商品の閲覧回数が増える
・進行敵探索を行っている消費者は多様な商品ページを見る
などなど。これらの情報はこれからのマーケティング関連の研究でどんどん増えていくであろう情報なので、必要な時に情報をアップデートしておく必要があります。


そして、購買魅力度を次のように規定します。
 U_{it}^{c_{it}} = \beta_i^{(c_{it})} X_{it}^{(c_{it})} + \epsilon_{it}^{(c_{it})}

更に購買閾値についてですがDoni 2017では
(Moe and Fader 2004,Johnson et al. 2003, Brucks 1985)
を参考に次のように特定化しています。

来店機会 tまでの累積購買回数 N_{it}を用いて、
 \tau_{it}^{c_{it}} = INIT_{it}^{c_{it}} +\beta_{0i}^{c_{it}} NPUR_{it} +E_{it}^{c_{it}}
ただし NPUR_{it} =1-exp(-N_{it})
 E_{it}は誤差項

情報探索行動目的のモデル化

マルコフ過程なので、遷移確率を考えます。
 Pr(c_{it}=s) =u_{it}(s)
 c_{it}が1次マルコフ過程に従うとすると、直前の状態にのみ依存するため、次のようになる。

 u_{it}(s) = \sum_{r=1}^2 Pr(c_{it}=s|c_{it-1}=r)u_{it-1}(r)

ここで遷移確率行列 \Gamma_{i}を次のように表す。

 \Gamma_{i} =\begin{pmatrix}
\lambda_i &&1-\lambda_i\\
1-\rho_i && \rho_i
\end{pmatrix}


初期状態の確率を \delta_i =(\delta_{i1},\delta_{i2})とおく。

選択確率

そして、 Pr(I_{it}=1|c_{it}=s) = \frac{exp(V_{it}^{s}-\tau_{it}^{s})}{1+exp(V_{it}^{s}-\tau_{it}^{s})}
ただし、 Vとは先ほどの購買魅力度の確定要素の部分とします。これはロジットですね。


こうして、長くなりましたが顧客iの購買有無のベクトルに関する尤度 l(\theta)を求めることが出来ました。
 l(\theta) =sum\{ \delta_i diag(Pr(I_{i1}|c_{i1}=1),Pr(I_{i1}|c_{i1}=2)) \Gamma_i \cdots \Gamma_i diag(Pr(I_{iT_i}|c_{iT_i}=1) Pr(I_{iT_i}|c_{iT_i}=2)) \}

階層構造

最後にパラメータ \betaに関して Z_iを使って階層構造を規定しておきます。

 \beta_i^{(s)} = MultiNormal(\Delta_s Z_i , \Sigma_s)


非常に雑なシミュレーションデータの作成

今回のメインはStanコードなので、形だけデータをそろえておきます。

I=20;MAXT=100
T=sample(1:MAXT,I,replace=TRUE)
Y=matrix(,I,MAXT)
for(i in 1:I){
	Y[i,1:T[which(1:I==i)]] <- sample(c(0,1),T[which(1:I==i)],replace=TRUE)
	Y[i,T[which(1:I==i)]:MAXT] <- 0
}
C=2
s1=5
X1=matrix(,I*MAXT,s1)

for(s in 1:s1){
index1=1
index2=MAXT
	for(i in 1:I){
		X1[index1:index2,s] <- rnorm(MAXT,mean=0.05*I,sd=sqrt(s))
		index1 = index1 + MAXT
		index2 = index2 + MAXT
	}
}

s2=4
X2=matrix(,I*MAXT,s2)

for(s in 1:s2){
index1=1
index2=MAXT
	for(i in 1:I){
		X2[index1:index2,s] <- rnorm(MAXT,mean=0.3*I,sd=sqrt(s))
		index1 = index1 + MAXT
		index2 = index2 + MAXT
	}
}

N=matrix(,I,MAXT)
for(i in 1:I){
	for(t in 1:MAXT){
		x=rpois(1,lambda=1)
		if(t>1){
			N[i,t] = N[i,t-1] + x
		}
		else{
			N[i,t] = x
		}
	}
}

zs=3
Z=matrix(,I,zs)
for(ss in 1:zs){
	for(i in 1:I){
		Z[i,ss]=rpois(1,lambda=sqrt(sqrt(i))*ss)
	}
}

data=list(I=I,MAXT=MAXT,T=T,Y=Y,C=C,s1=s1,X1=X1,s2=s2,X2=X2,N=N,zs=zs,Z=Z)

ちょっと工夫としては、Stanでunbalanced dataを扱うのは厄介なので、一旦 T_iを一番大きいところにそろえているところです。
Stanのモデルブロックで調整します。

Stanコード

では、今までの話をコードにしてみます。
事前分布は基本的にはDoni 2017を参考にしましたが、MCMCの初期値探索が難しかったため、かなり元の論文よりも制約をかけています。
一応回るようにはなってますが、この辺の制約は実データも交えて相談ということになると思います。

#unbalanced-dataの欠損部分に一時的に0を入れて、Stan内で0の部分は実行をかけないよう調整する
#調整をかけるのは購買の有無と説明変数行列
model ="
functions{
	//ある人iの購買の有無(T=1,...T[i])、T[i]、二種類の状況のVtベクトル配列、推移確率パラメータ、初期値を入れると、対数尤度を表す関数
	real lmm_log(vector Y,int T, vector Vt1,vector Vt2, real lambda, real rho, row_vector init){
		matrix[2,2] tr_mat;
		row_vector[2] prr;
		real pr;
		vector[2] vec1;
		vector[2] vec2;
		matrix[2,2] pr_mat;
		tr_mat[1,1] = lambda;
		tr_mat[1,2] = 1-lambda;
		tr_mat[2,1] = 1-rho;
		tr_mat[2,2] = rho;
		prr = init;
		for(t in 1:(T-1)){
			vec1[1] = 1-inv_logit(-Vt1[t]);
			vec1[2] = 1-inv_logit(-Vt2[t]);
			vec2[1] = inv_logit(-Vt1[t]);
			vec2[2] = inv_logit(-Vt2[t]);
			pr_mat = Y[t]*diag_matrix(vec1)+(1-Y[t])*diag_matrix(vec2);
			prr = prr*pr_mat*tr_mat;
		}
		vec1[1] = 1-inv_logit(-Vt1[T]);
		vec1[2] = 1-inv_logit(-Vt2[T]);
		vec2[1] = inv_logit(-Vt1[T]);
		vec2[2] = inv_logit(-Vt2[T]);
		pr_mat = Y[T]*diag_matrix(vec1)+(1-Y[T])*diag_matrix(vec2);
		pr = prr*pr_mat*rep_vector(1.0,2);
		return(log(pr));
	}

//  matrix seiho_Kronecker_product(matrix A,matrix B){
//  int Ad1;
//  int Bd1;
//  matrix[Ad1*Bd1,Ad1*Bd1] res_mat; 
 //   Ad1=num_elements(A[1]);
 //   Bd1=num_elements(B[1]);
 //   for(i in 1:Ad1){
 //     for(ii in 1:Ad1){
 //       res_mat[((i-1)*Bd1+1):i*Bd1,((ii-1)*Bd1+1):ii*Bd1]=A[i,ii]*B;
  //    }
  //  }
 //   return(res_mat);
//  }

//  vector vec_operator(matrix A){
//    int length;
//    int Ad1;
//    int Ad2;
//    vector[length] res_vec;
//    Ad1=num_elements(A'[1]);
//    Ad2=num_elements(A[1]);
//    length = num_elements(A);
//    for(i in 1:Ad2){
//      res_vec[((i-1)*Ad1+1):(i*Ad1)]=A'[i]';
//    }
//    return(res_vec);
//  }

}

data{
	//顧客数I,TのMAX値,実際のT,購買有無Y
	int I;
	int MAXT;
	int T[I];
	matrix[I,MAXT] Y;

	//状態数C=2
	int C;
	//状態1の説明変数s1種類
	int s1;
	matrix[I*MAXT,s1] X1;
	
	//状態2の説明変数s2種類
	int s2;
	matrix[I*MAXT,s2] X2;

	//閾値の為の変数
	int N[I,MAXT];

	//ベータに階層性を持たせるための顧客情報Z,zsは顧客情報の種類,
	int zs;
	matrix[I,zs] Z;
}

transformed data{
	real NPUR[I,MAXT];
	//NPUR = 1-exp(-N)
	for(i in 1:I){
		for(t in 1:MAXT){
			NPUR[i,t] = 1-exp(-N[i,t]);
		}
	}
}

parameters{
	matrix[I,s1] beta1;
	matrix[I,s2] beta2;
	matrix[I,C] beta0;
	//遷移確率
	real<lower=0,upper=1> lambda[I];
	real<lower=0,upper=1> rho[I];
	//初期値
	vector<lower=0,upper=1>[I] delta1;
	real INIT[I,C];
	//ベータとINITの階層パラメータ
	matrix[zs,s1+2] b1;
	matrix[zs,s2+2] b2;
	vector[s1+2] b01;
	vector[s2+2] b02;
	cov_matrix[s1+2] Sigma1;
	cov_matrix[s2+2] Sigma2;
	
}

transformed parameters{
	//効用と閾値
	matrix[I,MAXT] tau[C];
	matrix[I,MAXT] U[C];
	matrix[I,MAXT] Vt[C];
	matrix[I,2] D;
	//ベータの階層パラメータ
	matrix[I,s1+2] B1;
	matrix[I,s2+2] B2;
  //ベータの階層パラメータの連結
  matrix[zs+1,s1+2] bv1;
  matrix[zs+1,s2+2] bv2;
	//パラメータのまとめ
	vector[s1+2] theta1[I];
	vector[s2+2] theta2[I];

  ////ベータの階層パラメータの連結を行う
  bv1[1]=b01';
  bv2[1]=b02';
  for(ss in 2:(zs+1)){
    bv1[ss]=b1[ss-1];
    bv2[ss]=b2[ss-1];
  }  

  //U,tau,Vt,Dを作る
	for(i in 1:I){
		for(t in 1:MAXT){
			U[1][i,t] = dot_product(beta1[i],X1[(i-1)*MAXT + t]); 
			U[2][i,t] = dot_product(beta2[i],X2[(i-1)*MAXT + t]); 
			for(c in 1:C){
				tau[c][i,t] = INIT[i,c] - beta0[i,c]*NPUR[i,t];
				Vt[c][i,t] = U[c][i,t]-tau[c][i,t];
			}
		}
   	D[i,1] = delta1[i];
	  D[i,2] = 1-delta1[i];   
	}

  //Bを作る
	for(i in 1:I){
		for(ss in 1:(s1+2)){
			B1[i,ss] = b01[ss] + dot_product(b1[1:zs,ss],Z[i]);
		}
		for(SS in 1:(s2+2)){
			B2[i,SS] = b02[SS] + dot_product(b2[1:zs,SS],Z[i]);
		}
	}

  //thetaを作る
	for(i in 1:I){
		theta1[i][1] = INIT[i,1];
		theta1[i][2] = beta0[i,1];
		theta2[i][1] = INIT[i,2];
		theta2[i][2] = beta0[i,2];
    for(ss in 1:s1){
		  theta1[i][ss+2] = beta1[i,ss];
    }
    for(ss in 1:s2){
  		theta2[i][ss+2] = beta2[i,ss];
    }
	}

}	

model{
  //vec_operator(bv1) ~ multi_normal(rep_vector(0,(zs+1)*(s1+2)),seiho_Kronecker_product(Sigma1,diag_matrix(rep_vector(100,zs+1));
  //vec_operator(bv2) ~ multi_normal(rep_vector(0,(zs+1)*(s2+2)),seiho_Kronecker_product(Sigma2,diag_matrix(rep_vector(100,zs+1));
	for(i in 1:(zs+1)){
    for(ii in 1:(s1+2)){
      bv1[i,ii] ~ student_t(4,0,3);
    }
    for(ii in 1:(s2+2)){
      bv2[i,ii] ~ student_t(4,0,3);
    }
	}
  inverse(Sigma1) ~ wishart(10,diag_matrix(rep_vector(1,s1+2)));
	inverse(Sigma2) ~ wishart(10,diag_matrix(rep_vector(1,s2+2)));
	for(i in 1:I){
		D[i,1] ~ beta(0.001,0.001);
		lambda[i] ~ beta(0.001,0.001);
		rho[i] ~ beta(0.001,0.001);
		theta1[i] ~ multi_normal(B1[i],Sigma1);
		theta2[i] ~ multi_normal(B2[i],Sigma2);
	}
  for(i in 1:I){
    Y[i,1:T[i]]' ~ lmm(T[i],Vt[1][i]',Vt[2][i]',lambda[i],rho[i],D[i]);
  }
}"
library(rstan)
fit=stan(model_code=model,data=data,seed=1234,par=c("beta0","beta1","beta2"))

出力結果に関しては、そもそもデータがあれなので置いておくとして、こんな感じですね。
恐らく間違ってはいないはず。
ログデータを取れるweb屋さんとかだと使いようによっては便利かもしれないですね。

多項プロビットモデルの最尤推定をRで行う方法と数理、気を付けるべき点

今回は多項プロビットモデルです。少々高度な内容になります。

目次

スポンサーリンク



多項プロビットモデル

具体例

皆さんは離散選択モデルというモデルをご存知でしょうか。

主に人間の意思決定をモデル化する際に利用されるモデルです。
例えば、コンビニに行って、何を購入するか決めるとします。
あくまで例ですが仮に「何か甘いものが食べたい」という漠然とした目的をもってコンビニに入った人だと考えましょう。

ここでこの人は普通「甘いものが売られている売り場」へ向かうはずです。
そして、パッケージ、値段、ラインナップ等々を吟味したうえでたくさんある選択肢の中から一つの商品を選び、レジへと持っていったとします。

この選択の過程を今回は多項プロビットモデルというモデルを使って解析します。

数式

さっきの話を簡単にまとめると要は

パッケージ、値段等の情報を適当な重みを付けて比較し、「最もその人にとって良いと思える商品を1つ手に取った」という状況ですね。

つまり、その人にとっての何らかの選択基準が存在し、その選択基準に基づいて最も良いと思える商品を手に取ったと考えられます。

この時の個人 iの商品 jに対する選択基準 Y_{i,j}^{*}を以後効用と呼ぶことにします。
この効用とは経済学の概念で、例えばAさんが商品1を買うか商品2を買うか迷っているとした場合、Aさんはその商品に関する情報を使って効用 Y_{A,1}^{*},Y_{A,2}^{*}を頭の中で作り上げます。

そしてもし、 Y_{A,1}^{*}>Y_{A,2}^{*}なら商品1を購入し、 Y_{A,1}^{*}>Y_{A,2}^{*}ならば商品2を購入するといった状況を考えるわけです。


ここで、「いやいや、僕はそんなよく分からない効用なんてものを考えて行動していないし、そんなモデル意味が無いよ」という批判は一旦置いておきます。*1


この効用 Y_{i,j}^{*}は、誤差項 \epsilon_{i,j}とその人が持ちうる情報 I_{i,j}によって次のような線形関係があると考えます。

 Y_{i,j}^{*} = V_{i,j} + \epsilon_{i,j}
ただし、 \epsilon_{i,j}は期待値 0、分散共分散行列 \Sigmaの多変量正規分布に従うと仮定します。

更に、 V_{i,j} = \beta_{0}^{(i)} + \beta_{1}^{(i)} X_{(i,j),1} +\beta_{2}^{(i)} X_{(i,j),2} +,\cdots ,+\beta_{k}^{(i)} X_{(i,j),k}= X^{(i,j)}\beta^{(i)}
と仮定します。ここで Xとは値段や商品配置、パッケージ等の購買に関する判断基準を指します。
 \beta^{(i)}は個人 iが暗黙的もしくは無意識的に設けているパラメータベクトルで、 X^{(i,j)}は行列です。

当然添え字のつけ方で、各人共通のパラメータを設けたりすることも可能になります。
よって、 V_{i,j}を条件付けた際の Y_{i}^{*} の分布 P(Y_{i}^{*} | V_{i,j},j=1,2,\cdots , k)は次のようになることが分かりますね。


 P(Y_{i}^{*} | V_{i,j},j=1,2,\cdots , k) = \frac{1}{(2 \pi)^\frac{k}{2}\sqrt{ |\Sigma| }}exp(-\frac{1}{2}(Y_{i}^{*}-V_{i})' \Sigma^{-1} (Y_{i}^{*}-V_{i}))

ただし、 Y_{i}^{*},V_{i}はそれぞれ効用と判断基準を並べた j次列ベクトルとします。

与えられているデータ

ここで問題になってくるのは

効用 Y_{i}^{*}の値なんて観測できなくない??ってことです。
 iさんにとってのチ〇ルチョコの効用は10,チョコパイの効用は20だから、 iさんはチョコパイを購入したんだな~なんてことはわからなくて、実データでは、その辺の話を全部すっ飛ばして

 i さんは「チョコパイを購入した」という情報だけ得られます。

もし、効用まで観測できれば P(Y_{i}^{*} | V_{i,j},j=1,2,\cdots , k) に対して最尤法を適用しておしまいですが、観測出来ていないので簡単にはいきません。そこで、少し工夫する形で別の尤度を考えます。


今回の iさんは商品 Jを購入したとしましょう。

購買情報を踏まえた尤度

 iさんが商品 Jを購入したということは、次の状況が生じているはずです。


 Y_{i,J}^{*} ≧ Y_{i,j}^{*} for  j = 1,2,...k
こう書き換えても問題ありません。
 Y_{i,J}^{*} = max\{Y_{i,j}^{*} |j=1,2,...k \}


更に、 iさんの購買情報 \omega_{i}を次のように考えます。
 \omega_{i} = (\omega_{i,1},\cdots,\omega_{i,k})
ここで、 \omega_{i,J} = 1 \omega_{i,t} = 0
ただし、 t≠J

つまり、購入した商品には 1、購入しなかった商品には 0を与えたダミー変数ベクトルを用意するわけです。


そうすると、購買情報 \omega_{i}における分布 Prob(\omega_{i} | V_{i,j},j=1,2,\cdots , k)は次のようになるはずですね。

 Prob(\omega_{i} | V_{i,j},j=1,2,\cdots , k)
 = \int_{\{Y_{i,1}^{*}<Y_{i,J}^{*}\}},\cdots \int_{\{Y_{i,J}^{*}=-\infty\}}^{\{Y_{i,J}^{*}=+\infty\}},\int_{\{Y_{i,J+1}^{*}<Y_{i,J}^{*}\}}\cdots, \int_{\{Y_{i,k}^{*}<Y_{i,J}^{*}\}} P(Y_{i}^{*} | V_{i,j},j=1,2,\cdots , k) dY_{i,1}^{*},\cdots,dY_{i,k}^{*}

何かえらいこっちゃなってますが、要は多変量正規分布の密度を切断している分布になります。

まあ、当然ですがこの Prob(\omega_{i} | V_{i,j},j=1,2,\cdots , k)はオープンフォームになっているため、普通の方法で最尤推定量を求めることは出来ません。

そのため、プロビットモデルは長く日の目を拝むことなく、より計算が簡単なロジットモデルが発展しました。

しかし、今日ではコンピュータの計算速度の向上により、このプロビットモデルの最尤推定値を探索的に求めることが可能になっています。

GHKアルゴリズム

最尤推定値を探索的に求める方法は様々ありますが、例えばフィッシャーのスコア法では対数尤度の二回微分に関する期待値を利用します。
しかし、肝心の尤度がゴツすぎるため、微分したものの期待値なんてもってのほか。求めることが出来ません。
そこで、一旦乱数を取り出してやり、期待値の近似値を求めます。その近似値を利用して、探索解を求めるという手順を取ります。

その際、乱数の発生させ方に一工夫必要で、その取り出し方をGHKアルゴリズムと呼びます。

もう積分記号なんか見たくないとは思いますが、再度先ほどの尤度を書いてみます。

 Prob(\omega_{i} | V_{i,j},j=1,2,\cdots , k)
 = \int_{\{Y_{i,1}^{*}<Y_{i,J}^{*}\}},\cdots \int_{\{Y_{i,J}^{*}=-\infty\}}^{\{Y_{i,J}^{*}=+\infty\}},\int_{\{Y_{i,J+1}^{*}<Y_{i,J}^{*}\}}\cdots, \int_{\{Y_{i,k}^{*}<Y_{i,J}^{*}\}} P(Y_{i}^{*} | V_{i,j},j=1,2,\cdots , k) dY_{i,1}^{*},\cdots,dY_{i,k}^{*}


ここで、
 Y_{i,J}^{*} ≧ Y_{i,j}^{*} for  j = 1,2,...k
 Y_{i,j}^{*} = V_{i,j} + \epsilon_{i,j}

のため、
 V_{i,J} + \epsilon_{i,J} ≧ V_{i,j} + \epsilon_{i,j}
更に、
 V_{i,J} -V_{i,j} ≧   \epsilon_{i,j} -  \epsilon_{i,J}

ここで、 \epsilon_{i,j} -  \epsilon_{i,J} =E_{i,j^*} for  j^* = 1,2,...k-1
 E[E_{i,j^*}] =0, V[E_{i,j^*}] =\Sigma^*であるとします。

多変量正規分布には加法性があるので、 E_{i,j^*}の組も多変量正規分布に従います。
つまり E_{i}は多変量正規分布に従います。
さてここで、分散共分散行列 \Sigma^*に対してコレスキー分解を施し、
 \Sigma^{*} = \Gamma \Gamma^{'}と置くと、

 E_{i} =\Gamma \xi
ただし、 \xiは各要素がi.i.dに標準正規分布に従う確率変数ベクトル。

 \Gammaは上三角行列なので次のように書くことが出来ますね。

 E_{i,1} =\Gamma_{1,1} \xi_1
 E_{i,2} =\Gamma_{1,2} \xi_1+ \Gamma_{2,2} \xi_2
 E_{i,2} =\Gamma_{1,3} \xi_1+ \Gamma_{2,3} \xi_2+ \Gamma_{3,3} \xi_3
 \vdots

つまり、長くなりましたが何が言いたいかというと、多変量正規分布に従う確率変数は、期待値ベクトル(今回は0)と分散共分散行列さえわかっていれば、普通の標準正規分布から得られるわけです。


こうして、対象の多変量正規分布からサンプリングを行い、更に必要な期待値計算をモンテカルロ近似によって計算し、それを用いて適当な最適化アルゴリズムを使うというのがGHKアルゴリズムの概要です。

このアルゴリズム、上の説明を見ても分かる通り
ある一つの選択肢を基準に他の

Rによる多項プロビット回帰

「GHKアルゴリズム、、なんて大変なんだ」と思った方もご安心。Rではmlogitパッケージを使って簡単に推定することが出来ます。

使用するデータ
data("TravelMode", package = "AER")
> head(TravelMode)
  individual  mode choice wait vcost travel gcost income size
1          1   air     no   69    59    100    70     35    1
2          1 train     no   34    31    372    71     35    1
3          1   bus     no   35    25    417    70     35    1
4          1   car    yes    0    10    180    30     35    1
5          2   air     no   64    58     68    68     30    2
6          2 train     no   44    31    354    84     30    2

AERパッケージのTravelModeデータです。
mlogitパッケージのExampleで使っていたので、そのまま使うことにします。
https://cran.r-project.org/web/packages/mlogit/mlogit.pdf

ちなみに各変数は
individual : 個人iの識別番号(i=1,...200)
mode : 乗り物の種類(選択肢)
choice : その乗り物を選んだか否か
wait : 乗り場における待ち時間(当然自家用車なら0)
vcost : 乗り物を使うのにかかるコスト
travel : かかる時間
income : 個人の収入
http://127.0.0.1:23970/library/AER/html/TravelMode.html

今回のデータはchoice毎に行があり、individualが重複する形の縦に長いデータとなっています。
このようなデータフレームをlong型と呼びます。ちなみに1行で個人データをひとまとめにしているデータをwide型と呼び、例えば次のようなものがあります。

> head(Fishing)
     mode price.beach price.pier price.boat price.charter catch.beach catch.pier catch.boat catch.charter   income
1 charter     157.930    157.930    157.930       182.930      0.0678     0.0503     0.2601        0.5391 7083.332
2 charter      15.114     15.114     10.534        34.534      0.1049     0.0451     0.1574        0.4671 1250.000
3    boat     161.874    161.874     24.334        59.334      0.5333     0.4522     0.2413        1.0266 3750.000
4    pier      15.134     15.134     55.930        84.930      0.0678     0.0789     0.1643        0.5391 2083.333
5    boat     106.930    106.930     41.514        71.014      0.0678     0.0503     0.1082        0.3240 4583.332
6 charter     192.474    192.474     28.934        63.934      0.5333     0.4522     0.1665        0.3975 4583.332
mlogit.data関数

mlogitを使用する前に、解析で使うデータフレームをmlogitに適用しやすい形に変換しておくことが無難です。
その際の変換に用いるのがmlogit.data関数です。

ただし、別にmlogitを使用する際に関数内で適切な指定をしてやれば変換していないデータフレームに実行をかけることも出来ます。

TM <- mlogit.data(TravelMode, choice = "choice", shape = "long",
alt.levels = c("air", "train", "bus", "car"),id.var="individual",)

choiceにはlong型の場合は選択肢を選んだか否か(yes or no ; 1 or 0)に関する変数を指定します。wide型の場合は選んだ選択肢の入った変数を指定します。
shapeは"long" か "wide" か使用するデータの型を選択します。
alt.levelsで選択肢を改めて指定しておく必要があります。long型の場合は必須ですが、wide型の場合は必要ありません。
id.varでlong型の場合に個人識別をしておく必要があります。

mlogit関数

mlogit関数で実行します。デフォルトはただのlogitなのでprobitで行う場合はprobit=TRUEとします。

pr <- mlogit(choice ~ wait + vcost | income |travel , TM,
probit = TRUE)

formulaを|で区切っていますが、これはそれぞれ次のような意味になっています。
個人iの選択jに関する効用 U_{i,j}について
 U_{i,j} = \alpha_{(0,j)} + \beta_w Wait_{(i,j)} +\beta_v Vcost_{(i,j)} + \gamma_{j} Income_{i} + \omega_{j} Travel_{(i,j)}

添え字に注目してください!

実行&結果
> summary(pr)

Call:
mlogit(formula = choice ~ wait + vcost | income | travel, data = TM, 
    probit = TRUE)

Frequencies of alternatives:
    air   train     bus     car 
0.27619 0.30000 0.14286 0.28095 

bfgs method
22 iterations, 0h:0m:11s 
g'(-H)^-1g = 6.11E-07 
gradient close to zero 

Coefficients :
                     Estimate  Std. Error z-value  Pr(>|z|)    
train:(intercept)  0.11400922  0.50036804  0.2279 0.8197623    
bus:(intercept)   -0.25750600  0.46140096 -0.5581 0.5767789    
car:(intercept)   -1.70532857  0.67252768 -2.5357 0.0112223 *  
wait              -0.02591688  0.00596662 -4.3436 1.401e-05 ***
vcost             -0.00571812  0.00332269 -1.7209 0.0852629 .  
train:income      -0.02256556  0.00677919 -3.3287 0.0008727 ***
bus:income        -0.01125839  0.00862617 -1.3051 0.1918437    
car:income        -0.00637486  0.00788242 -0.8087 0.4186627    
air:travel        -0.01473926  0.00360134 -4.0927 4.264e-05 ***
train:travel      -0.00291735  0.00079889 -3.6517 0.0002605 ***
bus:travel        -0.00313288  0.00089814 -3.4882 0.0004863 ***
car:travel        -0.00296051  0.00089036 -3.3251 0.0008840 ***
train.bus          0.94852133  0.38768124  2.4467 0.0144190 *  
train.car          1.06025359  0.37869969  2.7997 0.0051147 ** 
bus.bus            0.38464898  0.16055087  2.3958 0.0165838 *  
bus.car            0.27405472  0.37816033  0.7247 0.4686330    
car.car            0.57842389  0.24356722  2.3748 0.0175584 *  
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Log-Likelihood: -172.56
McFadden R^2:  0.39187 
Likelihood ratio test : chisq = 222.39 (p.value = < 2.22e-16)

地味に"gradient close to zero" が重要です。これが"last step couldn't find higher value"となっている場合は問題有です。
coefficientsの最後の方に並んでいる(選択肢).(選択肢)はそれぞれの分散共分散です。一個分散が欠けているのはスケールを揃えるのに用いているためだと思われます。
McFadden R^2は一応経験的な適合度指標として使われています。あまり紋切型で考えるのもイマイチですが、0.5くらいが目安なので今回は少し小さめ。

各個人の選択確率予測値はprobabilitiesに入っています。

> head(pr$probabilities)
           air       train         bus       car
[1,] 0.1463981 0.230907011 0.069925608 0.6030067
[2,] 0.3550429 0.194959842 0.012583377 0.4531929
[3,] 0.4270799 0.025918716 0.052778939 0.4559133
[4,] 0.4102070 0.007009033 0.004554783 0.5633459
[5,] 0.2017710 0.479339802 0.114605173 0.1839117
[6,] 0.1113763 0.478200196 0.122358769 0.3097929

多項ロジットとの違いは?(IIA特性)

何故よくわからんGHKアルゴリズムを使ってまで(計算に時間をかけてまで)probitを使うのだ。多項ロジットで良いではないかと思う人もいると思うので違いを説明しておくと

多項ロジットの場合は、誤差項の共分散が無いので、各選択肢の代替関係が同等であるという仮定が入っています。
例えば車か飛行機かを選ぶにあたって、その選択肢の中に追加でバスが入ってきたとします。この場合、例えば陸路で安価に行けた方が嬉しいなら車の選択肢を選んでいた人がバスを選ぶようになる確率よりも飛行機の選択肢からバスに移る確率の方が多いように感じますし、旅行先が渋滞まみれのような場所では車で目的地に向かうのは不便なので、飛行機よりも車を選んでいた人の方がバスに流れていくはずです。このような変化に多項ロジットは対応出来ません。

これをIIA特性(無関係な選択肢からの独立)と呼びます。

これを解決するためにはネストロジットという手法もありますが、多項プロビットでも解決することが出来ます。

IIA仮定が妥当か否かの検定(ハウスマン・マクファーデン検定)

IIA仮定が妥当であれば多項ロジットを使う方が計算の点からも良く、またロジットは選択肢が実際には観測されていない場合でも機能するため、多項ロジットだって無下には出来ません。

そこで仮説検定を行います。
Hausman and McFadden (1984)
EconPapers: Specification Tests for the Multinomial Logit Model
で提案されたハウスマン・マクファーデン検定では

帰無仮説: IIA特性を持つ
対立仮説 : IIA特性は持たない

という検定を行います。検定においてはフルモデルでの推定結果と、選択肢を一個取り除いた場合の推定結果を比較します。
もしIIA特性があれば、直接関係のない選択肢が消えてもオッズ比は(誤差の範囲で)自然な値になるはずです。

mlogitパッケージではこの検定をhmftest関数として与えています。

> x <- mlogit(choice ~ wait + vcost | income |travel , TM,reflevel="car")
> g <- mlogit(choice ~ wait + vcost | income |travel , TM,reflevel="car",alt.subset=c("car","bus","train"))
> hmftest(x,g)

        Hausman-McFadden test

data:  TM
chisq = 27.572, df = 9, p-value = 0.001124
alternative hypothesis: IIA is rejected

> g2 <- mlogit(choice ~ wait + vcost | income |travel , TM,reflevel="car",alt.subset=c("car","bus","air"))
> g3 <- mlogit(choice ~ wait + vcost | income |travel , TM,reflevel="car",alt.subset=c("car","air","train"))
> hmftest(x,g2)

        Hausman-McFadden test

data:  TM
chisq = 178.24, df = 9, p-value < 2.2e-16
alternative hypothesis: IIA is rejected

> hmftest(x,g3)

        Hausman-McFadden test

data:  TM
chisq = 41.746, df = 9, p-value = 3.656e-06
alternative hypothesis: IIA is rejected

選択肢を除いたmlogitとフルモデルのmlogitで比較を行っています。
今回はどの場合でも棄却(p-value <0.05) されたため、IIA特性は妥当とは言えないということです。
よって、今回は多項プロビットモデルを利用しています(当然他のIIA仮定を緩めたロジットを利用しても良い)

離散選択におけるその他の選択肢

ここは、これから解説記事を書く中で少しずつリンクを貼ってこうと思います。
とりあえずはざっと羅列。

①多項プロビット
今回説明したもの

②多項ロジット回帰
IIA特性があるが、解釈のしやすい離散選択モデル
クローズドフォームなので計算が楽


②ネステッドロジット
IIAを緩和すべく、選択肢をネストしたロジット
選択の順序関係は考慮していないことに注意

③ネステッドロジットのIIAを更に緩めたモデル
paired combinatorial logit モデル
ordered generalized extreme value モデル
cross-nested logit モデル

④generalized nested logit
上のネストロジットを一般化したモデル
IIAはかなり緩和されるが、
パラメータの置き方や解釈の仕方が難解で計算量も多い

⑤C-logit モデルとpath size logit モデル
確定項の方で相関を表現したモデル

②~⑤までを総称してGEVモデルと呼ぶ。
GEVと多項プロビットの大きな違いは、確率項の大きさが異なっていても多項プロビットは利用できることで、GEVでは確率項の大きさが同一であることを仮定している。

⑥離散連続選択モデル
離散選択を行った後に連続値の選択を行う場合には離散連続選択モデルが利用される。これについてはmultiple discrete-continuous extreme value モデルとそれを一般化した generalized multiple discrete-continuous extreme value モデルがある。

このまとめについては下pdf参照
http://www.trans.civil.nagoya-u.ac.jp/~yamamoto/papers/TE2012.pdf


また、このpdfにはありませんが
⑦多変量プロビットモデル

というものもあります。
これについてもいつか記事にします。

*1:この記事を全部読めば、その手の批判は的外れであるということが分かるかと思います。ただ、この時点でそういう批判が出てしまう気持ちはよくわかります。

多項プロビットモデルだけで1冊書いてる本を見つけた

今日大学の図書館の書庫を彷徨っていたら、面白そうな本を見つけました。

Multinomial Probit: The Theory and Its Application to Demand Forecasting

なんとこの本、multinomial probit(多項プロビットモデル)だけで1冊の本になっているらしく、多項プロビット関連の話でそんなに書くことあるのだろうか....と思いつつ興味本位で借りてみました。

スポンサーリンク


一応目次だけ書いておくと
chapter1 An introduction to Disaggregate Demand Modeling in the Transportation Field
chapter2 Maximum-Likelihood Estimation: Computational Aspects
chapter3 Statistical Aspects of Multinomial Probit Model Calibration
chapter4 Prediction: Mechanical Aspects
chapter5 The Statistical Interpretation of Predictions
Appendix A~D

という構成。

とりあえず、一通り読んだらまた追加でレビューでも書こうと思います。

サービスの4つの特性と統計モデリングでの注意点

今日は、有名なサービスの4特性を通じて、サービスに関するデータ解析で注意すべき点を記事にしようと思います。

スポンサーリンク



サービス

まず、サービスとはという話ですが、例えばコンビニや銀行、ファミレス等での接客や企業に対するコンサルタント、美容院のカットなどを指します。

※経済経営学部に入るとまず初めに「サービス」と「製品」という商品の区別を教わりますが、製品とは私たちがすぐ想像するところの商品のことだと思ってもらって問題ないかと思います。例えばテレビやスマホ、冷蔵庫などといった形があるような商品のことですね。

サービスの性質

無形性

まず、サービスと製品の大きな違いに「無形性」があります。
その名の通りサービスは製品と違って物質的な性質が無い、簡単に言ってしまえばサービスには形がないということを指しています。

「そんなのあたりまえじゃないか!」と思うわけですが、これは意外と重要な性質で、

まず、形が無いので顧客が事前にその内容や質をチェックすることが出来なく(難しく)なります

例えばこれが製品であれば、掃除機を買いたいから実際の店舗で吸引力を比較するだとか、マッサージチェアの購入検討のために電気屋さんで30分ほど試しに使ってみるなんてことが出来るわけですし、ニ〇リでお皿を購入したいとなれば、実店舗で手にとってその質感や色み等々を確認することが出来ます。

一方これが例えばある企業がコンサルタント会社さんに新しい事業に関する戦略を立ててもらうことになったとして、そのコンサルタントの質は実際にそのサービスを購入し、消費しない限りわかりません。何故なら形がなく、事前のチェックが不可能だからです。

他にも例えば、福岡県に旅行に行こう!となってホテルを予約するとします。このホテルが実際どのようなサービスを提供してくれて、どの程度の質かは実際に泊まってみないとわかりませんよね。


この意味で、サービスと製品は明確に性質が異なるわけです。
故に、一般に顧客にとっては「製品を購入することよりも、サービスを購入することの方がリスキー」に感じるわけです。
これを、顧客の知覚リスクと呼びます。

そのためサービスを提供する企業は顧客の知覚リスクの低減を目的に、例えばホテル会社であれば自社のホテルの内装やサービス内容などをホームページ上に公開したり、権威ある賞や資格を保有していることを示したりするわけです。学習塾などであれば体験入塾制度なんてのもありますよね。

スポンサーリンク



同時性

次にサービスの同時性という性質について説明します。
同時性とは、生産と消費がほとんど同時に発生するという性質を指します。

製品でこれはまず起こりませんよね。工場で作られた冷蔵庫が作られた瞬間使用され始めるなんてことはないはずです。

一方、ファミレスの接客なんかを考えると、店員は実際にお客さんと関わることで初めて接客を行います。つまり生産を行うわけです。そしてそれと同時にお客さんは店員の接客を受けることで消費していますよね。

このような性質をサービスの同時性と呼んだりします。

異質性

製品は、工場で規格を管理しておけば、たまに不良品は出るかもしれませんが、ほとんど同じ質の商品が出来上がります。

一方でサービスはかなり対応するサービス提供者ありきのものになります。これをサービスの異質性と呼びます。
これを管理するのは中々難しいわけですが、極力同じになるよう事前研修を徹底したりするわけですね。

消滅性

普通製品は、在庫管理を行います。ちょっと多めに作っておいて注文が入る度に出荷するなんてことをよくやるわけですが、
サービスに関してはそんなことは出来ません。

生産と同時に消滅するのがサービスです。これを消滅性と呼んだりします。
これがサービス提供の難しい所で、サービスは生産(コストとして発生)されていても、消費が伴わない場合があり得ます。
具体的に言うと、客の入らない美容院だと、美容師の給料は時給単位で発生しているのに、お客はこないからコストだけが増えるとかそういう話です。

統計解析における注意点

無形性

顧客はサービスを購入する際、リスクを極力減らすために数々の情報を利用して、購買行動を決定します。
例えば初めて入る飲食店の口コミを事前に調べておくだとか、価格帯からそのサービスの質を推定するといった行動です。

そして、そうした情報を元に購買を決定するわけですが、この際顧客の購買行動をモデリングするにあたっては次のような問題点があります。

1. 使用経験とマーケティング変数
まず、初回の購買において、口コミ等の情報を集めたりするのはそうですが、例えばラーメン屋A、ラーメン屋Bのどちらを選ぶかといった場合に、ラーメン屋Aはネットの情報しか得られないが、ラーメン屋Bは一度実際に行ったことがあり、「おいしいかった」という経験的な事実があるとします。この場合、「今日は確実においしいラーメンを食べたい」と思えば、ネット上の情報は頼りにせず(というか、見ることもせず)、一度入ったことのあるラーメン屋Bに入るはずですし、今回は「比較のためにラーメン屋Aに入ってみるか」と思うかもしれません。

つまり、マーケティング変数は、顧客の初期購買の行動分析においては信用たり得ますが、それ以降の購買は使用経験が影響を及ぼしている可能性が高く、マーケティング変数だけでは説明出来ない結果になりやすいということになるわけです。


2. リスクの顧客依存
また、どの程度サービスの不透明さに対してリスクを感じるかは当然顧客に依存します。
全くリスクに感じない人であれば、何も調べず、何の情報もなくサービスを購買することもあるかもしれません。

「髪さえ切ってもらえればどこでもいいや」

なんて考えて適当に美容院を選んでいる男の人とかはそんな感じかもしれません。
こんな感じの顧客であれば、価格帯だけを見てあとは自分の家に近いかどうかでサービスを選ぶなんてこともしているかもしれませんよね。

スポンサーリンク



同時性

1. 顧客と提供者の相互関係
サービスは、サービス提供者と顧客が同じ場所に居合わせるという少々特殊な商品です。
そのため、当然サービス提供者と顧客の相互作用も発生していると考えられます。
その相互作用の結果が事後的なサービス品質の評価に影響し、リピート購買に繋がるかが変わって来ている可能性も否めません。

異質性

1. サービス提供者の質
サービスの異質性を考えれば、当然サービス品質はそのサービス提供者による影響が大きいと言わざるを得ません。
そのため、単純なモデリングで異質性を排してしまうと実態にそぐわない結果を導く恐れがあります。


このように、サービスは普通の製品とは大きく異なり、実際にサービス関連のモデリングをするにあたっても複雑なモデルを作ってやる必要が生じます(当然、得られているデータのサンプルサイズとの兼ね合いも考えた上で、妥協するところは妥協する。)

ミクロ経済学の勉強のための参考書紹介~初級レベルから上級レベルまで

寒いですね。いつの間にか年末になってしまいました。今年も早かった....

さて、年末といえば(?)、そろそろ大学院入試があるころでは無いでしょうかね!?

そこで今回は、ミクロ経済学を学ぶための参考書を次の3つのレベルに分けて紹介していこうと思います

☆初級レベル
経済学初学者の学部生が、無理なく勉強を進めることの出来る参考書

☆中級レベル
初級レベルは大丈夫だが、今一歩大学院入試を受けるレベルには届いていない人向けの参考書

☆上級レベル
経済系の大学院で勉強をしていくにあたって、本格的なミクロ経済学を学ぶための入門書

正直、上級レベルを突破すれば、後は自身の関心に合わせて個別のトピックについて深く調べたり、論文を読んだりした方が良いので、今回はこの3段階に分けて参考書を紹介したいと思います!

スポンサーリンク



初級レベル

・ミクロ経済学の第一歩

ミクロ経済学が何を扱っている学問なのか、イマイチピンと来ていない!いきなり難しい話をされてもついていけない!そんな人にオススメ。
具体例や絵がふんだんに入っていて、イメージを持ちながら読み進めることが出来る本です。


・マンキュー ミクロ経済学

定番のマンキューです。個人的にはあまり好きでは無いですが、とても丁寧に書かれているので、大学の入門授業でよく使われている良書です!


・今までで一番やさしいミクロ経済学

その名の通り、かなり易しく丁寧に話を進めてくれる本です。


・レヴィット ミクロ経済学基礎編

幅広く丁寧にミクロ経済学の概要を学べる本です!こちらもかなり分かりやすい良書です!

中級レベル

・レヴィット ミクロ経済学発展偏

先ほど紹介したミクロ経済学基礎編の続きです。少し理論的な話が出てきます。基礎編と合わせて読めば、かなりミクロ経済学に対する理解が深まります。


・ミクロ経済学の力

話題の神取ミクロです!難しいトピックスも多いですが、比較的わかりやすく取っつきやすい書かれ方がされていて、入門書としてとても良いです!


・ミクロ経済学

所謂林ミクロです。理論に関する議論が丁寧になされているので、特に理系の方にお勧め!当然文系でも読めます!


・ミクロ経済学

私が学部3年生~4年生の時の勉強で使った本です。所謂西村ミクロってやつですね。個人的にはオススメです!とても丁寧に理論の説明が書かれていて、かなり網羅的です。地味に注釈に重要なポイントが書かれてあったりする。


・入門ミクロ経済学

ヴァリアンのミクロ経済学ですね。こちらの本は私は大学院に入ってから読んだのですが、試験前に読んでおくと良いだろうな~と思ったのでのせました。

上級レベル

・Microeconomic theory

大学院レベルのミクロ経済学本として真っ先に出てくるのはこれですね。これを避けて経済学修士を終えることはまず無いと思います。それくらい有名な本です。実際内容もかなり網羅的、しかも各トピックスの説明も非常に丁寧。経済学修士1年目ではまず、この本の内容を1年かけて追うことが大事な勉強になります。


・ゲーム理論 新版

ミクロ経済学といえばゲーム理論ですよね。ゲーム理論のトピックをしっかり網羅的に知りたいのであれば、恐らくこの本がベストです!少々集合論関連の知識が無いと辛いですが、全部読み切ればかなり力がつく本になっています。



こんな感じでしょうか。このあたりの内容を一通り抑えて置けば、あとは個別トピックに勉強を移していけるのではないかな~と思います!


【その他参考書記事紹介】

計量経済学について学びたい方はこちらも参考にしてみてください!
www.bananarian.net

統計学についてはこちら!
www.bananarian.net

SQLについてはこちら!
www.bananarian.net

順序選択モデルの数理と解釈~ベイズモデルの構成法~

今日は順序選択モデルというモデルについて書くことにします。

スポンサーリンク



具体例

例えば焼肉屋さんで、ライスを頼みたいと思ったとします。
この時に小ライスか、中ライスか、大ライスが選べるとします。

これ、この順序選択モデルを知らない人だと真っ先に「多項選択モデルを当てはめてみよう」という思考になるはずです。

そんでもって小ライスを0,中ライスを1,大ライスを2として、その選択を考えることになります。

しかし、これは少し不自然で、例えばこれがカルビを選ぶかサガリを選ぶか、それともミノを選ぶかみたいな選択だと0,1,2と割り当てれば良いわけですが、今回って選択肢に順序がありますよね。 小ライス<中ライス<大ライス といった具合です。

この場合何が起こるかというと、量に順序があるということは、人間の選択に順序が影響し始めます
「今は小腹が空いている程度だから、大ライスは頼めないなー」なんて考えて、大ライス未満の順位にあるライス、つまり中ライスか小ライスを選択するだとか、「とてもお腹が空いているから小ライスや中ライスのような下の順位のライスではなく、大ライスを頼もう」なんて考えたりだとか。

そういった状況なのだから、モデルを使うにあたっても順序を考慮したモデルを考えた方がより状況を反映した結果が得られるはずです。


「ん?0,1,2を割り当てた時点で、数字上0が1番小さくて2が1番大きいのだから、順序は考慮できているのでは?」

と思う方もいらっしゃるかもしれませんが、これは少々まずくて、数字上の距離を仮定してしまっています。
どういうことかというと、要は小ライスと中ライスの距離、中ライスと大ライスの距離が1であり、小ライスと大ライスの距離が2であるという仮定のもと大小を決めていますよね。量的変数ならまだしも、質的変数でこの仮定はあまり説得力がありません。

数理モデル

それではモデル化していきます。
J+1個の順序付けられた選択、0,1,....Jがあるとします。

ここで、選択を行う主体iが、何かしらの情報(説明変数、特徴量)を元に、主観的に自身の効用を最大にするような値 s_iを決定すると考えます。

つまり、次のような状況です。
 s_i = \beta X_i +\epsilon_i

更に、iが取る選択を y_iとし、適当な境界を持って選択判断をしていると考え、その境界に関するパラメータ \gammaを考えます。

 s_i <\gamma_1ならば y_i = 0
 \gamma_1< s_i <\gamma_2ならば y_i = 1

 \gamma_J < s_iならば y_i = J

といった具合ですね。
ただし、 \gamma_1<…<\gamma_Jです。


この時、例えば y_i = 0を取る確率は次のように考えることが出来ます。
 P(y_i = 0 |\beta ,X_i)=P(s_i <\gamma_1 | \beta,X_i)=P(\epsilon_i <\gamma_1-\beta X_i | X_i,\beta)

同様に、 y_i =1であれば、
 P(y_i = 1|\beta,X_i)=P(\gamma_1-\beta X_i <\epsilon <\gamma_2 - \beta X_i|\beta,X_i )

と得ることが出来ますね。

後は尤度としてPに正規分布を割り当ててやればプロビットとなりますし、ロジット分布を割り当てればロジットとなります。
標準正規分布の分布関数を \Phiと書くことにすると


 P(y_i = 0 |\beta ,X_i)=\Phi (\gamma_1-\beta X_i)
 P(y_i = 1 |\beta ,X_i)=\Phi (\gamma_2-\beta X_i)-\Phi (\gamma_1-\beta X_i)

と書けますね。


尤度

よって、尤度 l(\beta,\gamma,y,x)は次のように表せます。

 l(\beta,\gamma,y,x) \propto \sum_{j=0}^J 1_{y_i =j}(\Phi (\gamma_{j+1} -\beta X_i )-\Phi(\gamma_j - \beta X_i) )


まあ、比例記号で書いているところからお察しの通り、制約等々中々複雑なのでベイズを使うのが無難でラクチンです。


事前分布

ベイズを使うということで、 \gammaの事前分布はどうしましょうという問題が発生します。
制約として \gamma_1<…<\gamma_Jがあるので、例えば次のように考えてやります。


 \gamma_{j+1}-\gamma_jが非負であればよいと考え、各 \gamma_{j+1}-\gamma_jに対して独立にガンマ分布を仮定する。



こんな感じですね。ベイズだとそんなに難しくなく実装できるかなと思います。また今度、試せるデータもあれば、Stanでの実装例も書きたいと思います。

相関関係から因果関係を判断するために注意すべきこと

今日は因果関係の話でも。

スポンサーリンク


「収入が増えれば消費が増える」とか「広告を打てば売上が上がる」とかいったことを考えるにあたって、簡単に回帰分析をしてみましょう~!みたいな話はよく聞くと思います。

で、中級者くらいになると、そこから一歩ステップアップして「回帰分析で見ているのは、相関であって因果じゃない」という話を聞き始めると思います。回帰分析に限らず、因子分析や構造方程式モデリングでも同じです。

しかし、そうはいっても「まあ、注意はするけど、、、相関関係は考えるのがラクチンだし、、じゃあ因果関係を厳密に示すにはどうすればいいかっていうのも思い浮かばないし、、、注意しつつも回帰分析すればいいってこと?」みたいな話で有耶無耶にしている方も多いかもしれません。


そこで今回は、相関関係から因果関係にジャンプするためには何に注意すればいいかを書いていきたいと思います。

前後関係

まず、これはよく聞く話ですが、最低限前後関係はわかっていなければいけません。
例えば「雨が降る」と「地面が濡れる」という因果関係があったとします。
これは雨が降ったことが原因で地面が濡れたという結果があるわけですが、地面が濡れたから雨が降るわけではないですよね。

このように、因果には常に前後関係が伴います。よって、前後関係がおかしいのであれば因果が成り立つわけないですよね。

他に原因が無いか確認する

例えば、「地面が濡れる」という結果は、近所のオバサンが水をまいたから濡れている可能性もあるわけです。朝雪が降って、解けて濡れているだけかもしれません。このように、結果だけ見れば他に原因がある可能性があり、一面的に原因っぽいものを見つけてくるだけでは因果関係を特定したと言うことは出来ません。

トレンド

因果関係を考えるにあたって、事前にトレンドを排除しておかなければなりません。
例えば、「広告を打った」から「売上が上がった」という因果関係を考える際に、もしかしたら、全体的な市場のトレンドとして、売上が上がる傾向があっただけで、広告を打ったことは直接関係が無いかもしれない。

こうした状況を考えるために、例えば同時期に広告を打った商品と似ている、広告を打っていない商品であったり、その競合を比較するということが考えられます。

ただ、その時セレクションバイアス(元々広告を打っていない企業は弱っている企業でお金が無かったなど)が無いかは考えなければいけません。

カウンターファクター

本当にその原因が無ければその結果は発生しないのかといったことも考えなければなりません。
雨が降ったことで地面が濡れるということを確認したいのであれば、晴れの日に地面が濡れていないことも確認しなければいけませんよね。



当然実際問題、因果関係を考えるたびにこれらの条件を全て調べ上げてとなると骨が折れるし、ビジネスの現場で毎度こんなことを考慮する必要はありませんが、知ってるかどうかは結構行う分析の質に繋がってくると思います。

また、こういうことを知っておくと、まあ少なくとも常識的に前後関係はしっかりあって、他にも影響はありそうだけど、そこまで大きく影響を及ぼすことはなさそうだからこういう分析をしました~といった説明をすることも出来ます。

数学・統計初心者から始めて計量経済学で中上級者になるためのオススメ本

前に書いた統計本紹介記事(下リンク参照)が、とても評判がよかったので、計量経済学も書いてみようかと思います。
www.bananarian.net

スポンサーリンク


目次

計量経済学とは?

多分経済学部以外の方だと「計量経済学って何?」って人もいるかと思います。そこで簡単に計量経済学の説明をしようと思います。余談ですが経済学部だと科目名が計量経済学ではなく、「エコノメトリックス」になってたりするので、「エコノメ」と略称で呼んでる人も居たりします。

計量経済学を語るには、簡単に経済学を説明しなければいけません。(あまり適当なことを書くと有識者の方々に怒られそうですが、)長々書いても読者が退場するので簡単に。

経済学は、基本的には人間の行動・営み、またそれに伴って生じた経済事象全般をモデルを通して理解・解釈し、場合によってはその後の変化や動き介入による効果を分析する学問です。

人間の行動・営みとは、例えばコンビニのチョコレートコーナーに行って、どのチョコレートを購入するか選ぶといったミクロな視点から、企業や政府が大規模なお金を動かすようなマクロの視点まで全て含みます。その後の変化・動き・介入とは例えば、じゃあチョコレートを値上げしたら今まで購入していた人がどれくらい購入を控えるのだろうかだとか、政府が税金を増やしたら経済効果はどの程度だろうかとかそういうものを指します。

そういったものをモデル(数式)から考えるわけですが、当然実際に出てくるデータはモデル通りにはなりません。
所詮人間の営みなので、分散が大きく、誤差も増えます。しかし、誤差はあれど、一定の法則性を見出すことが出来るのであれば、モデルとして利用する価値はありそうですよね。


そこで、経済学の中でも特に計量経済学という学問は主に次のような役割を持ちます。

① モデルの妥当性を評価するために計量方法を考える⇒それにあたり、モデルを確率的なモデルに拡張し、統計学を利用する
② ①とは逆のアプローチで、データから確率的な影響を考慮した上で法則性を見出し、それをモデルとして提案する


非研究者志望が計量経済学を学ぶメリット

で、それじゃあ計量経済学って経済学を研究する人しか使わないんじゃないの?勉強する意味あるの?と思われる方もいると思うので簡単に宣伝しておきます。

・文系対象の本が多い
結構数学わからないんでしょ?のスタンスで説明してくれる本があるので、数学を習熟していなくても計量経済学を学びながら適宜各種教科書を参照することで、勉強することが出来ます。ついでに数学・統計学力も増します。実際、私も計量経済学を通して、行列演算の有用性に気が付きました。

・経済学徒が取っつきやすく統計学を学べる、そして活かせる
経済学部出身者が、統計学を学ぼうとすると結構苦労するんです。急にパラメータ空間とか言い出して混乱するし、とっつきやすくてまともな本は洋書ばっかりだし、ゴリゴリ証明が入って、(初学者にとっては)結局何に役立つのかわかりにくくハードルも高かったりします。何より、いざ統計学を学んでも、で、経済モデルはどうやって計量すればいいの?といった場面で手が止まりがちです。

・普通の統計学的な知見+αの分析が出来るようになる
ここ、多少賛否あって、思考停止で計量経済学モデルが優れているんだと考えちゃう方には向かないかもしれませんが、例えばお仕事をするにあたって顧客データに対して統計分析を行うにしても、昨今の統計・機械学習ブームの結果、やれることが飽和しつつある状況にあります。そんな時に、ある程度妥当な仮定をもって計量経済学的な手法も提案できるようになれば、お仕事の幅が広がるかもしれません(そういったお仕事の現場では数学科・物理出身の方が多く、計量経済学的な考え方に反発する方も多数いらっしゃいます。彼らの発想も理解しつつ説得力のある提案が出来るとベストかと。)。

・純粋な統計学と比較してみるのも面白い
これはメリットというより面白さですが、統計学と計量経済学ですが、やってること、考え方が微妙に(大きく?)違います。両者の考え方を学ぶと、両者の良い所悪い所なんかが見えてきて、より両者の学びを深めることが出来たりもしますよ!


そんな計量経済学ですが、まあお察しの通り少々ハードルが高く、最低限、統計学・微積分・行列計算 が出来ないと何をやっているのかよくわかりません。そこで今回は、高校の数学ⅡBくらいまでなら学習したけど、それ以降の数学は全く知りませんという状態から、大学院修士レベルにたどり着くために必要な本と、勉強プロセスを紹介していこうと思います。
読者の方は自分のレベルに合わせて参考にしてみてください。


高校以降数学やった記憶がありませんという方

まずは、簡単に数学をやりましょう。本当に簡単な教科書で大丈夫です。

・穴埋め式 線形代数 らくらくワークブック

・やさしく学べる微分積分

・コアテキスト 統計学


この3冊、どれも数学よくわかりません、知りませんって人を対象にメチャメチャ簡単なところから説明を始めるので、かなり取っつきやすい3冊になっています。オススメです。ちなみにうちの大学(経済学部)の最初の数学の授業で使われているものです。


経済学をよく知らない人向け

計量経済学をやるのは良いけど、まず、経済学をよく知らないんだけどって人のために、メチャメチャ簡単な入門書を紹介しておきます。
経済学をガチでやる人はそれはそれで別に勉強すればいいので、とりあえずそうじゃない人向けにはこの辺りで問題ないかと。

・ミクロ経済学の力

・マンキュー経済学(マクロ編)

※当然マンキューはミクロ編もあります。
※読みやすさとしてはミクロ経済学の力の方が読みやすいかな~と思ったのでそっちを載せました。

スポンサーリンク


計量経済学の入門

はい、とりあえず計量経済学をやる前の準備としてはこれで完了です。
当然難易度の高い計量経済学を学ぶにあたっては、それ相応の事前勉強をしなければいけないのですが、とりあえず入門としてはこの辺で大丈夫かと。

・計量経済学の第一歩 実証分析のススメ

和書の入門書としては、私的には一番オススメです。計量経済学って何をやるの?どんなことを大事にしてるの?といったことを実践を通して学ぶことが出来ます。


・計量経済学

個人的にはあまり好きではないのですが、大学の教科書としてよく使われています。数式展開が丁寧に書かれていて、どういう計算をしているのかはとてもわかりやすく、文系向け!という気もするのですが、どういう目的で数式展開しているのかを読んでいて見失いがちになるので、個人的には△。人によっては合うかもしれません。


・Introductory Econometrics

個人的に、入門書として一番オススメなのはこれです。ただ、洋書なので英語がちょっと...って人はとっつきにくいかもしれません。
英語を読むのが苦にならないのであれば、この本が一番丁寧。しかもデータを取ってきて実践することも出来ます。かなり学べる一冊です。
ちょっとお高いのですが、これ1冊で何とかなると思えばそれほど高い買い物でもないかなとは思います。


ベイズ計量経済学入門

当然、ベイズ版の計量経済学もあります。
これについては断然次の本がお勧めです。

・ベイジアン計量経済学

割と初心者でもいけるかな?確率を使った議論が少々厄介なので、確率論の本を一冊事前にやっておけば難なく読めると思います。


中級の計量経済学本

先ほどの入門書が終われば、かなり数学力・計量経済学力共についてるかと思いますので、とっかかりとしては結構簡単に次の本が読めると思います。

・計量経済分析

後で貼るグリーン本の和書版です。元の洋書版だと基礎的な数学の説明を省いてますが、和書版だと数学わかんないでしょ?という想定で始めてくれるので、かなり取っつきやすいです。

ちなみに洋書だとこちら

全体的な読みやすさとしてはやはり、洋書の方が良いとは思います。日本語だと、結構説明を省いてる個所や、強引に訳した箇所があり、読みにくさがあります。



ちなみに、この本が経済学部の大学院の講義で使われる参考書になるので、この本で大体一般的な経済学部の大学院水準まで到達できます!
(当然、計量経済学専門の人としては足りない)


ここまでくればかなり計量経済学中級者と言えるかも!
後は各種トピックについてもう少し詳しく勉強するか~といった方向へ向かうので、もっとニッチに専門書をあさるなり、論文を読むなりすることが出来ます。

前回の統計学本紹介の時と比べるとかなり、短めに収まりましたね(私自身が、統計学の方で苦労した経験の方が多いため笑)


先ほど言ったことの繰り返しになりますが、計量経済学と、統計学を比べながら、両者にどのような考え方の違いがあるかを考えてみるのも勉強になって面白いですよ。

追記

オススメ本記事、他にも書いたので、お時間があればご覧ください!

ミクロ経済学のオススメ本
www.bananarian.net

統計学のオススメ本
www.bananarian.net

SQLのオススメ本
www.bananarian.net

人間の選択をモデル化するプロビット・ロジットモデルの違いと経済学的解釈法

今日は人間の選択をモデル化する方法について書いてみようと思います。


一応事前にプロビット・ロジットについてネット記事が無いか漁ってみたのですが、

・機械学習の文脈で説明されていることが多い
・プロビット・ロジットの選択基準が説明変数の観点から説明されていることが多い(実はもう一個ある)

ということで

「じゃあ私は経済学的な観点から定式化を行って、統計学的な観点からこのモデルの説明をして、モデルにおける2種類の選択基準に関して説明をすることで差別化して書いてみよう~」

なんて考えてみました。
それでは記事を書いていきます。


ちなみに、プロビット・ロジットや経済学の話はとりあえず良いから、この二つのモデルの違いだけ知りたい!という方は

↓下の方にあるこの見出しから読むと良い感じに理解できるかと。

プロビット・ロジットの特徴



スポンサーリンク


人間の選択とは

まず人間の選択について少し経済学の観点から考えます。
例えば今、夕食の献立を考えているとしてとりあえずザックリ「肉」か「魚」かみたいなものを選択したいとします。

人によっては「その日の気分」のみで決めるかもしれないし、とりあえずスーパーへ行ってみて「値段」を見て決めようと思うかもしれないし、はたまたスーパーで実際に見て「おいしそうな方」にしようと思うかもしれないし、今の基準を全て総合して決めるかもしれません。

とにかく人によって何かしらの選択の「基準」があって、その「基準」に従って

「今日は肉より魚だな」

とか決めているわけです。


この考え方は妥当でしょうか?

①「いやいや私は、どっちでも良いと思ってスーパーに赴いて先に目に入った方を取るから何も考えてないよ」
②「肉とか魚だから今の例が成り立つけど、例えば何か甘いものを食べたいなと思ってコンビニに行き、何も考えずにほとんどランダムに近い形でチョコレートを手に取り、買うことだってあるよ」

みたいな反応もあるかもしれません。

ただ、これは解釈の拡張次第でどうにでもなって
①は「一番初めに目に入ったものを買う」という基準に従って選択をした結果であり、
②は確かに何も考えて居なかったかもしれませんが、ということは最も手に取りやすい配置にあった商品が取られる可能性が高かったはずです。つまり無意識に配置によって選ばされていると考えれば「配置」という基準に従って選択をした結果です。


つまり、こんな風に人間の選択を考えると実際にその人が基準に従って選択をしているか否かに関わらず、柔軟に当てはめることが可能になります。


選択の定式化

さて、このままでは数式に出来そうにないのでもう少し考えてみます。
先ほどの選択で「肉より魚だな」と選ぶ時、人は何らかの基準に従って「肉を買った時の嬉しさよりも魚を買った時の嬉しさの方が大きいな」と思ったから魚を選んだと考えましょう。

"嬉しさ"と言ってしまうと多少語弊があるので"メリット"と言い換えてもらっても問題ありません。

この嬉しさ(メリット)は測ることは出来ませんが、各人が自分の頭の中で大きさを比べて、その大きい方を選んでいると考えます。
この嬉しさ(メリット)を経済学の用語で効用と呼びます。


肉を選ぶ際の効用を U_{肉},魚を選ぶ際の効用を U_{魚}とし、

「肉より魚を選んだ」というのを次のように表現します。

 U_{肉}< U_{魚}



そして、この効用の大きさ自体は分かりませんが、選んだ基準(価格や配置等々)を数値化することは可能です。

そこで回帰分析の要領で次のように考えてやろうというわけですね。



 U_{肉}=X' \beta_{肉} +\epsilon_{肉}
 U_{魚}=X' \beta_{魚} +\epsilon_{魚}

そんでもって U_{肉}< U_{魚}という式に放り込んでやれば

 X' \beta_{肉} +\epsilon_{肉}<X' \beta_{魚} +\epsilon_{魚}
つまり
 X' (\beta_{魚}-\beta_{肉}) +(\epsilon_{魚}-\epsilon_{肉})>0

よってこの式は
 X' \beta +\epsilon>0

と変形出来ました。この左側は回帰分析でよく見る形ですね。


ただ、このままでは回帰分析としては扱えません。だって Uの値はわからないわけですから。

そこで次のように考えます。

 U_{肉}< U_{魚} Y=1,
 U_{魚}< U_{肉} Y=0と表現することにします。

そして、 Xという情報(先ほどの価格等の基準)が与えられた時に Y=1を取る確率を Prob(Y=1|X)と書くことにすると

 Prob(Y=1|X)=Prob(X' \beta +\epsilon>0|X)


ですね。

ちなみにこれは次のように解釈することも出来ます。
 Y^*=X' \beta +\epsilonという線形関係を考えます。この時点では Y^*\in \mathcal{R}なので、

 Y^*≧0ならば Y=1
 Y^*<0ならば Y=0

つまり、実数範囲を取る潜在的な変数 Y^*が取る値に対して観測される変数 Y=0,1を結びつけるわけです。
統計や機械学習の文脈でこの話を学ぶ人はこちらの解釈の方が慣れてるかもしれませんね。

とにもかくにもこれでモデル化が出来ました。このプロセスが経済学的な観点から人間の選択を数式化するプロセスになります。



もっと背後のモデル

次に、個々の選択確率がわかったので、もっと背後にある確率を考えましょう。単純に考えると、0か1かという選択はベルヌーイ型の分布を考えるべきですよね。そこで次のような質量関数を考えます。

 f(y_i |X_i)=\pi_i^{y_i}(1-\pi_i)^{1-y_i}
 y_i=0,1

そしてこの \pi_iを先ほど考えた個々の選択確率に置き換えてやるわけです。
 \pi_i=Prob(Y_i=1|X_i)=Prob(X_i' \beta +\epsilon_i>0|X_i)


プロビット・ロジットモデル

ところで確率 Prob(Y=1|X)と確率 Prob(Y=0|X)はどうやって決めればよいでしょうか。

今のままだと漠然と Probとかいうよくわかんない記号でごまかしているにすぎませんよね。


そこでこの Probに対して、もう少し具体的な確率モデルを与えてやろうという方針が生まれ、そのモデルとしてプロビットモデル、ロジットモデルがあります。

それでは一つずつ見てみましょう。


プロビットモデル

別にロジットから説明しても良いんですけど、"プロビット"モデルの名付け親が日本人なので愛着をこめてこちらから説明します笑
先ほどの統計学的な解釈を思い出してみます。


 Y^*=X' \beta +\epsilonという線形関係を考えます。
この時点では Y^*\in \mathcal{R}なので、

 Y^*≧0ならば Y=1
 Y^*<0ならば Y=0

という話でした。線形回帰ではよく正規分布が考えられるので、プロビットでも正規分布を考えます。
使い勝手がよいので標準正規分布で考えましょう。つまり

 \epsilon ~Normal(0,1)

更に Y^*≧0とはつまり X' \beta +\epsilon≧0,

よって -X' \beta ≦\epsilonです。ということは

 -X' \beta ≦\epsilonならば Y=1
 -X' \beta >\epsilonならば Y=0

と書きなおすことが出来ます。
 \epsilonは標準正規分布に従うと仮定しているので図に表すとこんな感じ。

 -X' \beta ≦\epsilonならば Y=1
f:id:bananarian:20181115151037p:plain


この赤色部分を全部足してやれば Y=1の確率が出るよねということです。よって
 (標準正規分布の密度関数)=\phi(t)とおくことにすると

 Prob(Y=1|X)=\int_{-X'\beta}^{\infty}\phi(t)dt=1-\int_{-\infty}^{-X'\beta}\phi(t)dt

また、同様にして
 Prob(Y=0|X)=\int_{-\infty}^{-X'\beta}\phi(t)dt


これで分布による表現が出来ました。あとは先ほどのベルヌーイ型関数に放り込んで最尤法を行えば完了です。
ベイズの文脈で語るのであれば事前分布を設定して事後分布を見てやればおしまいです。

ロジットモデル

次、ロジットですね。こちらはどちらかというと数学的な扱いやすさから、よく使用されるモデルです。
 Prob(Y=1|X)=\frac{exp(X'\beta)}{1+exp(X'\beta)}

によって定義されます。これは何なんだというとロジスティック分布関数と呼ばれる分布です。少しマイナーですかね。


どういう点で数学的に扱いやすいかというと、オッズ比として解釈出来る点が強いです。
こちらのサイトさんのロジットの紹介がわかりやすかったので載せて置きます。是非参考にしてみてください。
ロジットとプロビットの使い分け - アドファイブ日記


※ちなみに、他にもモデルは存在します。記事下に補足として書いておくので興味のある方はご覧ください。


プロビット・ロジットの特徴

プロビットは標準正規分布を使っているので当たり前ですが、ロジットの方も左右対称になります。
というか、裾の重さを除いて標準正規分布とロジスティック分布はほぼ同じです。
実際-1.2から1.2の範囲においてこの二つの分布はほぼ同じ確率を弾きます。
ただ、ロジスティック関数の方が裾の重い分布となるのが特徴です。


「ほぼ同じならば、オッズ比を使えるロジットの方が扱いやすそうだし、ロジットの方がいいじゃ~ん。」

ということで考えるのが面倒な場面ではロジットを使う人が多かったりしますが、一応しっかり考えると次の二つの観点からモデル選択を行う必要があります。

①使っている基準(独立変数、Xのこと)の取りうる値の範囲が広いか否か(変数が多くある場合は、そのうち特に重要なものを見る)

②0,1の選択(従属変数、Yのこと)のデータでどちらか片方のデータがほとんど無いような場合



①はよく言われる話で、先ほど載せたサイトさんでも"サチる"といった言葉で説明していました。
要はロジットの方が裾が重いので、X'\betaがメチャメチャ大きい値を取った時の Y=1を取る確率、小さい値を取ったりした時の Y=0を取る確率がある程度大きいく、独立変数が広い範囲で動く場合はロジットを使うのが良いということです。
イメージとしては0から1へのシフトが緩やかって感じですね。


で、今回説明したかったのが②の方で、

例えば10000人のデータを取ってみて Y=0を選んでいる人が20人とかしかいなかったーなんてデータがあったとします。

この時、もしプロビットの方を適用しようとすると、標準正規分布の両端はほぼゼロに近いのでプロビットは使えません
これは別の記事で説明した「正規分布を適用するモデリングでは、外れ値(?)のようなデータに弱い」という話に似ていますね。

ただ、外れ値モデリングの話でもやりましたが、突飛な特徴がなく、0も1も万遍なくあるならより無難な分布を考慮できるという背景からプロビットを使った方が良さそうです。


プロビット・ロジットの使い分け

よって、次のような事が言えるかと思います。

まず、 Y=0, Y=1の内どちらかのデータが極端に少ない場合はロジットを考慮に入れ、
更に(重要な)説明変数の取りうる値の範囲が広いなら、ロジット使用が濃厚。

また、モデルを使用した後にオッズ比を考えたいような場面でもロジットが良い。

上の状況に該当しなければ、無難にプロビットの方が良い。


もし、他にも使い分けが考えられるよーという方がいれば是非コメントを頂けるとありがたいです。


オマケ

ちなみに、上の話では次のような疑問に答えていません。

対称性

ロジットもプロビットも左右対称な分布になっていました。この仮定は妥当なのか?という疑問があるかもしれません。そんな時はこのモデル

ガンベル分布を使ったガンベルモデルです。
 Prob(Y=1|X)=exp(-exp(-X'\beta))

ちなみに補対数対数モデルというのもあります。
 Prob(Y=1|X)=1-exp(exp(X'\beta))

独立性

プロビットもロジットも各選択が独立だという仮定を置いて、同時分布では各確率の積を取っているが、誤差の共分散は本当にゼロなのか?という疑問もあり得ます。そんな時は誤差を多変量に拡張した、多変量プロビットモデルというモデルがあります。

多項選択

今回は「肉」と「魚」の2種類の選択だったから Y=0, Y=1のモデルで良かったけど、サーロインステーキか焼肉か焼き鳥かケバブかどれにするみたいな選択だったらどうするんだという疑問も考えられますね。

そんな時はベルヌーイ型を多項分布に拡張し、多項プロビットモデル・多項ロジットモデルを考えたりします。