モデル
情報探索の目的(購買前探索 or 進行的探索)が不明であるというのを隠れマルコフモデルに落とし込むことで定式化を行います。また、ある時点の消費者の購買魅力度(効用)の確定項はそれ以前の具体的な情報探索行動をもって定式化し、購買魅力度のパラメータに対して階層構造を与え、この魅力度のパラメータは消費者の属性が影響していると考えます。
以下詳しく説明していきます。
データ
~与えられている情報~
消費者
消費者が店(又はECサイト)に訪れた機会
消費者iが各機会tにおいて購買をしたか否か
⇒商品を購入
⇒商品未購入
各消費者iの属性行列
ex) 性別などなど
各消費者iの各機会tにおける情報探索行動に関する行列
ex)アクセスログなどなど
~与えられていない未知の情報~
消費者iが機会t時点でどのような目的
で情報探索行動を行っていたか
情報探索行動目的のモデル化
マルコフ過程なので、遷移確率を考えます。
が1次マルコフ過程に従うとすると、直前の状態にのみ依存するため、次のようになる。
ここで遷移確率行列
を次のように表す。
初期状態の確率を
とおく。
選択確率
そして、
ただし、
とは先ほどの購買魅力度の確定要素の部分とします。これはロジットですね。
こうして、長くなりましたが顧客iの購買有無のベクトルに関する尤度
を求めることが出来ました。
階層構造
最後にパラメータ
に関してを使って階層構造を規定しておきます。
非常に雑なシミュレーションデータの作成
今回のメインは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を扱うのは厄介なので、一旦
を一番大きいところにそろえているところです。
Stanのモデルブロックで調整します。
Stanコード
では、今までの話をコードにしてみます。
事前分布は基本的にはDoni 2017を参考にしましたが、MCMCの初期値探索が難しかったため、かなり元の論文よりも制約をかけています。
一応回るようにはなってますが、この辺の制約は実データも交えて相談ということになると思います。
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屋さんとかだと使いようによっては便利かもしれないですね。
