Edited at

時系列顧客ロイヤリティの算出

More than 5 years have passed since last update.


概要

標題のとおり顧客ロイヤリティを以下の書籍を参考に実装してみました。

といっても、過去に作ったモノを繋ぎ合わせているので多少おかしな部分があるかもしれません...

(私のローカル環境ではまともに動きました。)

おや?っと思われる部分がある場合はコメント頂けますと助かります。


コード

以下にコードを記載していきます。

多少分量があるので注意してください。

まずはmain側のファイルのコードです。


one_to_one_mrk_run.py

# coding: utf8


import time
import sys
import numpy.random as npr
from numpy.oldnumeric.linear_algebra import inverse
import scipy
from scipy import linalg
import random
import scipy.stats.distributions as dis
from scipy.stats import uniform,chi2,invgamma
import scipy.stats as ss
import numpy as np
import math
from scipy.linalg import inv,cholesky,det
import scipy.optimize as so
import one_to_one_mrk as hbm
from multiprocessing import *
from numpy import genfromtxt

# データインポートファイルパスの定義
inpt_D_filename = 'D_sliced.csv'
inpt_Z_filename = 'Z_sliced.csv'
inpt_payment_filename = 'pymData_sliced.csv'
inpt_C_filename = 'C_data_sliced.csv'
inpt_visit_filename = 'visit_sliced.csv'

### データinput-------------------
# Dデータの取得(共通性を使って個人毎のパラメータを推定するのに必要なデータ)
D = genfromtxt(inpt_D_filename, delimiter=',')
# Cデータの取得(日平均消費金額)
C = genfromtxt(inpt_C_filename, delimiter=',')
# PMデータの取得( INVを求めるための購買日別の金額 )
PM = genfromtxt(inpt_payment_filename, delimiter=',')
# yデータの取得( 訪問flagによって切断正規分布のa,bを求める )
y = genfromtxt(inpt_visit_filename, delimiter=',')
print "done input"
### ------------------------------

# ランダムシードの発生
random.seed(555)

##--定数の定義-------
TIMES = len(PM[0,:])
nhh = len(PM[:,0])
SEASONALYTY = 7 #季節変動の周期
RP = 50000 # MCMCでのサンプリング数
keep = 10 # 表示のデータ間隔
nz = 2 # 来店効用を表す式の説明変数の数(Z_iの数)
nD = len(D[0,:]) #Dのベクトル長
m = 1 + 1 + nz
nvar = 1 # 個人属性で求める対象の数(生鮮食品とそれ以外の日曜品とかだと2)
limy = 1e20 # 欠測とみなす数値の境界
k = 1 # トレンド成分モデルの次数
zc = k+1+ SEASONALYTY-2 + nz
pr = 4 #プロセス数
##-------------------

### data発生プロセス-----------------------------------
## 前回来店からの日数の対数,イベントの有無(ダミー変数)------
## 値引きの有無(ダミー変数)
Ztld=np.zeros((zc,TIMES,nhh))
## ----------------------
# Ztldの整形
Ztld[0,:,:] = [[1]*nhh]*TIMES
Ztld[1,:,:] = [[1]*nhh]*TIMES
Ztld[zc-2,:,:] = genfromtxt(inpt_Z_filename, delimiter=',').T #サイト訪問間隔

Ztld[(k-1)+SEASONALYTY,:,:] = hbm.standardization(Ztld[(k-1)+SEASONALYTY,:,:])
Ztld[(k-1)+SEASONALYTY+nz-1,:,:] = hbm.calc_INV(TIMES, PM, C, nhh)
Ztld[(k-1)+SEASONALYTY+nz-1,:,:] = hbm.standardization(Ztld[(k-1)+SEASONALYTY+nz-1,:,:])
print "done data shape"
###------------------------------------------------------------

## 事前分布のパラメタ-----------------
# '''step1:潜在効用サンプリング用'''
A = 0.01 * np.identity(zc) ## AはB_0の逆数
b0 = np.zeros((zc,nvar),dtype=np.float)

# step3:システムノイズの分散サンプリング用
mu0 = 0; kaps0 = 25
nu0 = 0.02; s0 = 0.02
Sita_sys0 = np.array([np.float(10)]*m)

# step5:消費者異質性の回帰パラメタHのデータ枠
m0 = np.zeros((nD,nvar),dtype=np.float)
A0 = 0.01 * np.identity(nD) #樋口本ではA

# step6:消費者異質性の回帰パラメタVのデータ枠
f0 = nvar+3
V0 = f0 * np.identity(nvar)
##------------------------------------

## 事後分布サンプリングに必要なデータの枠を作成-----------
# step1:潜在効用サンプリングのデータ枠
u = np.zeros((nhh,TIMES),dtype=np.float)

# step2:状態ベクトルの算出のデータ枠
# 処理の都合上、初期値の設定部分で設定

# step3:システムノイズの分散サンプリングのデータ枠
Sita_sys = np.array([np.float(10)*np.identity(m) for i in xrange(nhh)]) # 3次元のためmatrix化できないので、計算時にmatrix化

