流失預測模型實證-Pareto/NBD模型

       互聯網購物基本是一種非契約型協議，顧客的購買行為均具有隨機性和不可預測性，那如何在此激烈的網絡市場立於不敗之地，

那就應該盡可能的降低網絡顧客的流失率。

目前用於預測顧客流失率的模型有：

      SVM模型，Logistics模型，Pareto/NBD模型，BG/NBD模型以及引申的各類模型，通過結合分類模型評估方法，就可以檢驗模型的准確率，

從而進一步應用與實際用例中。

其中：Pareto/NBD模型介紹：

 

Pareto/NBD模型：

一：變量：

ID：用戶ID

x:在指定時間段內，重復購買次數(不包括第一次購買)

Tx:最后一次購買的時長，即在指定時間段內，第一次購買到最后一次購買的周期數，保留小數

T:第一購買到最后指定時間的最后一天的周期數，保留小數

T_star: 預測時間周期長度

 

 

二：4個參數的估計：

 

R code:

 

 

 

Python Code:

 

import pandas as pd
from pnbd import pnbd  # 自己開發的模塊

consume_sources=pd.read_excel('test.xls')
cal_cbs=consume_sources.loc[:,['x','t_x','T_cal']]

fun,x=pnbd.EstimateParameters(cal_cbs,par_start = [1, 1, 1, 1], max_param_value = 1000)
print fun,x
# 0.4846256 1.4887628 0.6738844 4.4011076
# 150346.098305 [ 0.48458758  1.48866284  0.67405456  4.40356461]

 

兩個軟件運行的結果一樣

 

 

三：預測用戶的活躍度
params=x
consume_sources['P']=consume_sources.apply(lambda x:
     pnbd.PAlive(params,x=x['x'], t_x=x['t_x'], T_cal=x['T_cal']),axis=1)

 

print consume_sources.head(20)

 

 

*** 閾值選取可以選擇算法選擇合理的閾值

 

活躍度計算指標：

樣本活躍百分比=3817/22731=16.79%

Accuracy准確率=(2421+15342)/22731=78.14%

Sensitivity覆蓋率=2421/3817=63.43%

Precision命中率=2421/5993=40.40%

 

流失率計算指標：

樣本流失百分比=3817/22731=83.21%

Accuracy准確率=(2421+15342)/22731=78.14%

Sensitivity覆蓋率=15342/18914=81.12%

Precision命中率=15342/16738=91.66%

 

四：預測未來一個月用戶的購買期望次數：

consume_sources['E']=consume_sources.apply(lambda x:
                 pnbd.ConditionalExpectedTransactions(params,T_star=5,
                    x=x['x'], t_x=x['t_x'], T_cal=x['T_cal']),axis=1)
consume_sources.to_excel('test.xls',index=False)

 

Python代碼：

pnbd.py

# -*-coding:utf-8-*-

'''
Pareto/NBD模型代碼：
仿照R語言：BTYD
params <- c("r", "alpha", "s", "beta")
# 分組求和
pnbd.compress.cbs(cbs, 0)
# 個體客戶在時刻T的活躍度模型
pnbd.PAlive(params, 10, 35, 39)
# 客戶在時間T后t*時間內發生交易的交易次數的期望模型
pnbd.ConditionalExpectedTransactions(params, T.star=2, x=10, t.x=35, T.cal=39)
# 計算對數極大似然值
pnbd.LL(params, cbs[, "x"], cbs[, "t.x"], cbs[, "T.cal"])
pnbd.cbs.LL(params, cbs)
# 評估4個參數，采用最小方差法
result<-pnbd.EstimateParameters(cal.cbs, par.start = c(1, 1, 1, 1),max.param.value = 100)
'''


from math import lgamma
from scipy.special import gammaln
import math
from scipy.optimize import minimize
import numpy as np


def LL(params, x, t_x, T_cal):
    x = np.matrix(x).T
    t_x = np.matrix(t_x).T
    T_cal = np.matrix(T_cal).T

    # z:為數值
    def h2f1(a, b, c, z):
        lenz = len(z)
        j = 0
        uj =[1 for i in range(1,lenz+1)]
        uj=np.matrix(uj).T
        y = uj
        lteps = 0

        while (lteps < lenz):
            lasty = y
            j = j + 1
            uj = np.multiply(uj , (a + j - 1))
            uj = np.multiply(uj,(b + j - 1))
            uj = np.multiply(uj,z)
            temp=np.multiply((c + j - 1) ,j)
            uj = np.multiply(uj,1/temp)
            y = y + uj
            lteps = sum(y == lasty)

        return y

    r = params[0]
    alpha = params[1]
    s = params[2]
    beta = params[3]
    maxab = max(alpha, beta)
    absab = abs(alpha - beta)
    param2 = s + 1

    if (alpha < beta):
        param2 = r + x

    part1 = r * np.log(alpha) + s * np.log(beta) - lgamma(r) + gammaln(r +x)
    part2 = np.multiply(-(r + x),np.log(alpha + T_cal)) - s * np.log(beta + T_cal)

    if (absab == 0):
        partF = np.multiply(-(r + s + x),np.log(maxab + t_x))\
              + np.log(1 - np.power(((maxab + t_x) / (maxab + T_cal)),r + s + x))

    else:
        F1 = h2f1(r + s + x, param2, r + s + x + 1, absab / (maxab +t_x))
        F2 = h2f1(r + s + x, param2, r + s + x + 1,absab / (maxab + T_cal))
        F2 = np.multiply(F2,np.power((maxab + t_x)/(maxab + T_cal), r + s + x))

        partF = np.multiply(-(r + s + x), np.log(maxab + t_x)) + np.log(F1 - F2)

    part3 = np.log(s) - np.log(r + s + x) + partF
    return (part1 + np.log(np.exp(part2) + np.exp(part3)))


