Wie der Titel schon sagt, habe ich die Kundenbindung anhand der folgenden Bücher implementiert.
Es kann jedoch einige seltsame Teile geben, da es sich um eine Kombination von Dingen handelt, die in der Vergangenheit hergestellt wurden ... (In meiner Umgebung hat es gut funktioniert.) Oh? Wenn es etwas gibt, das überraschend erscheint, kommentieren Sie bitte.
Ich werde den Code unten schreiben. Bitte beachten Sie, dass es eine kleine Menge gibt.
Erstens ist der Code der Datei auf der Hauptseite.
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
#Pfaddefinition für Datenimportdateipfade
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'
###Dateneingabe-------------------
#Erfassung von D-Daten (Daten, die zur Schätzung einzelner Parameter unter Verwendung von Gemeinsamkeiten erforderlich sind)
D = genfromtxt(inpt_D_filename, delimiter=',')
#Erfassung von C-Daten (tägliche Durchschnittsverbrauchsmenge)
C = genfromtxt(inpt_C_filename, delimiter=',')
#Erfassung von PM-Daten (Betrag nach Kaufdatum, um INV zu erhalten)
PM = genfromtxt(inpt_payment_filename, delimiter=',')
#y Daten abrufen (Normalverteilung nach Besuchsflag,Finden Sie b)
y = genfromtxt(inpt_visit_filename, delimiter=',')
print "done input"
### ------------------------------
#Vorkommen von zufälligen Samen
random.seed(555)
##--Definition der Konstante-------
TIMES = len(PM[0,:])
nhh = len(PM[:,0])
SEASONALYTY = 7 #Zeitraum saisonaler Schwankungen
RP = 50000 #Anzahl der von MCMC entnommenen Proben
keep = 10 #Datenintervall anzeigen
nz = 2 #Anzahl der erklärenden Variablen in der Formel für das Dienstprogramm für Ladenbesuche(Z_Anzahl von i)
nD = len(D[0,:]) #Vektorlänge von D.
m = 1 + 1 + nz
nvar = 1 #Anzahl der Objekte, die nach persönlichem Attribut berechnet werden sollen (2 für frische Lebensmittel und andere Sonntagsprodukte)
limy = 1e20 #Grenze der als fehlend anzusehenden Zahlenwerte
k = 1 #Reihenfolge des Trendkomponentenmodells
zc = k+1+ SEASONALYTY-2 + nz
pr = 4 #Anzahl der Prozesse
##-------------------
###Datengenerierungsprozess-----------------------------------
##Protokoll der Anzahl der Tage seit dem letzten Besuch,Vorhandensein / Nichtvorhandensein eines Ereignisses (Dummy-Variable)------
##Mit oder ohne Rabatt (Dummy-Variable)
Ztld=np.zeros((zc,TIMES,nhh))
## ----------------------
#Ztld Formung
Ztld[0,:,:] = [[1]*nhh]*TIMES
Ztld[1,:,:] = [[1]*nhh]*TIMES
Ztld[zc-2,:,:] = genfromtxt(inpt_Z_filename, delimiter=',').T #Besuchsintervall vor Ort
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"
###------------------------------------------------------------
##Vorherige Verteilungsparameter-----------------
# '''Schritt 1: Für latente Utility-Sampling'''
A = 0.01 * np.identity(zc) ##A ist B._Inverse von 0
b0 = np.zeros((zc,nvar),dtype=np.float)
#Schritt 3: Zur verteilten Abtastung von Systemrauschen
mu0 = 0; kaps0 = 25
nu0 = 0.02; s0 = 0.02
Sita_sys0 = np.array([np.float(10)]*m)
#Schritt 5: Datenrahmen des Regressionsparameters H der Verbraucherheterogenität
m0 = np.zeros((nD,nvar),dtype=np.float)
A0 = 0.01 * np.identity(nD) #In Higuchimoto, A.
#Schritt 6: Datenrahmen des Regressionsparameters V der Verbraucherheterogenität
f0 = nvar+3
V0 = f0 * np.identity(nvar)
##------------------------------------
##Erstellen Sie einen Datenrahmen, der für die Abtastung der posterioren Verteilung erforderlich ist-----------
#Schritt 1: Datenrahmen für latent Utility Sampling
u = np.zeros((nhh,TIMES),dtype=np.float)
# step2:Datenrahmen zur Berechnung des Zustandsvektors
#Zur Vereinfachung der Verarbeitung im Anfangswert-Einstellungsteil einstellen
#Schritt 3: Datenrahmen für die verteilte Abtastung von Systemrauschen
Sita_sys = np.array([np.float(10)*np.identity(m) for i in xrange(nhh)]) #Da es sich um 3D handelt, kann es nicht in eine Matrix konvertiert werden, sodass es zum Zeitpunkt der Berechnung in eine Matrix konvertiert wird.
#Schritt 4: Datenrahmen für die Parameterabtastung zur Angabe des Pseudo-Haushaltsinventars
#Wenn θ 2 Variablen hat, machen Sie es Matrix anstelle von Vektor
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)
##---------------------------------------------------------
##Anfangswerteinstellung------------------------------
#Für Schritt 1
Xs = np.zeros((nhh, zc, TIMES),dtype=np.float) #Mann,Anzahl der Variablen,time
sigma = 1.0
##Für Schritt 2
param = hbm.FGHset(0, 1, SEASONALYTY, 0, nz)
L = 1
R = np.identity(L)
F = np.array(param['MatF'])
G = np.array(param['MatG'])
#Ein Rahmen zum Speichern der Verteilung des Systemmodells für jede Person
Q0 =np.array(param['MatQ'])
#Für Schritt 3
mu = 0.
sigs = 1.
##-------------------------------------------
##Angeben des Schnittbereichs
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 Probenahme(Keine Schleifen, sondern individuelle Berechnungen)
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)):
##Berechnung der Systemmodellparameter in Schritt 2----------------------
#Numerische Berechnung, um die wahrscheinlichste Schätzung von TAU2 zu finden------------------------------
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)
#Höchstwahrscheinlich Schätzung von TAU2
#Kalman Filter
#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
#Glätten----------------------------------------------------------
LLF3 = hbm.SMO(XPS, XFS, VPS, VFS, F, GSIG2, 1, SEASONALYTY, 1, zc, TIMES)
Xs[hh,:,:] = np.array(LLF3['XSS']).T #Konvertieren Sie zwangsweise, um dem Typ zu entsprechen
return Xs[hh,:,:]
#------------------------------------------------------------
def calculate_difference((hh, Xs)):
#Berechnung der Differenz in Schritt 3--------------------------------------
#Initialisierung von Variablen, die bei der Berechnung der Summe der Differenzberechnung in Schritt 3 verwendet werden
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): #Kompatibel mit allen Menschen
#Berechnung der Differenz in Schritt 3--------------------------------------
#Initialisierung von Variablen, die bei der Berechnung der Summe der Differenzberechnung in Schritt 3 verwendet werden
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)): #Θ von Schritt 4_Berechnung von δ
### step4--------------------------------------
## '''Berechnung auf der dlt-Seite'''
#Posteriore Verteilung von θ in Schritt 4 Initialisierung der Variablen, die bei der Berechnung der Summe des ersten Terms des Kernels verwendet werden
Lsita = 0.
#Fehlerberechnung des Nutzwerts(Wahrscheinlichkeitsberechnung von θ)
Lsita = sum( (u[hh,:] - np.diag(np.dot(Ztld[:,:,hh].T, Xs[hh,:,:])))**2 )
#Sichern Sie den Strom θ
old_sita_dlt = Sita_dlt[hh]
#Sampling new θ (Drunken Walk Sampling)
new_sita_dlt = Sita_dlt[hh] + ss.norm.rvs(loc=0, scale=sigma_dlt,size=1)[0]
#Berechnung der Wahrscheinlichkeit (angepasst von Jacobian für die logarithmische Wahrscheinlichkeit)
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 Schritt
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-Multivariate Version
### step4--------------------------------------
## '''Berechnung auf der dlt-Seite'''
#Posteriore Verteilung von θ in Schritt 4 Initialisierung der Variablen, die bei der Berechnung der Summe des ersten Terms des Kernels verwendet werden
Lsita = np.zeros(nhh,dtype=np.float)
#Fehlerberechnung des Nutzwerts(Wahrscheinlichkeitsberechnung von θ)
Lsita = [np.sum( (u[hh,:] - np.sum(Ztld[:,:,hh]*Xs[hh,:,:], axis = 0) )**2 ) for hh in xrange(nhh)]
#Sichern Sie den Strom θ
old_sita_dlt = Sita_dlt[:]
#Sampling new θ (Drunken Walk Sampling)
new_sita_dlt = Sita_dlt[:] + ss.norm.rvs(loc=0, scale=sigma_dlt, size=nhh)
#Berechnung der Wahrscheinlichkeit (angepasst von Jacobian für die logarithmische Wahrscheinlichkeit)
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 Schritt
alpha = np.array(new_Lsita_dlt/old_Lsita_dlt)
alpha = alpha*(alpha<1) #Wenn es größer als 1 ist, wird hier 0 eingegeben.Wenn es kleiner als 1 ist, bleibt der Wert gleich.
alpha[alpha==0] = 1 #1 ist größer als 1,Lauf hier
alpha[alpha!=alpha] = -1 #Der Teil, der NaN war-Auf 1 setzen
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 Adoptionsbedingungen
Sita_dlt = Sita_tmp*(Sita_tmp>=0) + old_sita_dlt*(Sita_tmp<0) #Die Bedingung, dass δ in dieser Analyse nicht negativ ist
tmp = Sita_dlt - old_sita_dlt #
rej += 1*(tmp[0,:]==0) #Wenn das Ergebnis der Berechnung mit dem Anfangswert übereinstimmt, addieren Sie 1 zur Zurückweisung
#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)): #Θ von Schritt 4_Berechnung von λ
### step4--------------------------------------
## '''Berechnung auf der lmbd Seite'''
#Posteriore Verteilung von θ in Schritt 4 Initialisierung der Variablen, die bei der Berechnung der Summe des ersten Terms des Kernels verwendet werden
Lsita = 0.
#Fehlerberechnung des Nutzwerts(Wahrscheinlichkeitsberechnung von θ)
Lsita = sum( (u[hh,:] - np.diag(np.dot(Ztld[:,:,hh].T, Xs[hh,:,:])))**2 )
#Sichern Sie den Strom θ
old_sita_lmbd = Sita_lmbd[hh]
#Sampling new θ (Drunken Walk Sampling)
new_sita_lmbd = Sita_lmbd[hh] + ss.norm.rvs(loc=0, scale=sigma_lmbd, size=1)[0]
#Berechnung der Wahrscheinlichkeit (angepasst von Jacobian für die logarithmische Wahrscheinlichkeit)
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 Schritt
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-Multivariate
### step4--------------------------------------
## '''Berechnung auf der lmbd Seite'''
#Posteriore Verteilung von θ in Schritt 4 Initialisierung der Variablen, die bei der Berechnung der Summe des ersten Terms des Kernels verwendet werden
Lsita = np.zeros(nhh,dtype=np.float)
#Fehlerberechnung des Nutzwerts(Wahrscheinlichkeitsberechnung von θ)
Lsita = [np.sum( (u[hh,:] - np.sum(Ztld[:,:,hh]*Xs[hh,:,:], axis = 0) )**2 ) for hh in xrange(nhh)]
#Sichern Sie den Strom θ
old_sita_lmbd = Sita_lmbd[:]
#Sampling new θ (Drunken Walk Sampling)
new_sita_lmbd = Sita_lmbd[:] + ss.norm.rvs(loc=0, scale=sigma_lmbd, size=nhh)
#Berechnung der Wahrscheinlichkeit (angepasst von Jacobian für die logarithmische Wahrscheinlichkeit)
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 Schritt
alpha = np.array(new_Lsita_lmbd/old_Lsita_lmbd)
alpha = alpha*(alpha<1) #Wenn es größer als 1 ist, wird hier 0 eingegeben.Wenn es kleiner als 1 ist, bleibt der Wert gleich.
alpha[alpha==0] = 1 #1 ist größer als 1,Lauf hier
alpha[alpha!=alpha] = -1 #Der Teil, der NaN war-Auf 1 setzen
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 Adoptionsbedingungen
rej += 1*(tmp[0,:]<=0) #Wenn das Ergebnis der Berechnung mit dem Anfangswert übereinstimmt, addieren Sie 1 zur Zurückweisung
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 ist 1e-Wenn es 6 oder mehr ist, wird der INV-Wert übernommen, und wenn er kleiner als 6 ist, wird das Molekül zu klein (Warnung erscheint), also 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]Ist 1e-0 wenn 6 oder weniger
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"
##Abtastschleife
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--------------------------------------
##Berechnung auf der dlt-Seite----
#Berechnung von Parametern für die multivariate Normalverteilung
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 )
#Probenahme mit multivariater Normalverteilung
Hdlt = np.ndarray.flatten( hbm.randn_multivariate(rtld, np.matrix(sig), n=nvar) )
##-----------------
##Berechnung auf der lmbd Seite----
#Berechnung von Parametern für die multivariate Normalverteilung
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 )
#Probenahme mit multivariater Normalverteilung
Hlmbd = np.ndarray.flatten( hbm.randn_multivariate(rtld, np.matrix(sig), n=nvar) )
##-----------------
#--------------------------------------------
### step6--------------------------------------
##Berechnung auf der dlt-Seite
#Berechnung von Parametern der inversen Wishart-Verteilung
#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) #In der obigen Berechnung wird das div zu einem horizontalen Vektor, so dass T dahinter liegt, um S zu einem Skalar zu machen.
#Probenahme mit inverser Whishert-Verteilung
Vsita_dlt = hbm.invwishartrand(f0 + nhh, V0 + S)
##------------
##Berechnung auf der lmbd Seite
#Berechnung von Parametern der inversen Wishart-Verteilung
div = np.array(Sita_lmbd) - np.dot( D, Hlmbd.T ).T
S = np.dot(div, div.T)
#Probenahme mit inverser Whishert-Verteilung
Vsita_lmbd = hbm.invwishartrand(f0 + nhh, V0 + S)
##------------
#--------------------------------------------
## step7-------------------------------------------
##Aktualisierter Z-Ausgabenbetrag
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):
#In Ordnung lagern
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])]
#Schreiben
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")
Als nächstes folgt die Funktionsdatei für die Zweige und Blätter.
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
###Funktionsdefinition------------------------------------
#Bayesianische Schätzung des binären Probit-Modells(Für Schritt 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 =Q-Wert
for t in xrange(lngth)
]
)
result[result == -np.inf] = -100
result[result == np.inf] = 100
return result[:,0]
#Eine Funktion, die den Wert von Z standardisiert(Für Schritt 1)
def standardization(z):
return np.log( 0.001 + (z - np.min(z))/(np.max(z) - np.min(z))*(0.999 - 0.001) )
#Berechnung des Multiplikatorteils des Kernelteils der univariaten Normalverteilung(Für Schritt 4)
def Nkernel(sita, H, D, Vsita):
if Vsita == 0: Vsita=1e-6
return ((sita - np.dot(H,D.T))**2)/Vsita
#Berechnung des Multiplikatorteils des Kernelteils der multivariaten Normalverteilung(Für Schritt 4)
def NMkernel(sita, H, D, Vsita):
res = sita - np.dot(H.T, D)
return np.dot( np.dot(res.T, inverse(Vsita)), res )
###Matrixeinstellung der Zustandsraumdarstellung des saisonalen Anpassungsmodells
def FGHset(al, k, p, q, nz): #al ist der α-Vektor des AR-Modells, k;p;q ist Trend, Saison, AR
m = k + p + q + nz -1
if(q>0): G = np.zeros((m,3+nz)) #Wenn das Zustandsmodell drei enthält, Trend, Saison und AR
else: G = np.zeros((m,2+nz)) #Wenn keine AR-Komponente enthalten ist(q=0)
F = np.zeros((m,m))
#H = matrix(0,1,m) #Ztld wird anstelle von H verwendet, daher ist H nicht erforderlich
##Konstruktion der Blockmatrix des Trendmodells
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
##Konstruktion einer Blockmatrix von saisonalen Anpassungskomponenten
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
##Konstruktion der Z-Komponentenblockmatrix
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
##Berechnung des Rahmens der Dispersions-Co-Verteilungsmatrix Q des Systemmodells
Q = np.identity(NS+nz)
return {'m':m, 'MatF':F, 'MatG':G, 'MatQ':Q}
#Einstellen der Matrix Q in der Zustandsraumdarstellung------------------------
def Qset(Q0,parm):
NS = len(Q0)
Q = Q0
#Berechnung des Rahmens der Dispersions-Co-Verteilungsmatrix Q des Systemmodells
for i in xrange(NS): Q[i,i] = parm[i]
return np.array(Q)
#Kalman Filterfunktion------------------------------------------
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 Laufzeit voraus Prognose
XP = np.ndarray.flatten( np.dot(F, XF.T) ) #Da es ab der zweiten Woche ein vertikaler Vektor wird, wird es immer in einen horizontalen Vektor konvertiert
VP = np.dot( np.dot(F, VF), F.T ) + np.dot( np.dot(G, Q), G.T)
#Filter
#R ist ein vertikaler Vektor, wenn er nicht manipuliert wird. Beachten Sie, dass Python ein horizontaler Vektor ist!
if y[n] < limy:
NSUM = NSUM + 1
B = np.dot( np.dot(H[:,n], VP), H[:,n].T) + R #H ist mathematisch ein horizontaler Vektor
B1 = inverse(B) #Vertikaler Vektor der nvar-Dimension
K = np.matrix(np.dot(VP, H[:,n].T)).T * np.matrix(B1) #K wird ein vertikaler Vektor(matrix)
e = np.array(y[n]).T - np.dot(H[:,n], XP.T) #Vertikaler Vektor der nvar-Dimension
XF = np.array(XP) + np.array( K * np.matrix(e) ).T #Horizontaler Vektor
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] #Machen Sie es Matrix, so dass es auch in einer Dimension berechnet werden kann
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}
#Glättungsfunktion----------------------------------------------------
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}
#Definition der logarithmischen Wahrscheinlichkeitsfunktion von 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 #Da optimeze eine Minimierungsfunktion ist, wird die logarithmische Wahrscheinlichkeit multipliziert mit minus zurückgegeben.
#Definition der Generierungsfunktion der multivariaten Normalverteilung--------------------------------------
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
#Definition der inversen Whishert-Funktion----------------------------------------------
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)))
# -------------------------------------------------------
#Berechnung des Pseudohaushaltsinventarbetrags----------------------------------------------
def calc_INV(TIMES, PM, C, nhh):
##Kundenverbrauchstrend Daten
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)
# -------------------------------------------------------
Der Einfachheit halber haben wir die Anzahl der erklärenden Variablen bei der Berechnung der Lizenzgebühren reduziert. Bitte beachten Sie, dass der Inhalt der Bücher in dieser Hinsicht unterschiedlich ist. Obwohl das Buch keine spezifische Beschreibung enthielt, ist auch das Konvergenzurteil von MCMC enthalten. Das Urteil basiert auf Gewekes Diagnose. Eine Übersicht über MCMC finden Sie unter iAnalysis-Dads Analysetagebuch über MCMC (Markov-Kette Monte Carlo). Ich hätte es gerne.
Wenn Sie möchten, dass die Daten für den Test verwendet werden, können Sie sie von der Benutzer-Support-Seite des Buches herunterladen. Ich finde es einfach und gut.
Ich habe es unterlassen, die Ausgabe zu veröffentlichen, weil sie ärgerlich war </ del>.
Es ist Nachlässigkeit. Es tut mir Leid. .. ..
Wenn Sie interessiert sind, verwenden Sie es bitte in der Praxis. (Ich denke, meine Implementierung ist auch ziemlich schlecht, aber ich denke, es ist besser, sie in der Compilersprache zu schreiben, weil sie in Python langsam und praktisch schwierig ist ^^;)
Recommended Posts