# step4:擬似家庭内在庫を規定するためのパラメタサンプリングのためのデータ枠
# θが2変数を持つ時はベクトルではなくmatrixにすること
Lsita_dlt = np.zeros((nhh),dtype=np.float)
Lsita_lmbd = np.zeros((nhh),dtype=np.float)
Hdlt = np.zeros((nvar,nD),dtype=np.float)
Hlmbd = np.zeros((nvar,nD),dtype=np.float)
Vsita_dlt = 0.01*np.identity(nvar)
Vsita_lmbd = 0.01*np.identity(nvar)
Sita_dlt = np.zeros((nhh),dtype=np.float)
Sita_lmbd = np.zeros((nhh),dtype=np.float)
sigma_dlt = 1.5*np.identity(nvar)
sigma_lmbd = 1.5*np.identity(nvar)
rej_dlt = np.zeros((nhh),dtype=np.float)
rej_lmbd = np.zeros((nhh),dtype=np.float)
##---------------------------------------------------------

## 初期値の設定------------------------------
# step1用
Xs = np.zeros((nhh, zc, TIMES),dtype=np.float) #人,変数数,time
sigma = 1.0

## step2用
param = hbm.FGHset(0, 1, SEASONALYTY, 0, nz)
L = 1
R = np.identity(L)
F = np.array(param['MatF'])
G = np.array(param['MatG'])
# システムモデルの分散を個人ごとに格納する枠
Q0 =np.array(param['MatQ'])

# step3用
mu = 0.
sigs = 1.
##-------------------------------------------
## 切断範囲の指定
at = np.array([[-100 if y[hh,t]==0 else 0 for t in xrange(TIMES)] for hh in xrange(nhh)])
bt = np.array([[0 if y[hh,t]==0 else 100 for t in xrange(TIMES)] for hh in xrange(nhh)])
y = None

##-------------------
udraw = np.zeros((nhh,1,RP), dtype=np.float)

def calculate_utility((hh, Xs, Ztld)):
# step1--------------------------------------------------
# uのサンプリング(ループはさせていないが個人ごとの計算)
mu = np.sum(Ztld[:,:,hh] * Xs[hh,:,:], axis = 0)
u[hh,:] = hbm.rtnorm(mu, sigma, at[hh,:], bt[hh,:])
return u[hh,:]
#------------------------------------------------------------

def Kalman_filtering((hh, u, Sita_sys, Ztld)):
## step2のシステムモデルパラメータの計算----------------------
# TAU2の最尤推定を求める数値計算------------------------------
ISW = 0; XPS=0;
#mybounds=[(1e-4,1e2),(1e-4,1e2),(1e-4,1e2),(1e-4,1e2),(1e-4,1e2)]
#LLF1 = so.fmin_l_bfgs_b(hbm.LogL, x0=tau0,
# args=(np.array(u[hh,:]), F, np.array(Ztld[:,:,hh]), G, R, limy, ISW, zc, m, TIMES, Q0),
# bounds=mybounds, approx_grad=True)
# TAU2の最尤推定

# カルマンフィルタ
#Q = hbm.Qset(Q0 ,TAU2);
XF0 = np.zeros(zc)
VF0 = np.float(10) * np.identity(zc); OSW = 1
LLF2 = hbm.KF(u[hh,:], XF0, VF0, F, Ztld[:,:,hh], G, Sita_sys[hh,:,:], R, limy, ISW, OSW, zc, TIMES)
XPS = LLF2['XPS']; XFS = LLF2['XFS']
VPS = LLF2['VPS']; VFS = LLF2['VFS']
SIG2 = LLF2['Ovar']; GSIG2 = 1
# 平滑化 ----------------------------------------------------------
LLF3 = hbm.SMO(XPS, XFS, VPS, VFS, F, GSIG2, 1, SEASONALYTY, 1, zc, TIMES)
Xs[hh,:,:] = np.array(LLF3['XSS']).T #型を合わすために無理やり変換
return Xs[hh,:,:]
#------------------------------------------------------------

def calculate_difference((hh, Xs)):
# step3の階差の計算--------------------------------------
# step3の階差計算の和の計算で使用する変数の初期化
dift = 0.
difw = 0.
difbeta = np.array([np.float(0)]*nz)

dift = sum( (Xs[hh,0,1:TIMES] - Xs[hh,0,0:TIMES-1])**2 )
difw = sum( (Xs[hh,k,1:TIMES]+sum(Xs[hh, k:(k-1)+SEASONALYTY-1, 0:TIMES-1]))**2 )
for d in xrange(nz):
difbeta[d] = sum( (Xs[hh, (k-1)+SEASONALYTY+d, 1:TIMES]
- Xs[hh, (k-1)+ SEASONALYTY+d, 0:TIMES-1])**2 )
# step3--------------------------------------
Sita_sys[hh,0,0] = invgamma.rvs((nu0+TIMES)/2, scale=(s0+dift)/2, size=1)[0]
Sita_sys[hh,1,1] = invgamma.rvs((nu0+TIMES)/2, scale=(s0+difw)/2, size=1)[0]
for d in xrange(nz):
Sita_sys[hh,2+d,2+d] = invgamma.rvs((nu0+TIMES)/2, scale=(s0+difbeta[d])/2, size=1)[0]
return Sita_sys[hh,:,:]
#--------------------------------------------