# cbs:DataFrame , params:list
def cbs_LL(params, cal_cbs):

    x = cal_cbs['x'].values
    t_x = cal_cbs['t_x'].values
    T_cal = cal_cbs['T_cal'].values

    if  'custs' in cal_cbs.columns:
        custs = cal_cbs['custs'].values
    else:
        custs = [1 for i in range(len(cal_cbs))]

    custs = np.matrix(custs).T
    ll=sum(np.multiply(custs, LL(params,x,t_x,T_cal)))

    return ll[0,0]


# 個體客戶在時刻T的活躍度模型
# params=[0.55, 10.56, 0.61, 11.64]  PAlive(params, 10, 35, 39)
def PAlive(params, x, t_x, T_cal):
    # h2f1： F(a,b,c,z)為高斯超幾何函數
    def h2f1(a, b, c, z):
        j = 0
        uj = 1
        y = uj
        lteps = 0
        while (lteps < 1):
            lasty = y
            j = j + 1
            uj = uj * (a + j - 1) * (b + j - 1) / (c + j - 1) * z / j
            y = y + uj
            if y == lasty:
                lteps = lteps + 1
        return y

    r = params[0]
    alpha = params[1]
    s = params[2]
    beta = params[3]

    if (alpha > beta):
        F1 = h2f1(r + s + x, s + 1, r + s + x + 1, (alpha - beta) / (alpha + t_x))
        F2 = h2f1(r + s + x, s + 1, r + s + x + 1, (alpha - beta) / (alpha + T_cal))
        # 與R語言中BTYD.pnbd.PAlive有差異
        A0 = F1 / math.pow(alpha+t_x,r+x+s) - F2 / math.pow(alpha + T_cal,r + s + x)
        # A0 = F1 / (math.pow(alpha+T_cal,s)*math.pow(alpha + t_x,r+x)) - \
        #      F2 / math.pow(alpha + T_cal,r + s + x)

    elif (alpha < beta):
        F1 = h2f1(r + s + x, r + x, r + s + x + 1, (beta - alpha) / (beta + t_x))
        F2 = h2f1(r + s + x, r + x, r + s + x + 1, (beta - alpha) / (beta + T_cal))
        A0 = F1 / math.pow(beta + t_x,r + s + x) - \
             F2 / math.pow(beta + T_cal , r + s + x)
    else:
        return math.pow(1 + s / (r + s + x) * (math.pow((alpha+T_cal)/(alpha+t_x),r+x+s)-1),-1)

    return math.pow(1 + s / (r + s + x) * math.pow(alpha + T_cal,r + x) * math.pow(beta +T_cal,s)* A0,-1)



# 估計4個參數: c("r", "alpha", "s", "beta")
# 所以參數不超過：max_param_value值
# for example:
# EstimateParameters(data_Soures, par_start = [1, 1, 1, 1], max_param_value =100)
def EstimateParameters(cal_cbs, par_start = [1, 1, 1, 1], max_param_value = 10000):
    # 對數似然值
    def pnbd_ell(params,cal_cbs,max_param_value):

        params = np.exp(params)
        params[0] = min(params[0], max_param_value)
        params[1] = min(params[1], max_param_value)
        params[2] = min(params[2], max_param_value)
        params[3] = min(params[3], max_param_value)
        return -1 * cbs_LL(params=params,cal_cbs=cal_cbs)

    logparams = np.log(par_start)

    results = minimize(fun=pnbd_ell,x0=logparams,
                       args=(cal_cbs,max_param_value),
                       method="L-BFGS-B")

    return results['fun'],np.exp(results['x'])


# 單一值
# 客戶在時間T后t*時間內發生交易的交易次數的期望模型
def ConditionalExpectedTransactions(params, T_star, x, t_x, T_cal):
    r = params[0]
    alpha = params[1]
    s = params[2]
    beta = params[3]
    P1 = (r + x) * (beta + T_cal) / ((alpha + T_cal) * (s - 1))
    # 與R語言中BTYD.pnbd.ConditionalExpectedTransactions有差異
    # P2 < - (1 - ((beta + T.cal) / (beta + T.cal + T.star)) ^ (s - 1))
    P2 = (1 -math.pow(((beta + T_cal) / (beta + T_cal + T_star)) , (s - 1)) )
    P3 = PAlive(params, x, t_x, T_cal)

    return P1 * P2 * P3