Calcul de la fidélité des clients dans les séries chronologiques

Aperçu

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.

code

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.

Résumé

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

Calcul de la fidélité des clients dans les séries chronologiques
Différenciation des données de séries chronologiques (discrètes)
Analyse des séries chronologiques 3 Prétraitement des données des séries chronologiques
[Python] Accélère le chargement du fichier CSV de séries chronologiques
Calcul de la moyenne avant (temps de sommeil Piyolog)
Analyse des séries chronologiques 4 Construction du modèle SARIMA
mesure du temps
Décomposition des séries temporelles
Acquisition de données chronologiques (quotidiennes) des cours des actions
Lissage des séries temporelles et des données de forme d'onde 3 méthodes (lissage)
Voir les détails des données de séries chronologiques dans Remotte
Python: analyse des séries chronologiques
Détection d'anomalies des données de séries chronologiques par LSTM (Keras)
Question sur la série chronologique Python
Analyse des séries chronologiques RNN_LSTM1
Analyse des séries chronologiques 1 Principes de base
Mesure du temps d'exécution
Afficher les séries chronologiques TOPIX
Diagramme de séries chronologiques / Matplotlib
Reformatez l'axe des temps du graphique de la série chronologique des pandas avec matplotlib
Une histoire de regroupement de données de séries chronologiques d'échange
Mesure du temps de calcul à l'aide de maf
Calcul de similitude par MinHash
Opérateur de calcul de temps datetime timedelta
Analyse des séries chronologiques partie 4 VAR
Analyse de séries chronologiques Partie 3 Prévisions
À propos du calcul des coûts de MeCab
[Python] Tracer des données de séries chronologiques
Analyse de séries chronologiques Partie 1 Autocorrélation
Comment extraire des fonctionnalités de données de séries chronologiques avec les bases de PySpark
J'ai comparé le temps de calcul de la moyenne mobile écrite en Python
Comparaison de la prédiction des données de séries chronologiques entre le modèle SARIMA et le modèle Prophet