Comme le titre l'indique, j'ai mis en place la fidélisation de la clientèle en référence aux livres suivants.
Cependant, il peut y avoir des parties étranges car il s'agit d'une combinaison de choses faites dans le passé ... (Cela a bien fonctionné dans mon environnement local.) Oh? S'il y a quelque chose qui semble être une surprise, veuillez commenter.
J'écrirai le code ci-dessous. Veuillez noter qu'il y a un petit montant.
Le premier est le code du fichier sur le côté principal.
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
#Définition du chemin du fichier d'importation de données
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'
###Entrée de données-------------------
#Acquisition de données D (données nécessaires pour estimer les paramètres individuels en utilisant la communalité)
D = genfromtxt(inpt_D_filename, delimiter=',')
#Acquisition des données C (consommation moyenne journalière)
C = genfromtxt(inpt_C_filename, delimiter=',')
#Acquisition des données PM (montant par date d'achat pour obtenir INV)
PM = genfromtxt(inpt_payment_filename, delimiter=',')
#y Obtenir des données (une distribution normale coupée par indicateur de visite,Trouver b)
y = genfromtxt(inpt_visit_filename, delimiter=',')
print "done input"
### ------------------------------
#Occurrence de graines aléatoires
random.seed(555)
##--Définition de constante-------
TIMES = len(PM[0,:])
nhh = len(PM[:,0])
SEASONALYTY = 7 #Période de fluctuation saisonnière
RP = 50000 #Nombre d'échantillons prélevés par MCMC
keep = 10 #Afficher l'intervalle de données
nz = 2 #Nombre de variables explicatives dans la formule de l'utilitaire de visite en magasin(Z_nombre de i)
nD = len(D[0,:]) #Longueur de vecteur de D
m = 1 + 1 + nz
nvar = 1 #Nombre d'objets à calculer par attribut personnel (2 pour les aliments frais et autres produits du dimanche)
limy = 1e20 #Limite des valeurs numériques à considérer comme manquantes
k = 1 #Ordre du modèle de composant de tendance
zc = k+1+ SEASONALYTY-2 + nz
pr = 4 #Nombre de processus
##-------------------
###processus de génération de données-----------------------------------
##Journal du nombre de jours depuis la dernière visite,Présence / absence d'événement (variable fictive)------
##Avec ou sans remise (variable factice)
Ztld=np.zeros((zc,TIMES,nhh))
## ----------------------
#Mise en forme Ztld
Ztld[0,:,:] = [[1]*nhh]*TIMES
Ztld[1,:,:] = [[1]*nhh]*TIMES
Ztld[zc-2,:,:] = genfromtxt(inpt_Z_filename, delimiter=',').T #Intervalle de visite sur site
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"
###------------------------------------------------------------
##Paramètres de distribution antérieurs-----------------
# '''étape 1: pour l'échantillonnage d'utilité latente'''
A = 0.01 * np.identity(zc) ##A est B_Inverse de 0
b0 = np.zeros((zc,nvar),dtype=np.float)
#étape 3: pour l'échantillonnage distribué du bruit du système
mu0 = 0; kaps0 = 25
nu0 = 0.02; s0 = 0.02
Sita_sys0 = np.array([np.float(10)]*m)
#étape 5: base de données du paramètre de régression H de l'hétérogénéité des consommateurs
m0 = np.zeros((nD,nvar),dtype=np.float)
A0 = 0.01 * np.identity(nD) #À Higuchimoto, A
#étape 6: base de données du paramètre de régression V de l'hétérogénéité des consommateurs
f0 = nvar+3
V0 = f0 * np.identity(nvar)
##------------------------------------
##Créer une base de données requise pour l'échantillonnage de distribution postérieure-----------
#étape 1: base de données pour l'échantillonnage d'utilité latente
u = np.zeros((nhh,TIMES),dtype=np.float)
# step2:Trame de données pour le calcul du vecteur d'état
#Pour faciliter le traitement, définissez dans la partie de réglage de la valeur initiale
#étape 3: trame de données pour l'échantillonnage distribué du bruit du système
Sita_sys = np.array([np.float(10)*np.identity(m) for i in xrange(nhh)]) #Puisqu'il s'agit de 3D, il ne peut pas être converti en matrice, il est donc converti en matrice au moment du calcul.
#étape 4: base de données pour l'échantillonnage des paramètres pour spécifier l'inventaire des pseudo-ménages
#Lorsque θ a 2 variables, faites-en une matrice au lieu de vecteur
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)
##---------------------------------------------------------
##Réglage de la valeur initiale------------------------------
#Pour l'étape 1
Xs = np.zeros((nhh, zc, TIMES),dtype=np.float) #Homme,Nombre de variables,time
sigma = 1.0
##Pour 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'])
#Un cadre pour stocker la distribution du modèle de système pour chaque individu
Q0 =np.array(param['MatQ'])
#Pour step3
mu = 0.
sigs = 1.
##-------------------------------------------
##Spécification de la plage de coupe
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 échantillonnage(Pas de boucles, mais des calculs individuels)
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)):
##Calcul des paramètres du modèle système à l'étape 2----------------------
#Calcul numérique pour trouver l'estimation la plus probable de 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)
#Estimation la plus probable de TAU2
#Filtre de Kalman
#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
#Lissage----------------------------------------------------------
LLF3 = hbm.SMO(XPS, XFS, VPS, VFS, F, GSIG2, 1, SEASONALYTY, 1, zc, TIMES)
Xs[hh,:,:] = np.array(LLF3['XSS']).T #Conversion forcée pour correspondre au type
return Xs[hh,:,:]
#------------------------------------------------------------
def calculate_difference((hh, Xs)):
#Calcul de la différence à l'étape 3--------------------------------------
#Initialisation des variables utilisées dans le calcul de la somme du calcul de la différence à l'étape 3
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): #Compatible avec tous les humains
#Calcul de la différence à l'étape 3--------------------------------------
#Initialisation des variables utilisées dans le calcul de la somme du calcul de la différence à l'étape 3
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)): #Θ de l'étape 4_Calcul de δ
### step4--------------------------------------
## '''Calcul côté DLT'''
#Distribution postérieure de θ à l'étape 4 Initialisation des variables utilisées lors du calcul de la somme du premier terme du noyau
Lsita = 0.
#Calcul d'erreur de la valeur d'utilité(Calcul de probabilité de θ)
Lsita = sum( (u[hh,:] - np.diag(np.dot(Ztld[:,:,hh].T, Xs[hh,:,:])))**2 )
#Sécuriser le courant θ
old_sita_dlt = Sita_dlt[hh]
#Échantillonnage du nouveau θ (échantillonnage à pied ivre)
new_sita_dlt = Sita_dlt[hh] + ss.norm.rvs(loc=0, scale=sigma_dlt,size=1)[0]
#Calcul de la vraisemblance (ajusté par Jacobien pour la vraisemblance logarithmique)
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
#Étape 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-Version multivariée
### step4--------------------------------------
## '''Calcul côté DLT'''
#Distribution postérieure de θ à l'étape 4 Initialisation des variables utilisées lors du calcul de la somme du premier terme du noyau
Lsita = np.zeros(nhh,dtype=np.float)
#Calcul d'erreur de la valeur d'utilité(Calcul de probabilité de θ)
Lsita = [np.sum( (u[hh,:] - np.sum(Ztld[:,:,hh]*Xs[hh,:,:], axis = 0) )**2 ) for hh in xrange(nhh)]
#Sécuriser le courant θ
old_sita_dlt = Sita_dlt[:]
#Échantillonnage du nouveau θ (échantillonnage à pied ivre)
new_sita_dlt = Sita_dlt[:] + ss.norm.rvs(loc=0, scale=sigma_dlt, size=nhh)
#Calcul de la vraisemblance (ajusté par Jacobien pour la vraisemblance logarithmique)
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
#Étape MH
alpha = np.array(new_Lsita_dlt/old_Lsita_dlt)
alpha = alpha*(alpha<1) #S'il est supérieur à 1, 0 est entré ici.S'il est inférieur à 1, la valeur reste la même.
alpha[alpha==0] = 1 #1 est supérieur à 1,Courez ici
alpha[alpha!=alpha] = -1 #La partie qui était NaN-Définir sur 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) #Conditions d'adoption de MH
Sita_dlt = Sita_tmp*(Sita_tmp>=0) + old_sita_dlt*(Sita_tmp<0) #La condition que δ est non négatif dans cette analyse
tmp = Sita_dlt - old_sita_dlt #
rej += 1*(tmp[0,:]==0) #Si le résultat du calcul est le même que la valeur initiale, ajoutez 1 au rejet
#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)): #Θ de l'étape 4_Calcul de λ
### step4--------------------------------------
## '''Calcul côté lmbd'''
#Distribution postérieure de θ à l'étape 4 Initialisation des variables utilisées lors du calcul de la somme du premier terme du noyau
Lsita = 0.
#Calcul d'erreur de la valeur d'utilité(Calcul de probabilité de θ)
Lsita = sum( (u[hh,:] - np.diag(np.dot(Ztld[:,:,hh].T, Xs[hh,:,:])))**2 )
#Sécuriser le courant θ
old_sita_lmbd = Sita_lmbd[hh]
#Échantillonnage du nouveau θ (échantillonnage à pied ivre)
new_sita_lmbd = Sita_lmbd[hh] + ss.norm.rvs(loc=0, scale=sigma_lmbd, size=1)[0]
#Calcul de la vraisemblance (ajusté par Jacobien pour la vraisemblance logarithmique)
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
#Étape 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-Multivarié
### step4--------------------------------------
## '''Calcul côté lmbd'''
#Distribution postérieure de θ à l'étape 4 Initialisation des variables utilisées lors du calcul de la somme du premier terme du noyau
Lsita = np.zeros(nhh,dtype=np.float)
#Calcul d'erreur de la valeur d'utilité(Calcul de probabilité de θ)
Lsita = [np.sum( (u[hh,:] - np.sum(Ztld[:,:,hh]*Xs[hh,:,:], axis = 0) )**2 ) for hh in xrange(nhh)]
#Sécuriser le courant θ
old_sita_lmbd = Sita_lmbd[:]
#Échantillonnage du nouveau θ (échantillonnage à pied ivre)
new_sita_lmbd = Sita_lmbd[:] + ss.norm.rvs(loc=0, scale=sigma_lmbd, size=nhh)
#Calcul de la vraisemblance (ajusté par Jacobien pour la vraisemblance logarithmique)
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
#Étape MH
alpha = np.array(new_Lsita_lmbd/old_Lsita_lmbd)
alpha = alpha*(alpha<1) #S'il est supérieur à 1, 0 est entré ici.S'il est inférieur à 1, la valeur reste la même.
alpha[alpha==0] = 1 #1 est supérieur à 1,Courez ici
alpha[alpha!=alpha] = -1 #La partie qui était NaN-Définir sur 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) #Conditions d'adoption de MH
rej += 1*(tmp[0,:]<=0) #Si le résultat du calcul est le même que la valeur initiale, ajoutez 1 au rejet
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 est 1e-S'il est 6 ou plus, la valeur INV est adoptée, et si elle est inférieure à 6, la molécule devient trop petite (avertissement apparaît), donc 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]Est 1e-0 lorsque 6 ou moins
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"
##Boucle d'échantillonnage
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--------------------------------------
##Calcul côté DLT----
#Calcul des paramètres pour la distribution normale multivariée
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 )
#Échantillonnage avec distribution normale multivariée
Hdlt = np.ndarray.flatten( hbm.randn_multivariate(rtld, np.matrix(sig), n=nvar) )
##-----------------
##Calcul côté lmbd----
#Calcul des paramètres pour la distribution normale multivariée
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 )
#Échantillonnage avec distribution normale multivariée
Hlmbd = np.ndarray.flatten( hbm.randn_multivariate(rtld, np.matrix(sig), n=nvar) )
##-----------------
#--------------------------------------------
### step6--------------------------------------
##Calcul côté DLT
#Calcul des paramètres de la distribution inverse de Wishart
#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) #Dans le calcul ci-dessus, le div devient un vecteur horizontal, donc T est en retard pour faire de S un scalaire.
#Échantillonnage avec distribution de whishert inverse
Vsita_dlt = hbm.invwishartrand(f0 + nhh, V0 + S)
##------------
##Calcul côté lmbd
#Calcul des paramètres de la distribution inverse de Wishart
div = np.array(Sita_lmbd) - np.dot( D, Hlmbd.T ).T
S = np.dot(div, div.T)
#Échantillonnage avec distribution de whishert inverse
Vsita_lmbd = hbm.invwishartrand(f0 + nhh, V0 + S)
##------------
#--------------------------------------------
## step7-------------------------------------------
##Montant des dépenses Z mis à jour
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):
#Stocker dans l'ordre
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])]
#l'écriture
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")
Vient ensuite le fichier de fonction pour les branches et les feuilles.
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
###Définition des fonctions------------------------------------
#Estimation bayésienne du modèle probit binaire(Pour l'étape 1)
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 =Valeur Q
for t in xrange(lngth)
]
)
result[result == -np.inf] = -100
result[result == np.inf] = 100
return result[:,0]
#Une fonction qui standardise la valeur de Z(Pour l'étape 1)
def standardization(z):
return np.log( 0.001 + (z - np.min(z))/(np.max(z) - np.min(z))*(0.999 - 0.001) )
#Calcul de la partie multiplicatrice de la partie noyau de la distribution normale univariée(Pour l'étape 4)
def Nkernel(sita, H, D, Vsita):
if Vsita == 0: Vsita=1e-6
return ((sita - np.dot(H,D.T))**2)/Vsita
#Calcul de la partie multiplicatrice de la partie noyau de la distribution normale multivariée(Pour l'étape 4)
def NMkernel(sita, H, D, Vsita):
res = sita - np.dot(H.T, D)
return np.dot( np.dot(res.T, inverse(Vsita)), res )
###Définition de la matrice de la représentation de l'espace d'états du modèle de désaisonnalisation
def FGHset(al, k, p, q, nz): #al est le vecteur α du modèle AR, k;p;q est la tendance, la saison, la RA
m = k + p + q + nz -1
if(q>0): G = np.zeros((m,3+nz)) #Lorsque le modèle d'état comprend trois, tendance, saison et RA
else: G = np.zeros((m,2+nz)) #Lorsqu'il ne contient pas de composant AR(q=0)
F = np.zeros((m,m))
#H = matrix(0,1,m) #Ztld est utilisé à la place de H, donc H n'est pas nécessaire
##Construction de la matrice de blocs du modèle de tendance
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
##Construction d'une matrice par blocs de composantes de désaisonnalisation
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
##Construction d'une matrice de blocs de composants 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
##Calcul de la trame de la matrice de co-distribution de dispersion Q du modèle système
Q = np.identity(NS+nz)
return {'m':m, 'MatF':F, 'MatG':G, 'MatQ':Q}
#Définition de la matrice Q dans la représentation de l'espace d'états------------------------
def Qset(Q0,parm):
NS = len(Q0)
Q = Q0
#Calcul de la trame de la matrice de co-distribution de dispersion Q du modèle système
for i in xrange(NS): Q[i,i] = parm[i]
return np.array(Q)
#Fonction de filtre de Kalman------------------------------------------
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 trimestre à venir
XP = np.ndarray.flatten( np.dot(F, XF.T) ) #Puisqu'il devient un vecteur vertical à partir de la deuxième semaine, il est toujours converti en vecteur horizontal
VP = np.dot( np.dot(F, VF), F.T ) + np.dot( np.dot(G, Q), G.T)
#filtre
#R est un vecteur vertical s'il n'est pas manipulé. Notez que python est un vecteur horizontal!
if y[n] < limy:
NSUM = NSUM + 1
B = np.dot( np.dot(H[:,n], VP), H[:,n].T) + R #H est mathématiquement un vecteur horizontal
B1 = inverse(B) #vecteur vertical de dimension nvar
K = np.matrix(np.dot(VP, H[:,n].T)).T * np.matrix(B1) #K devient un vecteur vertical(matrix)
e = np.array(y[n]).T - np.dot(H[:,n], XP.T) #vecteur vertical de dimension nvar
XF = np.array(XP) + np.array( K * np.matrix(e) ).T #Vecteur horizontal
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] #Faites-en une matrice pour qu'elle puisse être calculée même dans une dimension
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}
#Fonction de lissage----------------------------------------------------
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}
#Définition de la fonction de vraisemblance logarithmique de 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 #Puisque l'optimeze est une fonction de minimisation, elle renvoie la vraisemblance logarithmique multipliée par moins.
#Définition de la fonction de génération de la distribution normale multivariée--------------------------------------
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
#Définition de la fonction de whishert inverse----------------------------------------------
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)))
# -------------------------------------------------------
#Calcul du montant d'inventaire des pseudo-ménages----------------------------------------------
def calc_INV(TIMES, PM, C, nhh):
##Données de tendance de consommation client
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)
# -------------------------------------------------------
Par souci de simplicité, nous avons réduit le nombre de variables explicatives lors du calcul des redevances. Veuillez noter que le contenu des livres est différent à cet égard. En outre, bien qu'il n'y ait pas de description spécifique dans le livre, le jugement de convergence de MCMC est également inclus. Le jugement est basé sur le diagnostic de Geweke. Pour un aperçu de MCMC, reportez-vous à iAnalysis-Dad's Analysis Diary-About MCMC (Markov Chain Monte Carlo). J'aimerais l'avoir.
Si vous souhaitez que les données soient utilisées dans le test, vous pouvez les obtenir à partir de la Page de support utilisateur du livre. Je pense que c'est facile et bon.
Je me suis abstenu de publier la sortie parce que c'était ennuyeux </ del>.
C'est de la négligence. Je suis désolé. .. ..
Si vous êtes intéressé, veuillez l'utiliser dans la pratique. (Je pense que mon implémentation est assez mauvaise, mais je pense qu'il vaut mieux l'écrire dans le langage du compilateur car c'est lent et pratiquement difficile en python ^^;)
Recommended Posts