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

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

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

参考文献
情報探索行動
- 購買前探索
- 進行的探索
モデル
非常に雑なシミュレーションデータの作成
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についてよく知らないという方は、この本がおすすめです。
<br />

情報探索行動

まず、情報探索行動について整理しておきます。情報探索行動は主に「自身の記憶」に対してアプローチする内的情報探索と、「外部」に対してアプローチする外的情報探索に分けられます。今回まとめる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屋さんとかだと使いようによっては便利かもしれないですね。