简述回归分析及其在车险定价中的应用

广义线性模型目前仍是财产保险，尤其是车险定价的主流模型，而大部分从网络上、书本上能查到的有关资料似乎都在通往实际应用的道路上给初学者留下了一些坑，比如为什么要关注被解释变量服从什么分布、连结函数的选择与被解释变量所属分布是否有关、为什么说最小二乘法是训练广义线性模型的一种特例、如何处理离散且无序的解释变量、为何要用最大似然估计而非矩估计等，在此一并做个系统阐述。

一、回归分析概述

（一）什么是回归分析

回归的英文是regression，可以理解成从现象“退回”本质，其实就是找规律。规律一般可以表述为“如果…那么…”，在数学上能给出这样表述的东西称为模型。其实就是一个多输入多输出的函数，“如果”后面的内容是入参，一般也称为解释变量或自变量；“那么”后面的内容是出参，一般也称为被解释变量或因变量。当出参可以用有序数值表示时，这个模型就称为回归模型，否则就称为聚类模型，本文只阐述回归模型。

（二）主流回归模型有哪些

当回归模型的入参与出参呈现单调关系时，大多数都可以考虑广义线性模型，否则就只能考虑非线性模型，比如幂多项式模型，或人工神经网络（一种模拟生物神经系统，理论上可以拟合任意非线性关系的黑箱模型，其基本原理可以参看《详解反向传播神经网络》）。本文重点阐述目前在财产险定价上被普遍使用的广义线性模型，主要是因为其适用性强，不需要消耗太多算力，且具有解释性。

（三）基本概念说明

为了方便后面的阐述，首先对几种不同的变量做区分说明如下：

1. 【普通变量】不具随机性的变量。用小写字母表示，如$x$。
2. 【随机变量】带有随机性的变量。用大写字母表示，如$X$。
3. 【标量】单一数量。用不加粗字母表示，如$x$或$X$。
4. 【向量】多维数量。用加粗字母表示，如$\mathbf{x}$或$\mathbf{X}$。
5. 【估计量】对某一个实际量做估计的样本函数，当样本被观测前，它是随机变量；被观测后，它成为普通变量。用小写字母上加“^”表示，如$\hat{\theta }$。

再说一下什么是样本。一对入参与出参即称为一个样本，已观测到的作为训练样本，未观测到的作为未知样本，回归模型就是由训练样本得到，用来由未知样本入参估计其出参的函数。对于一个线性模型，假设由$m$个训练样本，其所有入参记为$$\mathbf{x}=\left[\begin{array}{cccc}{\mathbf{x}}_1& {\mathbf{x}}_2& \cdots & {\mathbf{x}}_m\end{array}\right]=\left[\begin{array}{cccc}1& 1& \cdots & 1\\ x_{11}& x_{21}& \cdots & x_{m1}\\ x_{12}& x_{22}& \cdots & x_{m2}\\ ⋮& ⋮& ⋮& ⋮\\ x_{1n}& x_{2n}& \cdots & x_{mn}\end{array}\right]$$对应这$m$个样本入参，就有$m$个样本出参，本文仅讨论一元出参的情况，即所有出参都是标量。这些样本出参在未观测前都是随机变量，被观测后就是普通变量，由于训练样本都是已观测到的，故这$m$个训练样本的出参记为$$\mathbf{y}=\left[\begin{array}{c}{y}_1\\ y_2\\ ⋮\\ y_m\end{array}\right]$$这里面有一个很关键的概念要理解，就是每个出参在未被观测前都可视为一个随机变量，其入参的不同造成了该随机变量所属分布的不同，所以针对每个入参要预测最可能看到的出参结果，其实是要先得到该出参的具体分布，进而可以得到该出参的期望，也就是最有可能得到的观测结果。而所谓回归模型就是对每个入参${\mathbf{x}}_i$可以给出其出参期望的无偏估计的函数。之所以又是期望又是无偏估计说得这么绕口，就是因为任何观测结果本身就带有随机性，而回归模型中的参数本身就是估计量，也带有随机性，对于两个随机变量，只能说他们的期望相同。即$$E(Y|{\mathbf{x}}_i)=h_{\mathbf{\theta }}({\mathbf{x}}_i)=E[h_{\hat{\mathbf{\theta }}}({\mathbf{x}}_i)]$$其中$\mathbf{\theta }$为回归函数中的参数向量，简称回归参数。如果上面的概念觉得还是太绕，可以简单理解回归函数$h_{\mathbf{\theta }}()$的作用就是对于每个入参${\mathbf{x}}_i$给出最有可能观测到的出参结果$y_i$。如普通线性回归函数可表示为：$$h_{\hat{\mathbf{\theta }}}({\mathbf{x}}_i)={\mathbf{x}}_i^T\hat{\mathbf{\theta }}$$其中$$\hat{\mathbf{\theta }}=\left[\begin{array}{c}{\hat{\theta }}_0\\ {\hat{\theta }}_1\\ ⋮\\ {\hat{\theta }}_n\end{array}\right]$$
所谓做回归就是找到最能得到已知样本观测结果的一组回归参数$\hat{\mathbf{\theta }}$，回归参数之所以是估计量，是因为它其实是训练样本的函数，而训练样本是因为被观测了才从随机变量变成了普通变量。做回归的过程一般也称为模型训练，这是因为除普通线性回归模型外，大部分回归模型的参数是很难用样本函数解析式直接表达出来的，只能通过迭代算法找到数值解，这个过程就像是在做反复训练，至于如何训练，后面会讲到。

估计量的随机性来自于训练样本的观测误差；对于某一入参，其出参的随机性来自于对该出参的观测误差。

二、广义线性模型（Generalized Linear Models, GLM）

（一）由普通线性模型向广义线性模型的推广