def calculate_difference_m(Xs): # 人間すべてに対応
# step3の階差の計算--------------------------------------
# step3の階差計算の和の計算で使用する変数の初期化
dift = np.zeros(nhh, dtype=np.float)
difw = np.zeros(nhh, dtype=np.float)
difbeta = np.zeros((nhh,nz), dtype=np.float)

dift = np.sum( (Xs[:,0,1:TIMES] - Xs[:,0,0:TIMES-1])**2, axis=1 )
difw = np.sum( (Xs[:,k,1:TIMES]+np.sum(Xs[:, k:(k-1)+SEASONALYTY-1, 0:TIMES-1], axis=1))**2, axis=1 )
for d in xrange(nz):
difbeta[:,d] = np.sum( (Xs[:, (k-1)+SEASONALYTY+d, 1:TIMES]
- Xs[:, (k-1)+ SEASONALYTY+d, 0:TIMES-1])**2, axis=1 )
# step3--------------------------------------
Sita_sys[:,0,0] = invgamma.rvs((nu0+TIMES)/2, scale=(s0+dift)/2, size=nhh)
Sita_sys[:,1,1] = invgamma.rvs((nu0+TIMES)/2, scale=(s0+difw)/2, size=nhh)
for d in xrange(nz):
Sita_sys[:,2+d,2+d] = invgamma.rvs((nu0+TIMES)/2, scale=(s0+difbeta[:,d])/2, size=nhh)
return Sita_sys[:,:,:]
#--------------------------------------------

def calculate_spending_habits_param_delta((hh, u, Xs, Sita_dlt, Hdlt, Vsita_dlt, Ztld)): # step4のΘ_δの計算
### step4--------------------------------------
## '''dlt側の計算'''
# step4のθの事後分布カーネルの第一項の和計算時使用する変数の初期化
Lsita = 0.
# 効用値の誤差計算(θの尤度計算)
Lsita = sum( (u[hh,:] - np.diag(np.dot(Ztld[:,:,hh].T, Xs[hh,:,:])))**2 )
# 現状のθを確保する
old_sita_dlt = Sita_dlt[hh]
# 新しいθをサンプリング(酔歩サンプリング)
new_sita_dlt = Sita_dlt[hh] + ss.norm.rvs(loc=0, scale=sigma_dlt,size=1)[0]

# 尤度の計算(対数尤度の場合はヤコビアンで調整)
new_Lsita_dlt = Lsita + hbm.Nkernel(new_sita_dlt, Hdlt, D[hh,:], Vsita_dlt)
new_Lsita_dlt = math.exp(-0.5*new_Lsita_dlt)
#new_Lsita_dlt = new_Lsita_dlt*(1/math.exp(new_Lsita_dlt))
old_Lsita_dlt = Lsita + hbm.Nkernel(old_sita_dlt, Hdlt, D[hh,:], Vsita_dlt)
old_Lsita_dlt = math.exp(-0.5*old_Lsita_dlt)
#old_Lsita_dlt = old_Lsita_dlt*(1/math.exp(old_Lsita_dlt))
if old_Lsita_dlt == 0: old_Lsita_dlt=1e-6

# MHステップ
alpha = min(1, new_Lsita_dlt/old_Lsita_dlt)
if alpha==None: alpha = -1
uni = ss.uniform.rvs(loc=0 , scale=1, size=1)
if uni < alpha and new_sita_dlt > 0:
Sita_dlt[hh] = new_sita_dlt
else:
rej_dlt[hh] = rej_dlt[hh] + 1
return Sita_dlt[hh], rej_dlt[hh]

def calculate_spending_habits_param_delta_m((u, Xs, Sita_dlt, Hdlt, Vsita_dlt, Ztld, rej)): # step4-多変量版
### step4--------------------------------------
## '''dlt側の計算'''
# step4のθの事後分布カーネルの第一項の和計算時使用する変数の初期化
Lsita = np.zeros(nhh,dtype=np.float)
# 効用値の誤差計算(θの尤度計算)
Lsita = [np.sum( (u[hh,:] - np.sum(Ztld[:,:,hh]*Xs[hh,:,:], axis = 0) )**2 ) for hh in xrange(nhh)]
# 現状のθを確保する
old_sita_dlt = Sita_dlt[:]
# 新しいθをサンプリング(酔歩サンプリング)
new_sita_dlt = Sita_dlt[:] + ss.norm.rvs(loc=0, scale=sigma_dlt, size=nhh)

