数学模型

Wang Haihua

🍈 🍉🍊 🍋 🍌

Logistic 模型的参数解释

在流行病学中, 经常需要研究某一疾病发生与不发生的可能性大小, 如一个人得流行性感冒相对于不得流行性感冒的可能性是多少, 对此常用 赔率来度量。赔率的具体定义如下:

定义1 一个随机事件 $A$ 发生的概率与其不发生的概率之比值称为 事件 $A$ 的赔率, 记为odds $(A)$, 即odds $(A)=\frac{P(A)}{P(\bar{A})}=\frac{P(A)}{1-P(A)}$ 。

如果一个事件 $A$ 发生的概率 $P(A)=0.75$, 则其不发生的概率 $P(\bar{A})=1-P(A)=0.25$, 所以事件 $A$ 的赔率 $\operatorname{odds}(A)=\frac{0.75}{0.25}=3$ 。这就是说, 事 件 $A$ 发生与不发生的可能性是 3 比 1。粗略地讲, 即在 4 次观测中有 3 次事 件 $A$ 发生而有一次 $A$ 不发生。例如，事件 $A$ 表示“投资成功”，那么 $\operatorname{odds}(A)=3$ 即表示投资成功的可能性是投资不成功的 3 倍。又例如, 事件 $\boldsymbol{B}$ 表示“客户 理赔事件", 且已知 $P(B)=0.25$, 则 $P(\bar{B})=0.75$, 从而事件 $B$ 的赔率 $\operatorname{odds}(B)=\frac{1}{3}$, 这表明发生客户理赔事件的风险是不发生的三分之一。用赔率 可很好地度量一些经济现象发生与否的可能性大小。

仍以上述“客户理赔事件”为例, 有时我们还需要研究某一群客户相对于 另一群客户发生客户理赔事件的风险大小, 如职业为司机的客户群相对于 职业是教师的客户群发生客户理赔事件的风险大小, 这需要用到赔率比的 概念。 定义 $12.4$ 随机事件 $A$ 的赔率与随机事件 $B$ 的赔率之比值称为事件 $A$ 对事件 $B$ 的赔率比, 记为 $\operatorname{OR}(A, B)$, 即 $\operatorname{OR}(A, B)=\operatorname{odds}(A) / \operatorname{odds}(B)$ 。

若记 $A$ 是职业为司机的客户发生理赔事件，记 $B$ 是职业为教师的客户发 生理赔事件, 又已知 $\operatorname{odds}(A)=\frac{1}{20}, \operatorname{odds}(B)=\frac{1}{30}$, 则事件 $A$ 对事件 $B$ 的赔率 比 $\operatorname{OR}(A, B)=\operatorname{odds}(A) / \operatorname{odds}(B)=1.5$ 。这表明职业为司机的客户发生理赔的 赔率是职业为教师的客户的 $1.5$ 倍。 应用 Logistic 回归可以方便地估计一些事件的赔率及多个事件的赔率 比。下面仍以例 $12.7$ 为例来说明 Logistic 回归在这方面的应用。

例 $12.8$ (续例 12.7) 房地产商希望能估计出一个家庭年收入为 9 万元 的客户其签订意向后最终买房与不买房的可能性大小之比值, 以及一个家 庭年收入为 9 万元的客户其签订意向后最终买房的赔率是年收入为 8 万元 客户的多少倍。