首先回顾一下普通线性模型（Normal Linear Models, NLM）。它简单的给出了入参和出参的一个线性关系。但它其实是广义线性模型的一个特例，其作为广义线性模型时的训练方法恰好就是求使得各残差（residual，出参观测值与回归估计值之间的差异）平方和取得最小值时的回归参数的过程，因此叫最小二乘法（Least squares），英文直译过来其实叫最小平方更合适。使用NLM的条件可简述为：$$\{\begin{array}{l}{Y}_i=Y|{\mathbf{x}}_i\sim N({\mu }_i,{\sigma }^2)\\ {\mu }_i={\mathbf{x}}_i^T\mathbf{\theta }\end{array}$$即对于各入参${\mathbf{x}}_i$，对应的各出参${\mathbf{Y}}_i$相互独立，且服从同方差正态分布；同时对应的各期望与各入参线性相关。这个两个条件在实际应用中其实都很苛刻，首先有些变量本来就不可能是正态分布的，比如出险次数，只可能是正整数；还比如赔款金额，也只可能是正数，而正态分布则是要求该变量在整个实数域里都有取值可能，只是概率不同罢了。其次是要求出参的期望直接就等于入参的线性函数，就更不用说其局限性了。因此需要把普通线性模型的使用条件做一些拓展，从而得到广义线性模型，其使用条件可简述为：$$\{\begin{array}{l}{Y}_i=Y|{\mathbf{x}}_i\sim 均值为{\mu }_i的指数族分布\\ g({\mu }_i)={\mathbf{x}}_i^T\mathbf{\theta }\end{array}$$即对于各入参${\mathbf{x}}_i$，对应的各出参${\mathbf{Y}}_i$相互独立，且只需要服从指数族分布即可；与各入参线性相关的不再必须是各分布的期望了，而是期望的连结函数$g()$即可，这个连结函数（link function）如果等于其对应出参所属指数族分布的自然参数时，则称为正则连结函数（canonical link function），这时会使计算似然函数梯度时简便很多，这点后面在模型训练中会看到。但具体选择什么连结函数一般是通过观察训练样本散点图来决定的，只要该函数单调可微即可。由于连结函数的存在，GLM的回归函数就变成了：$$h_{\mathbf{\theta }}({\mathbf{x}}_i)=g^{-1}({\mathbf{x}}_i^T\mathbf{\theta })$$

之所以要考虑出参服从什么分布，是为了最大限度的利用已知信息。比如当已知出参只会取到正整数（即出参服从泊松分布）时，就不应再考虑取到非正整数的可能，从而影响出参的期望估计。

（二）指数族分布（Exponential Family Distribution）

前面提到GLM需要出参服从指数族分布，那什么是指数族分布，为何说正态分布、泊松分布、伽马分布等一众常见分布都属于指数族分布，指数族分布的期望和方差是否有计算通式，这点关系到做迭代训练时的梯度通式推导。此处将回答这些问题。

若随机变量的概率质量函数（Probability Mass Function, PMF）或概率密度函数（Probability Density Function, PDF）可以化为以下形式：$$f_Y(y)=\exp \left[\frac{\eta y-b(\eta )}{\phi }+c(y,\phi )\right]$$则称此随机变量服从指数族分布，记为$Y\sim EF(\eta ,\phi )$。

其中$\eta$为该分布的自然参数或正则参数，是与分布期望有函数关系的唯一参数；$\phi$称为离散参数，仅与分布方差有关；$b$和$c$是两个已知函数，$b(\eta )$决定了分布的具体形式，$c(y,\phi )$起到了归一化作用（即使得该PDF的全域积分或该PMF的全域和为1）。可以证明，指数族分布的期望和方差分别为：$$\{\begin{array}{l}E(Y)=b^{\prime }(\eta )\\ V(Y)=\phi b^{\prime \prime }(\eta )\end{array}$$可见指数族分布的方差可以看作是其期望的函数。以上这些我们只需要记住即可，至于为何是这样，且留给数学专业的同学去证明好了，就像会用电脑不一定需要会造电脑一样。下面我们可以用几个常见分布印证一下上述结论。

1、正态分布（Normal Distribution）

可将其PDF化成指数族分布通式：$$\begin{array}{rl}{f}_Y(y)& =\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp \left[-\frac{(y-\mu {)}^2}{2{\sigma }^2}\right]\\ & =\exp \left[\frac{-y^2+2y\mu -{\mu }^2}{2{\sigma }^2}+\ln (\frac{1}{\sqrt{2\pi {\sigma }^2}})\right]\\ & =\exp \left[\frac{\mu y-{\mu }^2/2}{{\sigma }^2}+\ln (\frac{1}{\sqrt{2\pi {\sigma }^2}})-\frac{y^2}{2{\sigma }^2}\right]\end{array}$$可见：$$\eta =\mu$$$$\phi ={\sigma }^2$$$$b(\eta )=\frac{{\eta }^2}{2}$$$$c(y,\phi )=\ln (\frac{1}{\sqrt{2\pi {\sigma }^2}})-\frac{y^2}{2{\sigma }^2}=\ln (\frac{1}{\sqrt{2\pi \phi }})-\frac{y^2}{2\phi }$$于是有：$$\{\begin{array}{l}E(Y)=b^{\prime }(\eta )=(\frac{{\eta }^2}{2}{)}^{\prime }=\eta =\mu \\ V(Y)=\phi b^{\prime \prime }(\eta )={\sigma }^2(\frac{{\eta }^2}{2}{)}^{\prime \prime }={\sigma }^2\end{array}$$

2、泊松分布（Poisson Distribution）

可将其PMF化成指数族分布通式：$$\begin{array}{rl}{f}_Y(y)& =\frac{e^{-\lambda }{\lambda }^y}{y!}\\ & =\exp \left(\ln \frac{e^{-\lambda }{\lambda }^y}{y!}\right)\\ & =\exp \left[y\ln (\lambda )-\lambda -\ln (y!)\right]\end{array}$$可见：$$\eta =\ln (\lambda )$$$$\phi =1$$$$b(\eta )=\lambda =e^{\eta }$$$$c(y,\phi )=-\ln (y!)$$于是有：$$\{\begin{array}{l}E(Y)=b^{\prime }(\eta )=(e^{\eta }{)}^{\prime }=e^{\eta }=\lambda \\ V(Y)=\phi b^{\prime \prime }(\eta )=(e^{\eta }{)}^{\prime \prime }=e^{\eta }=\lambda =\mu \end{array}$$

3、伽马分布（Gamma Distribution）

可将其PDF化成指数族分布通式： f Y ( y ) = β α y α − 1 e − β y Γ ( α ) = exp ⁡ { − β y + α ln ⁡ ( β ) + ( α − 1 ) ln ⁡ ( y ) − ln ⁡ [ Γ ( α ) ] } = exp ⁡ { β α y + ln ⁡ ( β α ) − 1 α + α ln ⁡ ( α ) + ( α − 1 ) ln ⁡ ( y ) − ln ⁡ [ Γ ( α ) ] } \begin{aligned} f_Y(y) &= {\beta^\alpha y^{\alpha-1} e^{-\beta y} \over \Gamma(\alpha)} \\ &= \exp \left \{ -\beta y + \alpha \ln(\beta) +(\alpha-1) \ln (y) - \ln[\Gamma(\alpha)] \right \} \\ &= \exp \left \{ {{\beta \over \alpha}y+\ln({\beta \over \alpha}) \over -{1 \over \alpha}}+ \alpha \ln(\alpha) +(\alpha-1) \ln (y) - \ln[\Gamma(\alpha)] \right \} \end{aligned} fY​(y)​=Γ(α)βαyα−1e−βy​=exp{−βy+αln(β)+(α−1)ln(y)−ln[Γ(α)]}=exp{−α1​αβ​y+ln(αβ​)​+αln(α)+(α−1)ln(y)−ln[Γ(α)]}​可见：$$\eta =\frac{\beta }{\alpha }$$$$\phi =-\frac{1}{\alpha }$$$$b(\eta )=-\ln (\frac{\alpha }{\beta })=\ln (\eta )$$$$c(y,\phi )=\alpha \ln (\alpha )+(\alpha -1)\ln (y)-\ln [\Gamma (\alpha )]=\frac{1}{\phi }\ln (-\phi )-(\frac{1+\phi }{\phi })\ln (y)-\ln [\Gamma (-\frac{1}{\phi })]$$于是有：$$\{\begin{array}{l}E(Y)=b^{\prime }(\eta )=\frac{1}{\eta }=\frac{\alpha }{\beta }\\ V(Y)=\phi b^{\prime \prime }(\eta )=-\frac{1}{\alpha }(-\frac{1}{{\eta }^2})=\frac{\alpha }{{\beta }^2}=\frac{{\mu }^2}{\alpha }\end{array}$$

（三）模型训练

回归分析的目的是找到一个回归参数向量$\hat{\mathbf{\theta }}$ ，使得对任一入参回归出的出参期望（${\mu }_i=g^{-1}({\mathbf{x}}_i^T\hat{\mathbf{\theta }})$ ）所对应的随机变量分布观测到当前样本出参的概率整体最大，即最大似然得到当前样本观测值。为衡量这个似然度，可定义似然函数为：$$L(\hat{\mathbf{\theta }})=\prod _{i=1}^m{f}_{Y|{\mathbf{x}}_i}(y_i;\hat{\mathbf{\theta }})$$

只有通过最大似然估计才能把包括出参所属分布类型在内的所有已知信息都充分利用。而矩估计（即寻找使得分布期望等于样本均值的一组回归参数）仅适用于样本量足够大的时候。而且对于入参都是分类变量时，使用矩估计得到的模型其实是一个训练样本统计表，运行模型的过程就是查表，这种模型其实根本没有预测性，所以除非样本量已经大到几乎覆盖全体才没有必要再去做预测了。而似然估计是用当前观测样本来估计出参最有可能服从的分布，而该分布的期望不一定就是当前观测样本的均值。

使上述似然函数取得最大值的$\hat{\mathbf{\theta }}$ 同样也可以使该似然函数的对数取得最大值，因此为简便计算，可将指数族分布的似然函数写为：$$L(\hat{\mathbf{\theta }})=\sum _{i=1}^m{l}_i(\hat{\mathbf{\theta }})$$根据指数族分布PDF通式，其中，$$l_i(\hat{\mathbf{\theta }})=\frac{{\eta }_i{y}_i-b({\eta }_i)}{\phi }+c(y_i,\phi )$$当把各训练样本出入参代入此似然函数后，如何找到使该似然函数取得最大值的$\hat{\mathbf{\theta }}$ ？其实也简单，只需要在$n$维空间（设$\hat{\mathbf{\theta }}$向量有$n$个维度）中先随便找到一个点，再求得该点的似然函数梯度（一个指向函数值变化最大方向，模长是该方向函数值变化率的向量），然后沿此梯度方向前进一小步到下一个位置，再计算这个点的梯度，以此类推，最后找到一个梯度足够小的位置，即得到了要找的目标$\hat{\mathbf{\theta }}$ 。这就是所谓的梯度下降法（虽然是沿着梯度上升的方向移动，但梯度大小是在下降的，所以叫梯度下降）。当然，如果可以直接求解出梯度为0点，就不用这样迭代了。但无论哪种方式，写出似然函数梯度表达式都是绕不开的。

前人已证明，指数族分布的似然函数曲面在$n$维空间只存在单峰。

由GLM前提假设，可令：$$g({\mu }_i)=s_i={\mathbf{x}}_i^T\hat{\mathbf{\theta }}$$又由前述指数族分布期望计算公式$${\mu }_i=b^{\prime }({\eta }_i)$$可得该似然函数梯度的一个分量：$$\frac{\partial L}{\partial {\hat{\theta }}_j}=\sum _{i=1}^m\left(\frac{\partial l_i}{\partial {\eta }_i}\frac{\partial {\eta }_i}{\partial {\mu }_i}\frac{\partial {\mu }_i}{\partial s_i}\frac{\partial s_i}{\partial {\hat{\theta }}_j}\right)$$其中，$$\frac{\partial l_i}{\partial {\eta }_i}=\frac{y_i-b^{\prime }({\eta }_i)}{\phi }=\frac{y_i-{\mu }_i}{\phi }$$$$\frac{\partial {\eta }_i}{\partial {\mu }_i}={\left(\frac{\partial {\mu }_i}{\partial {\eta }_i}\right)}^{-1}=\frac{1}{b^{\prime \prime }({\eta }_i)}$$$$\frac{\partial {\mu }_i}{\partial s_i}={\left(\frac{\partial s_i}{\partial {\mu }_i}\right)}^{-1}=\frac{1}{g^{\prime }({\mu }_i)}$$$$\frac{\partial s_i}{\partial {\hat{\theta }}_j}=x_{ij}$$代入整理可得似然函数完整梯度通式：$$\nabla L=\left[\begin{array}{c}\frac{\partial L}{\partial {\hat{\theta }}_0}\\ \frac{\partial L}{\partial {\hat{\theta }}_1}\\ ⋮\\ \frac{\partial L}{\partial {\hat{\theta }}_n}\end{array}\right]=\left[\begin{array}{c}{\sum }_{i=1}^m\left[\frac{(y_i-{\mu }_i)}{{\sigma }_i^2{g}^{\prime }({\mu }_i)}\right]\\ {\sum }_{i=1}^m\left[\frac{(y_i-{\mu }_i)x_{i1}}{{\sigma }_i^2{g}^{\prime }({\mu }_i)}\right]\\ ⋮\\ {\sum }_{i=1}^m\left[\frac{(y_i-{\mu }_i)x_{in}}{{\sigma }_i^2{g}^{\prime }({\mu }_i)}\right]\end{array}\right]=\mathbf{x}(\mathbf{y}-\mathbf{\mu })\ast \frac{1}{{\mathbf{\sigma }}^2}\ast \frac{1}{g^{\prime }(\mathbf{\mu })}$$其中，“$\ast$”在这里代表两向量对应位置元素相乘；$$\mathbf{x}=\left[\begin{array}{cccc}{\mathbf{x}}_1& {\mathbf{x}}_2& \cdots & {\mathbf{x}}_m\end{array}\right]=\left[\begin{array}{cccc}1& 1& \cdots & 1\\ x_{11}& x_{21}& \cdots & x_{m1}\\ x_{12}& x_{22}& \cdots & x_{m2}\\ ⋮& ⋮& ⋮& ⋮\\ x_{1n}& x_{2n}& \cdots & x_{mn}\end{array}\right]，\mathbf{y}=\left[\begin{array}{c}{y}_1\\ y_2\\ ⋮\\ y_m\end{array}\right]，\hat{\mathbf{\theta }}=\left[\begin{array}{c}{\hat{\theta }}_0\\ {\hat{\theta }}_1\\ ⋮\\ {\hat{\theta }}_n\end{array}\right]$$$$\mathbf{\mu }=\left[\begin{array}{c}{\mu }_1\\ {\mu }_2\\ ⋮\\ {\mu }_m\end{array}\right]=\left[\begin{array}{c}{g}^{-1}({\mathbf{x}}_1^T\hat{\mathbf{\theta }})\\ g^{-1}({\mathbf{x}}_2^T\hat{\mathbf{\theta }})\\ ⋮\\ g^{-1}({\mathbf{x}}_m^T\hat{\mathbf{\theta }})\end{array}\right]=g^{-1}({\mathbf{x}}^T\hat{\mathbf{\theta }})，\frac{1}{g^{\prime }(\mathbf{\mu })}=\left[\begin{array}{c}\frac{1}{g^{\prime }({\mu }_1)}\\ \frac{1}{g^{\prime }({\mu }_2)}\\ ⋮\\ \frac{1}{g^{\prime }({\mu }_m)}\end{array}\right]，\frac{1}{{\mathbf{\sigma }}^2}=\left[\begin{array}{c}\frac{1}{{\sigma }_1^2}\\ \frac{1}{{\sigma }_2^2}\\ ⋮\\ \frac{1}{{\sigma }_m^2}\end{array}\right]=\left[\begin{array}{c}\frac{1}{v({\mu }_1)}\\ \frac{1}{v({\mu }_2)}\\ ⋮\\ \frac{1}{v({\mu }_m)}\end{array}\right]$$其中，${\sigma }_i^2$为第$i$个入参对应的分布方差，由于指数族函数的方差可视为其期望的函数，因此可写成${\sigma }_i^2=v({\mu }_i)$。可见，要计算似然函数梯度，除知道各训练样本出入参和连结函数形式外，还需要知道包含回归参数的各出参所属分布的期望和方差表达式。

当连结函数取正则时有：$$g({\mu }_i)={\eta }_i=s_i={\mathbf{x}}_i^T\hat{\mathbf{\theta }}$$此时似然函数梯度分量为：$$\frac{\partial L}{\partial {\hat{\theta }}_j}=\sum _{i=1}^m\left(\frac{\partial l_i}{\partial {\eta }_i}\frac{\partial {\eta }_i}{\partial {\mu }_i}\frac{\partial {\mu }_i}{\partial s_i}\frac{\partial s_i}{\partial {\hat{\theta }}_j}\right)=\sum _{i=1}^m\left(\frac{\partial l_i}{\partial {\eta }_i}\frac{\partial {\eta }_i}{\partial {\mu }_i}\frac{\partial {\mu }_i}{\partial {\eta }_i}\frac{\partial {\eta }_i}{\partial {\hat{\theta }}_j}\right)=\sum _{i=1}^m\left(\frac{\partial l_i}{\partial {\eta }_i}\frac{\partial {\eta }_i}{\partial {\hat{\theta }}_j}\right)=\sum _{i=1}^m\left[\frac{(y_i-{\mu }_i)x_{ij}}{\phi }\right]$$因此似然函数梯度为：$$\nabla L=\left[\begin{array}{c}\frac{\partial L}{\partial {\hat{\theta }}_0}\\ \frac{\partial L}{\partial {\hat{\theta }}_1}\\ ⋮\\ \frac{\partial L}{\partial {\hat{\theta }}_n}\end{array}\right]=\frac{1}{\phi }\left[\begin{array}{c}{\sum }_{i=1}^m(y_i-{\mu }_i)\\ {\sum }_{i=1}^m(y_i-{\mu }_i)x_{i1}\\ ⋮\\ {\sum }_{i=1}^m(y_i-{\mu }_i)x_{in}\end{array}\right]=\frac{1}{\phi }\mathbf{x}(\mathbf{y}-\mathbf{\mu })$$NLM相当于出参服从同方差正态分布，连结函数取期望原值的GLM。而巧合的是，正态分布PDF写成指数族分布通式时的自然参数就是其期望，因此取期望原值的连结函数也恰好就是其正则连结函数，即$$\mathbf{\mu }={\mathbf{x}}^T\hat{\mathbf{\theta }}$$代入上式，并直接令梯度为0得到估计方程：$$\nabla L=\frac{1}{\phi }\mathbf{x}(\mathbf{y}-{\mathbf{x}}^T\hat{\mathbf{\theta }})=0$$通过矩阵运算求解可得：$$\hat{\mathbf{\theta }}=(\mathbf{x}{\mathbf{x}}^T{)}^{-1}\mathbf{x}\mathbf{y}$$上式恰好与最小二乘法的计算公式一样。

除了NLM的情况外，就很难直接从估计方程得到回归参数的解析表达式了。下面分别分析一下出参服从泊松分布和伽马分布两种情况的回归过程。假设连结函数都选择成对数函数，即$$g({\mu }_i)=\ln ({\mu }_i)={\mathbf{x}}_i^T\hat{\mathbf{\theta }}$$则有：$${\mu }_i=e^{{\mathbf{x}}_i^T\hat{\mathbf{\theta }}}$$$$g^{\prime }({\mu }_i)=\frac{1}{{\mu }_i}=\frac{1}{e^{{\mathbf{x}}_i^T\hat{\mathbf{\theta }}}}$$

1、出参服从泊松分布

由前述指数族分布论证可知泊松分布的方差表达式为：$${\sigma }_i^2=e^{{\eta }_i}={\mu }_i=e^{{\mathbf{x}}_i^T\hat{\mathbf{\theta }}}$$代入似然函数梯度各分量$$\frac{\partial L}{\partial {\hat{\theta }}_j}=\sum _{i=1}^m\left[\frac{(y_i-{\mu }_i)x_{ij}}{{\sigma }_i^2{g}^{\prime }({\mu }_i)}\right]$$由于泊松分布化成指数族分布通式后，其自然参数$\eta =\ln (\mu )$，可见对数函数恰好是其正则连结函数，这就使得上式中的${\sigma }_i^2$和$g^{\prime }({\mu }_i)$恰好抵消掉了，这就是为何前面说取正则连结函数时会大大简化梯度计算。此时似然函数梯度为：$$\nabla L=\left[\begin{array}{c}\frac{\partial L}{\partial {\hat{\theta }}_0}\\ \frac{\partial L}{\partial {\hat{\theta }}_1}\\ ⋮\\ \frac{\partial L}{\partial {\hat{\theta }}_n}\end{array}\right]=\left[\begin{array}{c}{\sum }_{i=1}^m(y_i-e^{{\mathbf{x}}_i^T\hat{\mathbf{\theta }}})\\ {\sum }_{i=1}^m(y_i-e^{{\mathbf{x}}_i^T\hat{\mathbf{\theta }}})x_{i1}\\ ⋮\\ {\sum }_{i=1}^m(y_i-e^{{\mathbf{x}}_i^T\hat{\mathbf{\theta }}})x_{in}\end{array}\right]$$显然难以求解出此梯度0点处的回归参数解析式。此时就只能用前面提到的梯度下降法来迭代出使似然函数取得最大值的回归参数数值解。

2、出参服从伽马分布

由前述指数族分布论证可知伽马分布的方差表达式为：$${\sigma }_i^2=-\frac{1}{\alpha }\left(-\frac{1}{{\eta }_i^2}\right)=\frac{{\mu }_i^2}{\alpha }=\frac{e^{2{\mathbf{x}}_i^T\hat{\mathbf{\theta }}}}{\alpha }$$代入似然函数梯度各分量$$\frac{\partial L}{\partial {\hat{\theta }}_j}=\sum _{i=1}^m\left[\frac{(y_i-{\mu }_i)x_{ij}}{{\sigma }_i^2{g}^{\prime }({\mu }_i)}\right]$$由于对数函数不是伽马分布的正则连结函数，所以此时的似然函数梯度比较复杂：$$\nabla L=\left[\begin{array}{c}\frac{\partial L}{\partial {\hat{\theta }}_0}\\ \frac{\partial L}{\partial {\hat{\theta }}_1}\\ ⋮\\ \frac{\partial L}{\partial {\hat{\theta }}_n}\end{array}\right]=\left[\begin{array}{c}{\sum }_{i=1}^m\frac{(y_i-e^{{\mathbf{x}}_i^T\hat{\mathbf{\theta }}})}{e^{{\mathbf{x}}_i^T\hat{\mathbf{\theta }}}}\\ {\sum }_{i=1}^m\frac{(y_i-e^{{\mathbf{x}}_i^T\hat{\mathbf{\theta }}})}{e^{{\mathbf{x}}_i^T\hat{\mathbf{\theta }}}}{x}_{i1}\\ ⋮\\ {\sum }_{i=1}^m\frac{(y_i-e^{{\mathbf{x}}_i^T\hat{\mathbf{\theta }}})}{e^{{\mathbf{x}}_i^T\hat{\mathbf{\theta }}}}{x}_{in}\end{array}\right]$$此梯度的零点也只能用梯度下降法迭代出数值解。具体怎么迭代，在后面的示例中会给出具体代码。

（四）模型检验

实践是检验真理的唯一标准。所有回归模型建立后，也都需要做检验，来定量的判断这个模型的好坏。许多文献都是通过计算决定系数来检验出参的绝大部分离差是否都被回归掉了，以此来判断模型有效性，但它其实检验的是入参对样本做分类的效果好坏，是否同类样本间的出参差异足够小，异类样本间的出参差异足够大。而我们真正需要检验的是模型预测出的结果是否与实际观测足够接近。此处给出一个定量检验模型预测效果的通用方法，即拟合优度的卡方检验。

首先要有一批检验样本，由于任何同分布随机变量的和都可以近似视为正态分布，所以对于每个检验样本的入参${\mathbf{x}}_i$所对应的所有出参采样$Y_{i1}\sim Y_{im}$，设${\hat{y}}_i=h_{\hat{\theta }}({\mathbf{x}}_i)$，如果这个回归模型效果很理想，那么此时出参的统计均值应该满足：$${\bar{Y}}_i\sim N({\hat{y}}_i,\frac{{\sigma }^2}{m})$$设$S_i^2$为入参${\mathbf{x}}_i$所对应的所有出参采样的样本方差，即：$$S_i^2=\frac{1}{m-1}\sum _{j=1}^m(Y_{ij}-{\bar{Y}}_i{)}^2$$于是根据卡方分布定义可得：$$\frac{(m-1)S_i^2}{{\sigma }^2}=\sum _{j=1}^m{\left(\frac{Y_{ij}-{\bar{Y}}_i}{\sigma }\right)}^2\sim {\chi }^2(m-1)$$所以根据$t$分布定义可得：$$\frac{{\bar{Y}}_i-{\hat{y}}_i}{\sigma /\sqrt{m}}/\sqrt{\frac{(m-1)S_i^2}{{\sigma }^2(m-1)}}=\frac{{\bar{Y}}_i-{\hat{y}}_i}{S_i/\sqrt{m}}\sim t(m-1)\approx N(0,1)$$当$m$足够大时，$t$分布可以近似为标准正态分布，而$n$个独立服从标准正态分布的随机变量的平方和又服从卡方分布，即对于$n$个检验样本入参，其对应的$n$个估计出参和$n$个出参观测均值应该满足$$\sum _{i=1}^n{\left(\frac{{\bar{Y}}_i-{\hat{y}}_i}{S_i/\sqrt{m}}\right)}^2\sim {\chi }^2(n)$$故可通过将各入参处的出参采样观测结果带入上式，看是否落入所设定的置信区间内。如落入$(0,{\chi }_{1-\alpha }^2]$，则可以$1-\alpha$的置信水平接受此回归模型（一般取$\alpha =5\%$）。

三、对分类变量的处理

在许多实际应用（如车险定价）中经常会发现入参是分类变量的情况，如地区、性别、车系，这种变量取值的特点是不连续、没顺序，但数量有限。而回归分析都是在寻找数值到数值的关系，而对于分类变量，给各类别赋个有序数值是没有意义的，而且更重要的是也没必要给出一个连续回归函数，因为类别是枚举的，不存在渐变的中间值。此时有意义的是针对每个类别给出对应被解释变量的分布。仍要以线性函数的形式实现这个目的的方法就是对此变量做独热（One-Hot）编码。即将一个有$n$个枚举值的分类变量用$n-1$个0-1变量表示。

（一）独热（One-Hot）编码

设有分类变量“地区”，其取值包括{北京，天津，上海}这三种。对其做独热编码的过程如下图所示：

首先“地区”这一个变量的三个取值变成三个取值只有0和1的三个“伪变量”，此时“北京”、“天津”、“上海”三个取值就可以被分别编码成“1，0，0”、“0，1，0”、“0，0，1”，特点是每一行和每一列都只能有一个1，即“独热”。但由于设定了“地区”的取值一定会是这三个地方中的一个，而排除掉所有不选的，就表示选择了最后剩下的，所以一个有$n$个取值的分类变量其实只需要$n-1$个伪变量即可完成独热编码。做独热编码的目的是方便以线性回归的形式表示分类估计，在此例中要实现的分类估计效果是：$$\{\begin{array}{l}E(Y|地区=北京)=g^{-1}({\mathbf{x}}_{北京}^T\mathbf{\theta })=g^{-1}({\theta }_0+{\theta }_{北京}\ast 1+{\theta }_{天津}\ast 0)=g^{-1}({\theta }_0+{\theta }_{北京})\\ E(Y|地区=天津)=g^{-1}({\mathbf{x}}_{天津}^T\mathbf{\theta })=g^{-1}({\theta }_0+{\theta }_{北京}\ast 0+{\theta }_{天津}\ast 1)=g^{-1}({\theta }_0+{\theta }_{天津})\\ E(Y|地区=上海)=g^{-1}({\mathbf{x}}_{上海}^T\mathbf{\theta })=g^{-1}({\theta }_0+{\theta }_{北京}\ast 0+{\theta }_{天津}\ast 0)=g^{-1}({\theta }_0)\end{array}$$可见，此时对独热伪变量的线性回归，就是对原分类解释变量的各类别所对应的回归参数做估计（上例中${\theta }_{上海}$相当于默认为0）。

（二）连续变量离散化

对于多个分类变量，以及由多个分类变量与多个连续数值变量形成混合解释向量，仍适用此方法。对于其中的连续解释变量可以被分箱离散化，并在该变量的各分箱因子与变量值之前也形成了一个回归函数。如车险定价模型中的新车购置价，存在线性关系的其实是各价位上的回归因子与当前价位。有时为了提升模型对于异常数据的鲁棒性、方便分析、以及兼容非线性回归，一般对于连续型解释变量或是状态过多的离散变量，都需要做离散化处理，即分箱。

分箱的理想效果是使得不同箱入参对应的出参分布差异尽可能大，而同箱入参对应的出参分布差异尽可能小。为达到这种效果，在此介绍两种方法，卡方分箱和最佳KS值分箱。

1、卡方分箱

先按最小颗粒度对该解释变量做分箱排序，然后将具有最小卡方值的相邻箱体合并在一起，直到满足确定的停止准则。相邻箱体间的卡方值为：$${\chi }^2=\sum _{i=1}^2\sum _{j=1}^k\frac{(A_{ij}-E_{ij}{)}^2}{E_{ij}}$$其中，$k$为样本的分类数量；$A_{ij}$表示第$i$个箱体中，属于第$j$类的样本个数； $E_{ij}$表示第$i$个箱体中第$j$类样本按照平均占比应有的个数，即$$E_{ij}=\frac{C_j}{N}{R}_i$$其中，$N$为两箱体中的所有样本总数；$C_j$为两箱体中第$j$类样本总数；$R_i$为第$i$个箱体中的样本总数。可见，如果两个箱体间的样本分布足够均匀，则卡方值几乎为零。

上式卡方分布统计量的推导过程此处不做详述，原理就是因为某类样本的出现次数可视为服从泊松分布，而泊松分布的期望与方差相同，同时泊松分布也可以近似看作是正态分布，而正态分布减期望后除标准差就得到标准正态分布，多个标准正态分布的平方和就是卡方分布。

2、最优KS值分箱

也是先按最小颗粒度对该解释变量做分箱排序，然后找到KS值最大的那个箱体作为划分点将所有箱体分成两组，在各组中反复重复此动作，直到达到目标分组数。设目前共有$n$个分箱，则第$i$个箱体上$KS$值为正负样本的累计出现频率之差，到当前箱体的累计正样本数在正样本总数的占比即为正样本累计出现频率，当前箱体的累计负样本数在负样本总数的占比即为负样本累计出现频率，即：$$K{S}_i=|cumGoodRat{e}_i-cumBadRat{e}_i|$$$$\{\begin{array}{l}cumGoodRat{e}_i=\frac{{\sum }_{bins=1}^{i}good{s}_{bins}}{{\sum }_{bins=1}^{n}good{s}_{bins}}\\ cumBadRat{e}_i=\frac{{\sum }_{bins=1}^{i}bad{s}_{bins}}{{\sum }_{bins=1}^{n}bad{s}_{bins}}\end{array}$$直观看就是正负样本的累计概率曲线间的最大差值。具体示例见下表：

可见，选择5号箱体作为划分点的效果是可以最大限度（82.75%）将负样本分出去，以最小限度（29.69%）同时将正样本也分出去。当然我们最希望的是找到这么个分箱位置，所有正样本都在它的一边，而所有负样本都在它的另一边。当然这实际是不可能的，能达到上表中第5个分箱处53.06%的KS值就很不错了。

四、以车险纯风险保费厘定举例

我国自2015年启动商业车险费率改革以来，各保险公司的车险自主定价能力愈显重要，用到的关键技术就是本文所述的广义线性模型。建模原理并不难，难的是拿到足够的数据去做维度细分，从而实现更精细化定价，准确挖掘优质业务，使承保资源得到最有效配置，进而在激烈的市场竞争中赢得先机。为更直观的阐明前面叙述的建模原理，将通过一个实际建模来说明。在实际建模前，一般先要通过检验各维度的因子显著性，来挑选出有效的风险定价维度。检验方法详见《因子显著性的方差检验》。

由前文可知，对于分类变量，模型训练的结果是给每个变量的每个取值赋一个回归因子，如《车损险纯风险保费定价模型示例》所示。为了简便，此示例只选择了车龄、座位数、车主性别、车系、NCD（即无赔优待系数）这5个分类维度来为车损险（商业车险中的一个险别）风险成本（即纯风险保费）定价。在实际应用中要用到的定价维度会更多。要计算某一投保单中车损险风险成本时，根据该投保单上的车龄、座位数等信息，得到各维度上对应因子值做连乘（因为连结函数选择对数函数），再乘基准风险成本即得到当前单的车损险风险成本。同理可得当前单上其他险别的风险成本，加总即得当前单的整体风险成本。再结合当前公司的非市场费用率、可接受的经营成本率以及当前市场情况，就可以评估出可以给到当前单的保费折扣和市场费用率。

（一）构造数据

我们就用《车损险纯风险保费定价模型示例》中的5个维度分箱因子来伪造一批保单数据用于建模，并假设基准出险频度（即每车年出险次数）为20%，基准案均赔款为5000元。每个业务板块生成20张保单，根据示例表中的因子和基准值来分别生成各保单的出险次数和各次出险的赔款金额。根据前人经验，出险次数一般服从泊松分布，赔款金额服从伽马分布，具体数据构造代码如下：

import os
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import numpy.matlib as matlib
import scipy.stats as st

BaseCF = 0.2                              # 基准出险频度（100%）
BasePPC = 5000                            # 基准案均赔款（元）
NUM = 100                                 # 各板块保单数
a = 4                                     # 设定伽马分布的形状参数
mypath = os.getcwd() + "\定价模型"        # 定位数据目录

# =============================================================================
# 按照分箱因子表给定的各维度各分箱系数计算出各业务板块的出险频度和案均赔款
# =============================================================================
dimensionTable = pd.read_excel(mypath + "/" + "车险分箱因子.xlsx") # 分箱因子表
DTcolumns = dimensionTable.columns.values # 分箱因子表各列名

# 获取各分箱维度名列表
dimensionNames = []
i = 0
while i < len(DTcolumns):
    dimensionNames.append(DTcolumns[i])
    i += 4

# 获取各维度的分箱数
dimensionSize = []
for i in dimensionNames:
    dimensionSize.append(sum(dimensionTable[i].notna()))

# 计算各因子的各分箱在笛卡尔组合表中每次出现时的连续出现次数
rep = []    # 第i列数值表示第i维中各分箱每次出现要连续复制的个数
for i in range(1, len(dimensionSize)):
    j = i + 1
    k = dimensionSize[i]
    while j < len(dimensionSize):
        k = k * dimensionSize[j]
        j += 1
    rep.append(k)
rep.append(1)

# 计算各因子的各分箱在笛卡尔组合表中的循环次数
loo = []
n = rep[0]*dimensionSize[0]             # 笛卡尔表总行数
for i, j in zip(rep, dimensionSize):
    loo.append(n / (i*j))

# 构造一个由各维度各分箱构造的笛卡尔组合(即由各维度划分的所有业务板块)表
bct = []                                # 业务板块表 
for r, l, dn, ds in zip(rep, loo, dimensionNames, dimensionSize):
    tt = []
    while l > 0:
        for i in dimensionTable[dn][:ds]:
            tt.extend([i]*r)
        l -= 1
    bct.append(tt)
# 通过矩阵化将原二阶列表转置后写入DataFrame
bctDF = pd.DataFrame(np.matrix(bct).T.tolist(), columns=dimensionNames)

# 根据分箱因子表中的分值填写各业务板块出险频度和案均赔款
KPl = []                                # 初始化出险频度因子乘积列表
KGl = []                                # 初始化案均赔款因子乘积列表
for i in bctDF.index:
    KP = 1
    for j in range(len(bctDF.columns)):
        KP = KP * dimensionTable[(dimensionTable[bctDF.columns[j]]==bctDF.iloc[
                i, j])][bctDF.columns[j] + '系数P'].iloc[0]
    KPl.append(KP)
    KG = 1
    for j in range(len(bctDF.columns)):
        KG = KG * dimensionTable[(dimensionTable[bctDF.columns[j]]==bctDF.iloc[
                i, j])][bctDF.columns[j] + '系数G'].iloc[0]
    KGl.append(KG)
# 分别乘以基准值后写入板块期望表
bctDF['出险频度'] = np.array(KPl) * BaseCF
bctDF['安均赔款'] = np.array(KGl) * BasePPC

# =============================================================================
# 根据板块期望表分别按泊松分布和伽马分布生成各单理赔数据（出险次数和各次出险的赔款金额）
# =============================================================================
# 生成出险次数并填入案均赔款
policyBlocks = []                         # 初始化各保单所属板块的描述清单
policyClaimNums = []                      # 初始化各保单的出险次数清单
policyMeanClaimPays = []                  # 初始化各保单的案均赔款清单
for i in bctDF.index:
    policyBlocks.extend([bctDF.loc[i][:-2].tolist()] * NUM)  # 扩充板块描述各列
    claimFre = bctDF.loc[i][-2]                          # 当前板块出险频度
    claimNum = st.poisson.rvs(mu=claimFre, size=NUM)     # 当前板块各保单的出险次数
    policyClaimNums.extend(claimNum.tolist())            # 扩充出险次数列
    policyMeanClaimPays.extend([bctDF.loc[i][-1]] * NUM) # 扩充案均赔款列

# 生成各保单号
policyNos = [str(n) for n in (np.arange(0, len(policyBlocks))+100000000).tolist()]

# 写入DataFrame
claimLDF = pd.DataFrame(policyBlocks, columns=dimensionNames) # 板块描述各列
claimLDF['出险次数'] = policyClaimNums                         # 添加出险次数列
claimLDF.insert(0, '保单号', policyNos)                        # 添加保单号列
claimLDF['赔款金额'] = policyMeanClaimPays                     # 添加案均赔款列

# 如果当前行出险大于1次，则在当前行下复制一行并减1次出险次数
for i in claimLDF.index:
    if claimLDF['出险次数'][i] > 1:  # 如果当前行出险次数大于1
        # 在当前行下面复制一行
        claimLDF = claimLDF.iloc[:i].append(claimLDF.iloc[i],
                                ignore_index=True).append(claimLDF.iloc[i:],
                                                 ignore_index=True)
        # 将新增的这行出险次数减1
        claimLDF.loc[i+1, '出险次数'] -= 1

# 对出险次数不为0的行，按案均赔款生成伽马随机数填入赔款金额列；其余赔款金额置0
for i in claimLDF.index:
    if claimLDF['出险次数'][i] > 0:      # 如果当前行出险次数大于0
        # 按赔款金额列中的值（案均赔款）生成伽马随机数
        claimPay = st.gamma.rvs(a=a, scale=claimLDF['赔款金额'][i]/a, 
                                size=10).tolist()[0]
        claimLDF.loc[i, '赔款金额'] = claimPay
    else:                                 # 如果出险次数为0 
        claimLDF.loc[i, '赔款金额'] = 0    # 则赔款金额也为0 
        
# 将保单理赔样本数据保存成excel
claimLDF.to_excel(mypath + "/" + '车险保单理赔样本.xlsx', index=False)

（二）训练出险次数因子

使用伪造出的保单理赔清单训练出险次数因子。根据前人经验，此时连结函数一般选择对数函数。具体代码如下：

# =============================================================================
# 根据保单理赔样本训练出险次数回归因子
# =============================================================================
# 从理赔列表中去重读取各单最大出险次数
ClaimNumDF = pd.read_excel(mypath + "/" + "车险保单理赔样本.xlsx"
                           ).drop_duplicates(subset=['保单号'], ignore_index=True)
# 只保留出参和入参各维度
ClaimNumDF.drop(['保单号', '赔款金额'], axis='columns', inplace=True)

# 准备出入参矩阵用于模型训练
dumCN = pd.get_dummies(ClaimNumDF)        # 对各维度分箱做独热编码
cnXX = dumCN.iloc[:,1:]                   # 从样本独热编码表中截取入参
cnXX.drop(['车龄_0年', '座位数_5座及以下', '性别_男', '车系_车系3', 'NCD_新车新保'
         ], axis='columns', inplace=True) # 删除多余自由度
keyPDN = cnXX.columns.tolist()            # 获取各非基准伪变量名
keyPDN.insert(0, '基准')                  # 添加一个基准变量
cnXX = np.matrix(cnXX).T                  # 矩阵化
# 添加x0=1项，形成完整入参矩阵
cnXX = np.r_[matlib.ones((1,len(ClaimNumDF.index))), cnXX]
cnYY = np.matrix(dumCN.iloc[:,0]).T       # 从样本独热编码表中截取出参矩阵

# 通过梯度下降迭代训练回归参数
iterations=2000                                   # 设定迭代次数
step = 0.01                                       # 设定迭代步长
B = matlib.ones((cnXX.shape[0], 1))               # 初始化回归参数向量
modGList = []                                     # 初始化似然梯度模长列表
while iterations > 0:
    MU = np.exp(cnXX.T * B)                       # 计算当前估计期望向量
    G = cnXX * (cnYY - MU)                        # 计算当前梯度
    modG = np.linalg.norm(G)                      # 计算梯度模长
    stG = G / modG                                # 梯度做标准化，即使其模长为1
    B = B + step* stG                             # 按梯度上升方向做一步迭代
    modGList.append(modG)                         # 记录当前梯度模长
    iterations -= 1

# 绘制梯度下降曲线
plt.figure('出险次数梯度下降曲线')
plt.plot(range(0, len(modGList)), modGList)

allPDN = dumCN.columns.tolist()                   # 构造全体伪变量名列表
allPDN[0] = keyPDN[0]                             # 将第一个变量名也改成“基准”
FF = np.exp(B).T.tolist()[0]                      # 计算各分箱对应因子列表
# 初始化全分箱因子表，先都填写成1（因为e^0=1）
cnBinFactors = pd.DataFrame([[1]*len(dumCN.columns)], columns=allPDN)
# 在全分箱因子表中，按对应位置填写回归出的各分箱因子
for i, j in zip(keyPDN, range(len(FF))):
    cnBinFactors[i] = FF[j]

可见其结果与示例表中的很接近。

（三）训练案均赔款因子

使用伪造出的保单理赔清单训练案均赔款因子。根据前人经验，此时连结函数一般也选择对数函数。具体代码如下：

# =============================================================================
# 根据保单理赔样本训练案均赔款回归因子
# =============================================================================
# 从理赔列表中去重读取各单历次出险金额
ClaimPayDF = pd.read_excel(mypath + "/" + "车险保单理赔样本.xlsx"
                           ).drop_duplicates(subset=['保单号'], ignore_index=True)
# 只保留出参和入参各维度
ClaimPayDF.drop(['保单号', '出险次数'], axis='columns', inplace=True)
# 删除未出险保单
ClaimPayDF.drop(ClaimPayDF.index[ClaimPayDF['赔款金额']==0], inplace=True)
# 对索引重新排序
ClaimPayDF.reset_index(drop=True, inplace=True)

# 准备出入参矩阵用于模型训练
dumCP = pd.get_dummies(ClaimPayDF)        # 对各维度分箱做独热编码
cpXX = dumCP.iloc[:,1:]                   # 从样本独热编码表中截取入参
cpXX.drop(['车龄_0年', '座位数_5座及以下', '性别_男', '车系_车系3', 'NCD_新车新保'
         ], axis='columns', inplace=True) # 删除多余自由度
keyPDN = cpXX.columns.tolist()            # 获取各非基准伪变量名
keyPDN.insert(0, '基准')                  # 添加一个基准变量
cpXX = np.matrix(cpXX).T                  # 矩阵化
# 添加x0=1项，形成完整入参矩阵
cpXX = np.r_[matlib.ones((1,len(ClaimPayDF.index))), cpXX]
cpYY = np.matrix(dumCP.iloc[:,0]).T       # 从样本独热编码表中截取出参矩阵

# 通过梯度下降迭代训练回归参数
iterations=2000                                   # 设定迭代次数
step = 0.01                                       # 设定迭代步长
B = matlib.ones((cpXX.shape[0], 1))               # 初始化回归参数向量
modGList = []                                     # 初始化似然梯度模长列表
while iterations > 0:
    MU = np.exp(cpXX.T * B)                       # 计算当前估计期望向量
    # 计算当前梯度。由于伽马分布的形状参数alpha不影响梯度方向，故可设为1
    G = cpXX * np.multiply((cpYY - MU), 1/MU)
    modG = np.linalg.norm(G)                      # 计算梯度模长
    stG = G / modG                                # 梯度做标准化，即使其模长为1
    B = B + step* stG                             # 按梯度上升方向做一步迭代
    modGList.append(modG)                         # 记录当前梯度模长
    iterations -= 1

# 绘制梯度下降曲线
plt.figure('案均赔款梯度下降曲线')
plt.plot(range(0, len(modGList)), modGList)

allPDN = dumCP.columns.tolist()                   # 构造全体伪变量名列表
allPDN[0] = keyPDN[0]                             # 将第一个变量名也改成“基准”
FF = np.exp(B).T.tolist()[0]                      # 计算各分箱对应因子列表
# 初始化全分箱因子表，先都填写成1（因为e^0=1）
cpBinFactors = pd.DataFrame([[1]*len(dumCP.columns)], columns=allPDN)
# 在全分箱因子表中，按对应位置填写回归出的各分箱因子
for i, j in zip(keyPDN, range(len(FF))):
    cpBinFactors[i] = FF[j]

# 通过出险次数因子与案均赔款因子相乘汇总成赔付成本因子
binFactor = cnBinFactors.T
binFactor.rename(columns={0:'出险次数因子'}, inplace=True)
binFactor['赔款金额因子'] = cpBinFactors.T
binFactor['赔付成本因子'] = binFactor['出险次数因子'] * binFactor['赔款金额因子']

可见其结果与示例表中的很接近。根据复合随机变量全期望公式，整体赔付成本期望就等于出险次数期望乘以赔款金额期望，因此各维度各分箱的总赔付成本因子就等于出险次数因子乘以赔款金额因子。

（四）模型检验

# =============================================================================
# 回归检验
# =============================================================================
# 读取理赔列表
ClaimDF = pd.read_excel(mypath + "/" + "车险保单理赔样本.xlsx")
# 透视出各保单总赔款
policyPayDF = ClaimDF.pivot_table(values='赔款金额', index='保单号', aggfunc=np.sum)
policyPayDF.insert(0, '保单号', policyPayDF.index)
policyPayDF.reset_index(drop=True, inplace=True)

# 对原理赔列表做保单号去重
ClaimDF.drop_duplicates(subset=['保单号'], ignore_index=True, inplace=True)
# 删除原出险次数和理赔金额列
ClaimDF.drop(['赔款金额', '出险次数'], axis='columns', inplace=True)
# 将透视出的各保单总赔款连结到去重后的理赔列表，得到保单赔款表
policyPayDF = pd.merge(ClaimDF, policyPayDF, on='保单号', how='inner')
policyPayDF.reset_index(drop=True, inplace=True)

# 计算各板块赔付金额的均值和方差
meanCP = []                                 # 初始化各入参对应的赔付金额均值列表
varianceCP = []                             # 初始化各入参对应的赔付金额方差列表
for i in np.arange(0, len(policyPayDF.index), NUM):
    claimPay = policyPayDF.iloc[i:i+NUM]['赔款金额']
    mCP = sum(claimPay) / len(claimPay)            # 生成当前各出参样本平均数
    vCP = (sum((claimPay-mCP)**2) / (NUM-1))**0.5  # 生成当前各出参样本标准差
    meanCP.append(mCP)
    varianceCP.append(vCP)

# 构造由各板块对应的入参向量构成的入参矩阵
del ClaimDF['保单号']
x = pd.get_dummies(ClaimDF.drop_duplicates())  # 初始化非重复入参独热编码表
x.insert(0, '基准', [1]* len(x.index))         # 编码表列首添加一列1
X = np.matrix(x).T                             # 矩阵化各入参点

# 计算卡方统计量
BF = np.matrix(binFactor['赔付成本因子']).T     # 分箱因子矩阵化
ZZ = []                                        # 初始化标准正态分布平方统计量列表
for i in range(0, len(x)):
    # 由各入参点的出参观测均值、观测方差、估计出参、观测样本量构造的卡方统计量
    if meanCP[i] > 0:
        ZZ.append(((meanCP[i] - np.exp(X[:,i].T * np.log(BF)).A[0][0]) * (NUM**0.5)\
              / varianceCP[i])**2)
cp = 1 - st.chi2.cdf(sum(ZZ), df=len(ZZ))      # 计算当前回归检验的显著水平
CL = 0.05                                      # 设定最小能接受的显著水平
if cp > CL:
    print(binFactor)
    # 将回归因子保存成excel
    binFactor.to_excel(mypath + "/" + '赔付成本回归因子.xlsx')
    print(f"不可以在{round(100*CL,2)}%的显著水平下拒绝此回归模型，回归向量为：")
else:
    print(f"可以在{round(100*CL,2)}%的显著水平下拒绝此回归模型")

事实上如果此检验统计量的观测结果能落在上30%分位数的位置就已经拟合的很完美了。如果用于模型训练的数据更多，效果会更好。