# 尤度の計算(対数尤度の場合はヤコビアンで調整)
new_Lsita_dlt = Lsita + hbm.Nkernel(new_sita_dlt, Hdlt, D[:,:], Vsita_dlt)
new_Lsita_dlt = np.exp(-0.5*new_Lsita_dlt)
#new_Lsita_dlt = new_Lsita_dlt*(1/math.exp(new_Lsita_dlt))
old_Lsita_dlt = Lsita + hbm.Nkernel(old_sita_dlt, Hdlt, D[:,:], Vsita_dlt)
old_Lsita_dlt = np.exp(-0.5*old_Lsita_dlt)
#old_Lsita_dlt = old_Lsita_dlt*(1/math.exp(old_Lsita_dlt))
old_Lsita_dlt[old_Lsita_dlt == 0] = 1e-6

# MHステップ
alpha = np.array(new_Lsita_dlt/old_Lsita_dlt)
alpha = alpha*(alpha<1) # 1より大きければここでは0が入る.1より小さければそのままの値.
alpha[alpha==0] = 1 # 1より大きいものは1とする,をここで実行
alpha[alpha!=alpha] = -1 # NaNであった部分を-1とする
uni = ss.uniform.rvs(loc=0 , scale=1, size=nhh)
tmp = alpha - uni
Sita_tmp = new_sita_dlt*(tmp>0) + old_sita_dlt*(tmp<=0) # MHの採択条件
Sita_dlt = Sita_tmp*(Sita_tmp>=0) + old_sita_dlt*(Sita_tmp<0) # この分析ではδは非負であるという条件
tmp = Sita_dlt - old_sita_dlt #
rej += 1*(tmp[0,:]==0) # 計算の結果が初期値と同じであればリジェクトに1をプラス
#return np.ndarray.flatten(Sita_dlt), np.ndarray.flatten(rej)
return Sita_dlt[0,:], rej

def calculate_spending_habits_param_lambd((hh, u, Xs, Sita_lmbd, Hlmbd, Vsita_lmbd, Ztld)): # step4のΘ_λの計算
### step4--------------------------------------
## '''lmbd側の計算'''
# step4のθの事後分布カーネルの第一項の和計算時使用する変数の初期化
Lsita = 0.
# 効用値の誤差計算(θの尤度計算)
Lsita = sum( (u[hh,:] - np.diag(np.dot(Ztld[:,:,hh].T, Xs[hh,:,:])))**2 )
# 現状のθを確保する
old_sita_lmbd = Sita_lmbd[hh]
# 新しいθをサンプリング(酔歩サンプリング)
new_sita_lmbd = Sita_lmbd[hh] + ss.norm.rvs(loc=0, scale=sigma_lmbd, size=1)[0]

# 尤度の計算(対数尤度の場合はヤコビアンで調整)
new_Lsita_lmbd = Lsita + hbm.Nkernel(new_sita_lmbd, Hlmbd, D[hh,:], Vsita_lmbd)
new_Lsita_lmbd = math.exp(-0.5*new_Lsita_lmbd)
old_Lsita_lmbd = Lsita + hbm.Nkernel(old_sita_lmbd, Hlmbd, D[hh,:], Vsita_lmbd)
old_Lsita_lmbd = math.exp(-0.5*old_Lsita_lmbd)
if old_Lsita_lmbd == 0: old_Lsita_lmbd=1e-6

# MHステップ
alpha = min(1, new_Lsita_lmbd/old_Lsita_lmbd)
if alpha==None: alpha = -1
uni = ss.uniform.rvs(loc=0, scale=1, size=1)
if uni < alpha:
Sita_lmbd[hh] = new_sita_lmbd
else:
rej_lmbd[hh] = rej_lmbd[hh] + 1
return Sita_lmbd[hh], rej_lmbd[hh]

def calculate_spending_habits_param_lambd_m((u, Xs, Sita_lmbd, Hlmbd, Vsita_lmbd, Ztld, rej)): # step4-多変量
### step4--------------------------------------
## '''lmbd側の計算'''
# step4のθの事後分布カーネルの第一項の和計算時使用する変数の初期化
Lsita = np.zeros(nhh,dtype=np.float)
# 効用値の誤差計算(θの尤度計算)
Lsita = [np.sum( (u[hh,:] - np.sum(Ztld[:,:,hh]*Xs[hh,:,:], axis = 0) )**2 ) for hh in xrange(nhh)]
# 現状のθを確保する
old_sita_lmbd = Sita_lmbd[:]
# 新しいθをサンプリング(酔歩サンプリング)
new_sita_lmbd = Sita_lmbd[:] + ss.norm.rvs(loc=0, scale=sigma_lmbd, size=nhh)

# 尤度の計算(対数尤度の場合はヤコビアンで調整)
new_Lsita_lmbd = Lsita + hbm.Nkernel(new_sita_lmbd, Hlmbd, D[:,:], Vsita_lmbd)
new_Lsita_lmbd = np.exp(-0.5*new_Lsita_lmbd)
old_Lsita_lmbd = Lsita + hbm.Nkernel(old_sita_lmbd, Hlmbd, D[:,:], Vsita_lmbd)
old_Lsita_lmbd = np.exp(-0.5*old_Lsita_lmbd)
old_Lsita_lmbd[old_Lsita_lmbd == 0] = 1e-6