\begin{aligned} &\text { 表 } 12.4 \text { 签订购房意向和最终买房的客户数据 }\\ &\begin{array}{cccccc} \hline & \text { 家庭年收入号 } & \text { 签订意向书 } & \text { 实际购房 } & \text { 实际购房比例 } & \text { 逻辑变换 } \\ & \text { (万元） } & \text { 人数 } n_{i} & \text { 人数 } m_{i} & \hat{p}_{i}=\frac{m_{i}}{n_{i}} & y_{i}^{*}=\ln \left(\frac{\hat{p}_{i}}{1-\hat{p}_{i}}\right) \\ \hline 1 & 1.5 & 25 & 8 & 0.32 & -0.7538 \\ 2 & 2.5 & 32 & 13 & 0.4063 & -0.3795 \\ 3 & 3.5 & 58 & 26 & 0.4483 & -0.2076 \\ 4 & 4.5 & 52 & 22 & 0.4231 & -0.3102 \\ 5 & 5.5 & 43 & 20 & 0.4651 & -0.1398 \\ 6 & 6.5 & 39 & 22 & 0.5641 & 0.2578 \\ 7 & 7.5 & 28 & 16 & 0.5714 & 0.2877 \\ 8 & 8.5 & 21 & 12 & 0.5714 & 0.2877 \\ 9 & 9.5 & 15 & 10 & 0.6667 & 0.6931 \\ \hline \end{array} \end{aligned}

解 由例 $12.7$ 中所得的模型 (12.29) 得 $$ \ln \left(\frac{\hat{p}}{1-\hat{p}}\right)=-0.8863+0.1558 x, $$ 因此 $$ \frac{\hat{p}}{1-\hat{p}}=e^{-0.8863+0.1558 x} \text {. } $$ 将 $x=x_{0}=9$ 代入上式, 得一个家庭年收入为 9 万元的客户其签订意向 后最终买房与不买房的可能性大小之比值为 $$ \text { odds(年收入9万) }=\frac{\hat{p}_{0}}{1-\hat{p}_{0}}=e^{-0.8863+0.1558 \times 9}=1.6752 \text {. } $$ 这说明一个家庭年收入为 9 万元的客户其签订意向后最终买房的可能性是 不买房可能性的 $1.6752$ 倍。

另外, 由 (12.30) 式还可得 OR(年收入9万元, 年收入 8 万元 $)=\frac{e^{-0.8863+0.1558 \times 9}}{e^{-0.8863+0.1558 \times 8}}=1.1686 .$ 所以一个家庭年收入为 9 万元的客户其签订意向后最终买房的赔率是年收 入为 8 万元客户的 $1.1686$ 倍。

一般地, 如果 Logistic 模型 $(12.23)$ 的参数估计为 $\hat{\beta}_{0}, \hat{\beta}_{1}, \cdots, \hat{\beta}_{m}$, 则在 $x_{1}=x_{01}, x_{2}=x_{02}, \cdots, x_{m}=x_{0 m}$ 条件下事件赔率的估计值为 $$ \frac{\hat{\boldsymbol{p}}_{0}}{\mathbf{1 - \hat { p } _ { 0 }}}=\boldsymbol{e}^{\hat{\beta}_{0}+\sum_{j=1}^{m} \hat{\beta}_{j} x_{0 j}} . $$ 如果记 $X_{A}=\left[1, x_{A 1}, x_{A 2}, \cdots, x_{A m}\right]^{T}, X_{B}=\left[1, x_{B 1}, x_{B 2}, \cdots, x_{B m}\right]^{T}$, 并将相应条 件下的事件仍分别记为 $X_{A}$ 和 $X_{B}$, 则事件 $X_{A}$ 对 $X_{B}$ 赔率比的估计可由下式获 得 $$ \operatorname{OR}\left(X_{A}, X_{B}\right)=e^{\sum_{j=1}^{m} \hat{\beta}_{j}\left(x_{i j}-x_{\left.B_{j}\right)}\right.} . $$

用 statsmodels 库函数求解

例 $12.9$ 企业到金融商业机构贷款, 金融商业机构需要对企业进行评估。 评估结果为 0,1 两种形式, 0 表示企业两年后破产, 将拒绝贷款, 而 1 表 示企业两年后具备还款能力, 可以贷款。在表 $12.5$ 中, 已知前 20 家企业的 三项评价指标值和评估结果, 试建立模型对其他 2 家企业（企业 21、22) 进行评估。

\begin{array}{cccccc} \hline \text { 企业编号 } & x_{1} & x_{2} & x_{3} & y & y \text { 的预测值 } \\ \hline 1 & -62.3 & -89.5 & 1.7 & 0 & 0 \\ 2 & 3.3 & -3.5 & 1.1 & 0 & 0 \\ 3 & -120.8 & -103.2 & 2.5 & 0 & 0 \\ 4 & -18.1 & -28.8 & 1.1 & 0 & 0 \\ 5 & -3.8 & -50.6 & 0.9 & 0 & 0 \\ 6 & -61.2 & -56.2 & 1.7 & 0 & 0 \\ 7 & -20.3 & -17.4 & 1 & 0 & 0 \\ 8 & -194.5 & -25.8 & 0.5 & 0 & 0 \\ 9 & 20.8 & -4.3 & 1 & 0 & 0 \\ 10 & -106.1 & -22.9 & 1.5 & 0 & 0 \\ 11 & 43 & 16.4 & 1.3 & 1 & 1 \\ 12 & 47 & 16 & 1.9 & 1 & 1 \\ 13 & -3.3 & 4 & 2.7 & 1 & 1 \\ 14 & 35 & 20.8 & 1.9 & 1 & 1 \\ 15 & 46.7 & 12.6 & 0.9 & 1 & 1 \\ 16 & 20.8 & 12.5 & 2.4 & 1 & 1 \\ 17 & 33 & 23.6 & 1.5 & 1 & 1 \\ 18 & 26.1 & 10.4 & 2.1 & 1 & 1 \\ 19 & 68.6 & 13.8 & 1.6 & 1 & 1 \\ 20 & 37.3 & 33.4 & 3.5 & 1 & 1 \\ 21 & -49.2 & -17.2 & 0.3 & & 0 \\ 22 & 40.6 & 26.4 & 1.8 & 1 \\ \hline \end{array}

解 对于该问题, 可以用 Logistic 模型来求解。建立如下的 Logistic 回 归模型 $$ p=P\{y=1\}=\frac{1}{1+e^{-\left(\beta_{0}+\sum_{j=1}^{3} \beta_{j} x_{j}\right)}} . $$ 记 $x_{i j}(i=1,2, \cdots, 20 ; j=1,2,3)$ 分别为变量 $x_{j}(j=1,2,3)$ 的 20 个观测值, $y_{i}(i=1,2, \cdots, 20)$ 是 20 个 $y$ 的观测值。我们使用最大似然估计法, 求模型中 的参数 $\beta_{0}, \beta_{1}, \beta_{2}, \beta_{3}$, 即求参数 $\beta_{0}, \beta_{1}, \beta_{2}, \beta_{3}$ 使得似然函数 $$ \ln L(\beta)=\sum_{i=1}^{20}\left[y_{i}\left(\beta_{0}+\sum_{j=1}^{3} \beta_{j} x_{i j}\right)-\ln \left(1+e^{\beta_{0}+\sum_{j=1}^{3} \beta_{j} x_{i j}}\right]\right] $$ 达到最大值。

利用 Python 的 statsmodels 库函数求得 $$ \beta_{0}=0, \beta_{1}=-0.3497, \quad \beta_{2}=3.2290, \beta_{3}=2.2372 . $$ 因而得到的 Logistic 回归模型为 $$ \left\{\begin{array}{l} p=\frac{1}{1+e^{-\left(-0.3497 x_{1}+3.2290 x_{2}+2.2372 x_{3}\right)}}, \\ y= \begin{cases}0, & p \leq \mathbf{0 . 5}, \\ 1, & p>0.5 .\end{cases} \end{array}\right. $$ 利用已知数据对上述 Logistic 模型进行检验, 准确率达到 100\%, 说明模型 的准确率较高, 可以用来预测新企业的还款能力。2 个新企业的预测结果见 表 $12.5$ 的最后 1 列, 即企业 21 拒绝贷款, 企业 22 可以贷款。

代码

import numpy as np
import statsmodels.api as sm
a=np.loadtxt("data/loan.txt")
n=a.shape[1] #提取矩阵的列数
x=a[:,:n-1]; y=a[:,n-1]
md=sm.Logit(y,x)
md=md.fit(method="bfgs")  #这里必须使用bfgs方法，使用默认牛顿方法出错
print(md.params,'\n----------\n'); print(md.summary2())
print(md.predict([[-49.2,-17.2,0.3],[40.6,26.4,1.8]]))  #求预测值

用 sklearn 库函数求解

例 12.10（续例 12.9）用 sklearn 库函数求解例 12.9。 解 求得的 Logistic 回归模型为 $$ \left\{\begin{array}{l} p=\frac{1}{1+e^{-\left(-0.3906-0.0507 x_{1}+0.6707 x_{2}+0.1051 x_{3}\right)}} \\ y= \begin{cases}0, & p \leq 0.5, \\ 1, & p>0.5 .\end{cases} \end{array}\right. $$ 利用已知数据对上述 Logistic 模型进行检验, 准确率达到 100\%, 说明模型 的准确率较高, 可以用来预测新企业的还款能力。 2 个新企业的预测结果为： 企业 21 拒绝贷款，企业 22 可以贷款。

代码

import numpy as np
from sklearn.linear_model import LogisticRegression
a=np.loadtxt("data/loan.txt")
n=a.shape[1]  #提取矩阵的列数
x=a[:,:n-1]; y=a[:,n-1]
md=LogisticRegression(solver='lbfgs')
md=md.fit(x,y)
print(md.intercept_,md.coef_)
print(md.predict(x))   #检验预测模型
print(md.predict([[-49.2,-17.2,0.3],[40.6,26.4,1.8]]))  #求预测值