# MHステップ
alpha = np.array(new_Lsita_lmbd/old_Lsita_lmbd)
alpha = alpha*(alpha<1) # 1より大きければここでは0が入る.1より小さければそのままの値.
alpha[alpha==0] = 1 # 1より大きいものは1とする,をここで実行
alpha[alpha!=alpha] = -1 # NaNであった部分を-1とする
uni = ss.uniform.rvs(loc=0 , scale=1, size=nhh)
tmp = alpha - uni
Sita_lmbd = new_sita_lmbd*(tmp>0) + old_sita_lmbd*(tmp<=0) # MHの採択条件
rej += 1*(tmp[0,:]<=0) # 計算の結果が初期値と同じであればリジェクトに1をプラス
return Sita_lmbd[0,:], rej

## step7-----------------------------------------------------------
def calc_INV((hh, dlt, lmbd)):
INV = np.zeros(TIMES,dtype = np.double)
for t in xrange(1,TIMES):
if abs(INV[t-1]) < 1e-6 : Cons = 0.0
else:
denom = dlt[hh]*C[hh] + (INV[t-1])**lmbd[hh]
if denom == 0.0: denom = 1e-6
Cons = INV[t-1]*dlt[hh]*C[hh]/denom
INV[t] = INV[t-1] + PM[hh, t-1]- Cons
return INV

def calc_INV_m((dlt, lmbd)):
INV = np.zeros((TIMES,nhh), dtype = np.double)
Cons = np.zeros(nhh, dtype = np.double)
for t in xrange(1,TIMES):
tmp = INV[t-1,:]*(abs(INV[t-1,:])>=1e-6) + 1e-6*(abs(INV[t-1,:])<1e-6) # INVが1e-6以上ではINVの値を採用、未満の場合は分子が小さくなりすぎる(warnigが出る)ので1e-6とする
denom = np.array(dlt[:]*C[:] + (tmp)**lmbd[:])
denom[denom == 0.0] = 1e-6
Cons = (INV[t-1,:]*dlt*C/denom)*(abs(INV[t-1,:])>=1e-6) # INV[t-1]が1e-6以下の時は0とする
INV[t] = INV[t-1,:] + PM[:,t-1]- Cons
return INV
## ----------------------------------------------------------------

def pool_map(func, itr, args):
pool = Pool(processes=pr)
result = pool.map( func, ((i, args) for i in range(itr)) )
pool.close()
pool.join()
return result
#--------------------------------------------
print "done ini"
## サンプリングのループ
if __name__ == "__main__":
rej_dlt = [0]*nhh
rej_lmbd = [0]*nhh

#start = time.clock()
for nd in xrange(RP):
# step1--calculate utility
pool = Pool(processes=pr)
u = np.array( pool.map(calculate_utility, ((hh, Xs, Ztld) for hh in xrange(nhh))) )
pool.close()
pool.join()
udraw[:,:,nd] = np.matrix(u[:,10]).T
# step2--calculate implicit system v variable=Xa(do kalman filter)
pool = Pool(processes=pr)
Xs = np.array( pool.map(Kalman_filtering, ((hh, u, Sita_sys, Ztld) for hh in xrange(nhh))) )
pool.close()
pool.join()
# step3--calculate difference
Sita_sys = calculate_difference_m(Xs)
# step4--calculate_spending_habits_param_delta=Sita_dlt
Sita_dlt, rej_dlt = calculate_spending_habits_param_delta_m((u, Xs, Sita_dlt, Hdlt, Vsita_dlt, Ztld, rej_dlt))
# step4--calculate_spending_habits_param_lambd=Sita_lmbd
Sita_lmbd, rej_lmbd = calculate_spending_habits_param_lambd_m((u, Xs, Sita_lmbd, Hlmbd, Vsita_lmbd, Ztld, rej_lmbd))
### step5--------------------------------------
## dlt側の算出----
# 多変量正規分布のパラメタの算出
D2 = np.dot(D.T, D)
D2pA0 = D2 + A0
Hhat_dlt = np.ndarray.flatten( np.dot(np.dot(inverse(D2), D.T) , Sita_dlt.T) )
Dtld = np.dot( inverse(D2pA0) , (np.dot(D2, Hhat_dlt) + np.dot(A0, np.ndarray.flatten(m0))) )
rtld = np.ndarray.flatten(Dtld)
sig = np.array( np.kron(Vsita_dlt, inverse(D2pA0)).T )
# 多変量正規分布でサンプリング
Hdlt = np.ndarray.flatten( hbm.randn_multivariate(rtld, np.matrix(sig), n=nvar) )
##-----------------
## lmbd側の算出----
# 多変量正規分布のパラメタの算出
Hhat_lmbd = np.ndarray.flatten( np.dot( np.dot(inverse(D2), D.T) , Sita_lmbd.T) )
Dtld = np.dot( inverse(D2pA0) , (np.dot(D2, Hhat_lmbd) + np.dot(A0, np.ndarray.flatten(m0))) )
rtld = np.ndarray.flatten(Dtld)
sig = np.array( np.kron(Vsita_lmbd, inverse(D2pA0)).T )
# 多変量正規分布でサンプリング
Hlmbd = np.ndarray.flatten( hbm.randn_multivariate(rtld, np.matrix(sig), n=nvar) )
##-----------------
#--------------------------------------------

### step6--------------------------------------
##dlt側の算出
# 逆ウィッシャート分布のパラメタの算出
#div = np.array(Sita_dlt) - np.dot( D, np.ndarray.flatten(Hdlt).T ).T
div = np.array(Sita_dlt) - np.dot( D, Hdlt.T ).T
S = np.dot(div, div.T) # 上の計算でdivは横ベクトルになるのでSをスカラにするためにTは後ろ
# 逆ウィッシャート分布でサンプリング
Vsita_dlt = hbm.invwishartrand(f0 + nhh, V0 + S)
##------------
##lmbd側の算出
# 逆ウィッシャート分布のパラメタの算出
div = np.array(Sita_lmbd) - np.dot( D, Hlmbd.T ).T
S = np.dot(div, div.T)
# 逆ウィッシャート分布でサンプリング
Vsita_lmbd = hbm.invwishartrand(f0 + nhh, V0 + S)
##------------
#--------------------------------------------

## step7-------------------------------------------
##Zの消費金額を更新
INV = calc_INV_m((Sita_dlt, Sita_lmbd))
Ztld[(k-1)+SEASONALYTY+nz-1,:,:] = hbm.standardization(INV)
##-------------------------------------------
print "nd=%d", nd
#stop = time.clock()
#print stop - start
##-------------------------------------------
# cProfile.run("worker(int(sys.argv[1]), int(sys.argv[2]))")
## output area---------------------------------------------------------
outf = open('test_udraw.csv','w')
nums = []
for hh in xrange(nhh):
prsd = np.std(udraw[hh,0, 0:math.ceil(0.1*RP)]); posd = np.std(udraw[hh,0,RP-math.floor(0.5*RP):])
prav = np.average(udraw[hh,0,0:math.ceil(0.1*RP)]); poav = np.average(udraw[hh,0,RP-math.floor(0.5*RP):])
prf = np.sum(abs(udraw[hh,0,0:math.ceil(0.1*RP)])); pof = np.sum(abs(udraw[hh,0,RP-math.floor(0.5*RP):]))
Z = math.sqrt(RP)*(prav - poav)/math.sqrt( prsd**2/0.1 + posd**2/0.5 )
outf.write("UsrNon="+ str(hh) + ",Z="+ str(Z) +"\n")
for t in xrange(1):
nums = [str(udraw[hh,t,p/keep]) for p in xrange(RP) if p%keep == 0 ]
outf.write((",".join(nums))+"\n")
outf.close()
## --------------------------------------------------------------------
## output area---------------------------------------------------------
outf = open('test_crm.csv','w')
nums = []
for hh in xrange(nhh):
outf.write("UsrNon="+ str(hh) +"\n")
for r in xrange(zc):
nums = [str(Xs[hh,r,t]) for t in xrange(TIMES)]
outf.write((",".join(nums))+"\n")
outf.close()
## --------------------------------------------------------------------
## output area---------------------------------------------------------
outf01 = open('test_u.csv','w')
outf02 = open('test_INV.csv','w')
outf03 = open('test_rej_dlt.csv','w')
outf04 = open('test_rej_lmbd.csv','w')
outf05 = open('test_sita_dlt.csv','w')
outf06 = open('test_sita_lmbd.csv','w')
nums01 = []; nums02 = []; nums03 = []; nums04 = []; nums05 = []; nums06 = []
for hh in xrange(nhh):
#順序付きで格納
nums01 = [str(u[hh,t]) for t in xrange(TIMES)]
nums02 = [str(INV[t,hh]) for t in xrange(TIMES)]
nums03 = [str(rej_dlt[hh])]
nums04 = [str(rej_lmbd[hh])]
nums05 = [str(Sita_dlt[hh])]
nums06 = [str(Sita_lmbd[hh])]

# 書き込み
outf01.write((",".join(nums01))+"\n")
outf02.write((",".join(nums02))+"\n")
outf03.write((",".join(nums03))+"\n")
outf04.write((",".join(nums04))+"\n")
outf05.write((",".join(nums05))+"\n")
outf06.write((",".join(nums06))+"\n")
outf01.close()
outf02.close()
outf03.close()
outf04.close()
outf05.close()
outf06.close()
## --------------------------------------------------------------------

print("done all")


次に枝葉の部分の関数ファイルです。


one_to_one_mrk.py

# coding: utf8


import numpy.random as npr
from numpy.oldnumeric.linear_algebra import inverse
import scipy
from scipy import linalg
import random
import scipy.stats.distributions as dis
from scipy.stats import uniform,chi2,invgamma
import scipy.stats as ss
import numpy as np
import math
from scipy.linalg import inv,cholesky,det

### 関数の定義------------------------------------
# 2項プロビットモデルのベイズ推定(step1用)
def rtnorm(mu, sigma, a, b):
lngth = len(mu)
FA = dis.norm.cdf(a, loc=mu, scale=sigma, size = lngth)
FB = dis.norm.cdf(b, loc=mu, scale=sigma, size = lngth)
result = np.array([
dis.norm.ppf(
np.dot(ss.uniform.rvs(loc=0, scale=1,size=len(np.matrix(mu[t]))),
(FB[t] - FA[t]) + FA[t]),
loc=mu[t], scale=sigma, size=len(np.matrix(mu[t])) ) #percent point = Q値
for t in xrange(lngth)
]
)
result[result == -np.inf] = -100
result[result == np.inf] = 100
return result[:,0]

# Zの値を基準化する関数(step1用)
def standardization(z):
return np.log( 0.001 + (z - np.min(z))/(np.max(z) - np.min(z))*(0.999 - 0.001) )

# 単変量正規分布のカーネル部分の乗数の部分の計算(step4用)
def Nkernel(sita, H, D, Vsita):
if Vsita == 0: Vsita=1e-6
return ((sita - np.dot(H,D.T))**2)/Vsita

# 多変量正規分布のカーネル部分の乗数の部分の計算(step4用)
def NMkernel(sita, H, D, Vsita):
res = sita - np.dot(H.T, D)
return np.dot( np.dot(res.T, inverse(Vsita)), res )

### 季節調整モデルの状態空間表現の行列設定
def FGHset(al, k, p, q, nz): #alはARモデルのαベクトル、k;p;qはトレンド、季節、AR
m = k + p + q + nz -1
if(q>0): G = np.zeros((m,3+nz)) # 状態モデルでトレンド、季節、ARの3つを含む場合
else: G = np.zeros((m,2+nz)) #AR成分を含まない場合(q=0)
F = np.zeros((m,m))
#H = matrix(0,1,m) #Hの代わりにZtldを使うので、Hは不要

## トレンドモデルのブロック行列の構築
G[0,0] = 1
if k==1: F[0,0] = 1
if k==2: F[0,0] = 2; F[0,1] = -1; F[1,0] = 1
if k==3: F[0,0] = 3; F[0,1] = -3; F[0,2] = 1; F[1,0] = 1; F[2,1] = 1
LS = k
NS = 2

## 季節調整成分のブロック行列の構築
G[LS, NS-1] = 1
for i in xrange(p-1): F[LS, LS+i] = -1
for i in xrange(p-2): F[LS+i+1, LS+i] = 1

## Z成分のブロック行列の構築
LS = LS + p -1
NS = 2
for i in xrange(nz): F[LS+i, LS+i] = 1
for i in xrange(nz): G[LS+i, NS+i] = 1

if q>0:
NS = NS +1
G[LS, NS-1] = 1
for i in xrange(q): F[LS, LS+i-1] = al[i]
if q>1:
for i in xrange(q-1): F[LS+i, LS+i-1] = 1

## シスムモデルの分散共分散行列Qの枠の算出
Q = np.identity(NS+nz)

return {'m':m, 'MatF':F, 'MatG':G, 'MatQ':Q}

# 状態空間表現における行列Qの設定------------------------
def Qset(Q0,parm):
NS = len(Q0)
Q = Q0
# シスムモデルの分散共分散行列Qの枠の算出
for i in xrange(NS): Q[i,i] = parm[i]
return np.array(Q)

# カルマンフィルタの関数 ------------------------------------------
def KF(y, XF0, VF0, F, H, G, Q, R, limy, ISW, OSW, m, N):
if OSW == 1:
XPS = np.zeros((N,m),dtype=np.float); XFS = np.zeros((N,m),dtype=np.float)
VPS = np.zeros((N,m,m),dtype=np.float); VFS = np.zeros((N,m,m),dtype=np.float)
XF = XF0; VF = VF0; NSUM = 0.0; SIG2 = 0.0; LDET = 0.0
for n in xrange(N):
# 1期先予測
XP = np.ndarray.flatten( np.dot(F, XF.T) ) #2週目から縦ベクトルになってしまうので、常に横ベクトルに変換
VP = np.dot( np.dot(F, VF), F.T ) + np.dot( np.dot(G, Q), G.T)
# フィルタ
# Rは操作しなければ縦ベクトル。pythonは横ベクトルになるので注意!
if y[n] < limy:
NSUM = NSUM + 1
B = np.dot( np.dot(H[:,n], VP), H[:,n].T) + R # Hは数学的には横ベクトル
B1 = inverse(B) # nvar次元の縦ベクトル
K = np.matrix(np.dot(VP, H[:,n].T)).T * np.matrix(B1) # Kは縦ベクトルになる(matrix)
e = np.array(y[n]).T - np.dot(H[:,n], XP.T) # nvar次元の縦ベクトル
XF = np.array(XP) + np.array( K * np.matrix(e) ).T # 横ベクトル
VF = np.array(VP) - np.array( K* np.matrix(H[:,n]) * VP)
SIG2 = SIG2 + np.ndarray.flatten(np.array( np.matrix(e) * np.matrix(B1) * np.matrix(e).T ))[0] # 1次元でも計算できるようにmatrixにする
LDET = LDET + math.log(linalg.det(B))
else:
XF = XP; VF = VP
if OSW == 1:
XPS[n,:] = XP; XFS[n,:] = XF; VPS[n,:,:] = VP; VFS[n,:,:] = VF
SIG2 = SIG2 / NSUM
if ISW == 0:
FF = -0.5 * (NSUM * (math.log(2 * np.pi * SIG2) + 1) + LDET)
else:
FF = -0.5 * (NSUM * (math.log(2 * np.pi) + SIG2) + LDET)
if OSW == 0:
return {'LLF':FF, 'Ovar':SIG2}
if OSW == 1:
return {'XPS':XPS, 'XFS':XFS, 'VPS':VPS, 'VFS':VFS, 'LLF':FF, 'Ovar':SIG2}

# 平滑化の関数 ----------------------------------------------------
def SMO(XPS, XFS, VPS, VFS, F, GSIG2, k, p, q, m, N):
XSS = np.zeros((N,m),dtype=np.float); VSS = np.zeros((N,m,m),dtype=np.float)
XS1 = XFS[N-1,:]; VS1 = VFS[N-1,:,:]
XSS[N-1,:] = XS1; VSS[N-1,:,:] = VS1
for n1 in xrange(N-1):
n = (N-1) - n1; XP = XPS[n,:]; XF = XFS[n-1,:]
VP = VPS[n,:,:]; VF = VFS[n-1,:,:]; VPI = inverse(VP)
A = np.dot( np.dot(VF, F.T), VPI)
XS2 = XF + np.dot(A, (XS1 - XP))
VS2 = VF + np.dot( np.dot(A, (VS1 - VP)), A.T )
XS1 = XS2; VS1 = VS2
XSS[n-1,:] = XS1; VSS[n-1,:,:] = VS1
return {'XSS':XSS, 'VSS':VSS}

# TAU2xの対数尤度関数の定義 ----------------------------------------
def LogL(parm, *args):
y=args[0]; F=args[1]; H=args[2]; G=args[3]; R=args[4]; limy=args[5]
ISW=args[6]; k=args[7]; m=args[8]; N=args[9]; Q0=args[10]
Q = Qset(Q0 ,parm)
XF0 = np.array([0]*k); VF0 = np.array(10 * np.identity(k)); OSW = 0
LLF = KF(y, XF0, VF0, F, H, G, Q, R, limy, ISW, OSW, k, N)
LL = LLF['LLF']
return -LL # optimezeが最小化関数なので、対数尤度にマイナスをかけたものを返す

# 多変量正規分布の発生関数の定義--------------------------------------
def randn_multivariate(mu,Sigma,n=1):
X = np.random.randn(n,len(mu))
A = linalg.cholesky(Sigma)
Y = np.dot(np.array(X),np.array(A)) + mu
return Y

# 逆ウィッシャート関数の定義----------------------------------------------
def invwishartrand_prec(nu,phi):
return inv(wishartrand(nu,phi))

def invwishartrand(nu, phi):
return inv(wishartrand(nu, inv(phi)))

def wishartrand(nu, phi):
dim = phi.shape[0]
chol = cholesky(phi)
foo = np.zeros((dim,dim))

for i in xrange(dim):
for j in xrange(i+1):
if i == j:
foo[i,j] = np.sqrt(chi2.rvs(nu-(i+1)+1))
else:
foo[i,j] = npr.normal(0,1)
return np.dot(chol, np.dot(foo, np.dot(foo.T, chol.T)))
# -------------------------------------------------------

# 疑似家庭内在庫金額の計算----------------------------------------------
def calc_INV(TIMES, PM, C, nhh):
## 顧客の消費傾向Data
INV = np.zeros((TIMES,nhh),dtype = np.double)
dlt = np.array([np.float(1.0)]*nhh)
lmbd = np.array([np.float(1.0)]*nhh)
for t in xrange(1,TIMES):
denom = np.array(dlt[:]*C[:] + (INV[t-1,:])**lmbd[:])
denom[denom==0.0] = 1e-6
Cons = INV[t-1,:]*dlt[:]*C[:]/denom
INV[t,:] = INV[t-1,:] + PM[:,t-1] - Cons
return np.array(INV)
# -------------------------------------------------------


簡単のためロイヤリティ算出の際の説明変数の数を減らしています。

この点、書籍の内容は異なっていますので注意してください。

また書籍には具体的な記載が無かったのですが、MCMCの収束判定も入れてあります。

判定はGewekeの診断によっています。

MCMCの概要に関しましてはiAnalysis ~おとうさんの解析日記~ MCMC(Markov Chain Monte Carlo、マルコフ連鎖モンテカルロ法)についてを参照していただきたく存じます。

テストで使うデータが欲しい、という方は当該書籍のユーザサポートページから取得されるのが手軽で良いかと思います。


まとめ

アウトプットの掲載はめんどくさかったので控えました。

怠慢ですね。すみません。。。

ご興味のある方は実務で使用してみてください。

(私の実装もかなりまずいと思うのですが、pythonでは遅くて実用は厳しいため、コンパイラ言語で書かれた方が良いかと思います^^;)