前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
圈层
工具
发布
首页
学习
活动
专区
圈层
工具
MCP广场
社区首页 >专栏 >生存分析——泊松回归(LightGBM)实现生存分析(四)

生存分析——泊松回归(LightGBM)实现生存分析(四)

作者头像
悟乙己
发布于 2021-12-07 08:00:19
发布于 2021-12-07 08:00:19
1.6K00
代码可运行
举报
文章被收录于专栏:素质云笔记素质云笔记
运行总次数:0
代码可运行

You can find the full article here

来看一个比较特殊的Survival分析建模的案例,利用的是半参模型:Poisson Regression 具体参考文章:Survival Analysis with LightGBM plus Poisson Regression 里面的建模思路非常有意思,不适合工业落地,不过咨询公司的data scientist看过来~

1 Poisson Regression

1.1 松泊分布与泊松回归

参考:什么是松泊分布?泊松回归可以用来做什么?

试想一下,你现在就站在一个人流密集的马路旁,打算收集闯红灯的人群情况(?)。 首先,利用秒表和计数器,一分钟过去了,有5个人闯红灯; 第二分钟有4个人;而下一分钟有4个人。 持续记录下去,你就可以得到一个模型,这便是“泊松分布”的原型。

除此以外,现实生活中还有很多情况是服从泊松分布的:

  • 10分钟内从ATM中取钱的人数
  • 一天中发生车祸的次数
  • 每100万人中患癌症的人数
  • 单位面积土地内昆虫的数目

Poisson模型(泊松回归模型)是用于描述单位时间、单位面积或者单位容积内某事件发现的频数分布情况, 通常用于描述稀有事件(即小概率)事件发生数的分布。

上述例子中都明显的一个特点: 低概率性,以及单位时间(或面积、体积)内的数量。 通常情况下,满足以下三个条件时,可认为数据满足Poisson分布:

  • (1) 平稳性:发生频数的大小,只与单位大小有关系(比如1万为单位,或者100万为单位时患癌症人数不同);
  • (2) 独立性:发生频数的大小,各个数之间没有影响关系,即频数数值彼此独立没有关联关系; 比如前1小时闯红灯的人多了,第2小时闯红灯人数并不会受影响;
  • (3) 普通性:发生频数足够小,即低概率性。

如果数据符合这类特征时,而又想研究X对于Y的影响(Y呈现出Poisson分布); 此时则需要使用Poisson回归,而不是使用常规的线性回归等。

1.2 LightGBM 实现泊松回归的案例

参考来源:https://github.com/Microsoft/LightGBM/issues/807

代码语言:javascript
代码运行次数:0
运行
AI代码解释
复制
import lightgbm as lgb
import numpy as np
import pandas as pd

n=100000
lam = .01
X = np.floor(np.random.lognormal(size=(n,2))).astype(int)
y = np.maximum(X[:,0],X[:,1])+np.random.poisson(lam=lam, size=n)

train_inds = np.arange(int(n/3))
val_inds =   np.arange(int(n/3), int(2*n/3))
test_inds =   np.arange(int(2*n/3), int(n))

X_test, y_test = X[test_inds,:], y[test_inds]

ds = lgb.Dataset(X,y, categorical_feature=[1])

ds_train = ds.subset(train_inds)
ds_val = ds.subset(val_inds)

params = {'objective':'poisson',
         'metric':'rmse',
         'learning_rate':.1
         }
gbm = lgb.train(params, ds_train, num_boost_round=300, early_stopping_rounds=20, valid_sets=[ds_val, ds_train], 
                verbose_eval=100, categorical_feature=[1])

yhat = gbm.predict(X_test)
print('neg obs:', len(yhat[yhat<0]))

y是需要服从poisson分布的

2 数据解读

2.1 数据样式

数据集解释:美国Washington, D.C.的一个共享单车公司 数据来源:bike-sharing-dataset 数据集的中文字段解释: 参考:Capital Bikeshare (美国Washington, D.C.的一个共享单车公司)提供的自行车数据上进行回归分析 具体实现 + 代码:Survival_LGBM-github

代码语言:javascript
代码运行次数:0
运行
AI代码解释
复制
# 字段说明Instant 记录号

Dteday:日期
Season:季节
    1=春天
    2=夏天
    3=秋天
    4=冬天
yr:年份,(0: 2011, 1:2012)
mnth:月份( 1 to 12)
hr:小时 (0 to 23) (只在 hour.csv 有,作业忽略此字段)
holiday:是否是节假日
weekday:星期中的哪天,取值为 06
workingday:是否工作日
    1=工作日 (非周末和节假日)
    0=周末
weathersit:天气
1:晴天,多云
2:雾天,阴天
3:小雪,小雨
4:大雨,大雪,大雾
temp:气温摄氏度
atemp:体感温度
hum:湿度
windspeed:风速

y值
-   casual:非注册用户个数
-  registered:注册用户个数
-  cnt:给定日期(天)时间(每小时)总租车人数,响应变量 y

现在数据变成:

也就是把第一行数据,拆分成330行,新增了两列:

  • count_so_far就是人员计数从[0,1,2,3,4,6,…,330],
  • stop,就是人数终止计数

这样构造的一个特殊的数据结构,造就了该模型object的独特,原文自评:

We solve a not classical survival problem where we have not to estimate the probability of surviving past time, but we estimate the probability that a specific event occurs at the end of the day. This is a simple example but it shows how to applicate survival modeling techniques, with classical instruments, in a not common scenario where our duty is to estimate a probability density function.

这里就把问题变成了,2011/1/1 这天, 在casual users = 0的时候,不会停止增长;在casual users = 1,不会停止增长;…;在casual users = 330,停止增长 最后,如何来预测临时用户的人数? 用的是一个“可能在什么数量上停止”的曲线,这里是把 hazard = 临时用户的增长量,会拟合一个达到预测数值的可能性趋势。

这里跟生存分析 以及常规回归的差异:

  • 常规回归模型,y~x1+x2… 根据特征直接预测临时用户人数,点估计;这里的Poisson有点类似加强版区间估计,就叫趋势估计?
  • 生存分析,
    • 生存函数使用的是KM曲线,计算不同生存时间区间下的存活率;
    • 风险函数,计算不同生存时间区间下的风险值

2.2 衡量模型指标

2.2.1 CRPS

我们对结果的评价是用这种任务的标准分数来进行的。连续排序概率评分(Continuous ranking Probability Score, CRPS)将MAE推广到概率预测的情况。

在涉及概率预测的情况下,CPRS是应用最广泛的精度指标之一。它需要对所涉及的整个概率函数进行评估,因此需要对每个每日样本的整个生存曲线进行估计。

PDF / CDF

  • 概率密度函数probability density function(pdf)
  • 累积分布函数(Cumulative Distribution Function),又叫分布函数

2.2 训练与解读

数据的延展:

代码语言:javascript
代码运行次数:0
运行
AI代码解释
复制
### EXPAND TRAIN DATA ###

X_train['count_so_far'] = X_train.apply(lambda x: np.arange(x.casual), axis=1)
X_train['stop'] = X_train.apply(lambda x: np.append(np.zeros(x.casual-1), 1), axis=1)
X_train = X_train.apply(pd.Series.explode)   # 数据扩充的方式

X_train['count_so_far'] = X_train.count_so_far.astype(int)
X_train['stop'] = X_train.stop.astype(int)

y_train = X_train['stop']
X_train = X_train[columns+['count_so_far']]

print(X_train.shape)
X_train

### EXPAND VALID DATA FOR PREDICTION ###

X_val_surv = df[(df.dteday >= datetime(year=2012, month=1, day=1))&(df.dteday < datetime(year=2012, month=5, day=1))].copy()

X_val_surv['count_so_far'] = X_val_surv.apply(lambda x: np.arange(max_count), axis=1)
X_val_surv['stop'] = X_val_surv.apply(lambda x:np.append(np.zeros(x.casual), np.ones(max_count-x.casual)), axis=1)
X_val_surv = X_val_surv.apply(pd.Series.explode)

X_val_surv['count_so_far'] = X_val_surv.count_so_far.astype(int)
X_val_surv['stop'] = X_val_surv.stop.astype(int)

y_val_surv = X_val_surv['stop']
X_val_surv = X_val_surv[columns+['count_so_far']]

print(X_val_surv.shape)

具体的已经在2.1 描述过了;这里要额外来看一下max_count,这个在train数据集中不会出现, 也就是今天有多少临时用户,就是多少X_train.apply(lambda x: np.arange(x.casual), axis=1) 但你会看到,X_val_surv是X_val_surv.apply(lambda x: np.arange(max_count), axis=1),这里就是临时用户的最大值设定为4000,这个值根据临时用户max值来取的,就是临时用户上线。

代码语言:javascript
代码运行次数:0
运行
AI代码解释
复制
params = {
         'objective':'poisson',
         'num_leaves':30, 
         'learning_rate': 0.001,
         'feature_fraction': 0.8,
         'bagging_fraction': 0.9,
         'bagging_seed': 33,
         'poisson_max_delta_step': 0.8,
         'metric': 'poisson'
         }

### FIT LGBM WITH POISSON LOSS ### 

trn_data = lgb.Dataset(X_train, label=y_train)
val_data = lgb.Dataset(X_val, label=y_val)

model = lgb.train(params, trn_data, num_boost_round=1000,
                  valid_sets = [trn_data, val_data],
                  verbose_eval=50, early_stopping_rounds=150)

这里X_train的shape为:(247252, 13),那么可以知道训练集是非常大的; 需要重置数据,可到:

代码语言:javascript
代码运行次数:0
运行
AI代码解释
复制
### PREDICT HAZARD FUNCTION ON VALIDATION DATA AND TRANSFORM TO SURVIVAL ###

p_val_hz = model.predict(X_val_surv).reshape(-1,max_count)
p_val = 1-np.exp(-np.cumsum(p_val_hz, axis=1))
X_val_surv.shape
>>> (484000, 12)
p_val.shape
>>> (121, 4000)

p_val就回归正常,代表着每一天,不同人群数量的概率,直接上图:

这天在500左右,达到峰值了,那么这天临时用户的预测值就在500左右了。

简单的模型检验:

代码语言:javascript
代码运行次数:0
运行
AI代码解释
复制
### CRPS ON VALIDATION DATA ###
crps(t_val, p_val)
>>> 0.17425293919055515
### CRPS ON VALIDATION DATA WITH BASELINE MODEL ###
crps(t_val, np.repeat(cdf, len(t_val)).reshape(-1,max_count))

3 同等lightGBM分类测试

https://blog.csdn.net/wang263334857/article/details/81836578

来看一下同一份数据测试出来的结果如何,后续不贴太多,只贴一下我测试的代码,放在了之前的一个项目下面:Survival_Poisson_Regression

本文参与 腾讯云自媒体同步曝光计划,分享自作者个人站点/博客。
原始发表:2021/08/13 ,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 作者个人站点/博客 前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
暂无评论
推荐阅读
编辑精选文章
换一批
机器学习第5天:逻辑回归
简单来说, 逻辑回归(Logistic Regression)是一种用于解决二分类(0 or 1)问题的机器学习方法,用于估计某种事物的可能性。比如某用户购买某商品的可能性,某病人患有某种疾病的可能性,以及某广告被用户点击的可能性等。
K同学啊
2019/01/22
5780
机器学习还能预测心血管疾病?没错,我用Python写出来了
全球每年约有1700万人死于心血管疾病,当中主要表现为心肌梗死和心力衰竭。当心脏不能泵出足够的血液来满足人体的需要时,就会发生心力衰竭,通常由糖尿病、高血压或其他心脏疾病引起。
CDA数据分析师
2020/09/14
2.9K2
机器学习还能预测心血管疾病?没错,我用Python写出来了
公交车总迟到?你大概掉进了“等待时间悖论&quot;
你到了车站,准备搭乘声称每10分钟一班的公交车。你盯着你的手表留意着时间,结果公交车终于在11分钟后到来。
磐创AI
2018/12/18
3920
机器学习第5天:逻辑回归
简单来说, 逻辑回归(Logistic Regression)是一种用于解决二分类(0 or 1)问题的机器学习方法,用于估计某种事物的可能性。比如某用户购买某商品的可能性,某病人患有某种疾病的可能性,以及某广告被用户点击的可能性等。
三猫
2021/10/27
6020
机器学习第5天:逻辑回归
剑指LightGBM和XGboost!斯坦福发表NGBoost算法
Stanford ML Group 最近在他们的论文中发表了一个新算法,其实现被称为 NGBoost。该算法利用自然梯度将不确定性估计引入到梯度增强中。本文试图了解这个新算法,并与其他流行的 boosting 算法 LightGBM 和 XGboost 进行比较,以了解它在实践中是如何工作的。
AI研习社
2019/11/20
1.1K0
剑指LightGBM和XGboost!斯坦福发表NGBoost算法
LightGBM大战XGBoost,谁将夺得桂冠?
如果你是一个机器学习社区的活跃成员,你一定知道 **提升机器**(Boosting Machine)以及它们的能力。提升机器从AdaBoost发展到目前最流行的XGBoost。XGBoost实际上已经成为赢得在Kaggle比赛中公认的算法。这很简单,因为他极其强大。但是,如果数据量极其的大,XGBoost也需要花费很长的时间去训练。
机器学习算法工程师
2018/07/27
1.6K0
LightGBM大战XGBoost,谁将夺得桂冠?
XGBoost和LightGBM
这两个模型都属于集成学习中的树模型,每个机器学习模型都有它特定的应用场景,不同的数据集适合用到的模型是不一样的。
算法之名
2021/10/09
1.3K0
XGBoost和LightGBM
刷分神器,使用hyperopt实现lightgbm自动化调参!
Hyperopt是最受欢迎的调参工具包,目前在github上已经获得star数量5.8k,在kaggle天池等比赛中经常出现它的身影。
lyhue1991
2021/08/06
2.6K0
刷分神器,使用hyperopt实现lightgbm自动化调参!
Lightgbm 回归预测
铁柱在2018年11月底发了一篇 LSTM 回归预测模型 ,现在改用Lightgbm模型。本篇文章偏工程,需要读者了解Python关于Class的语法,理论部分也会在后续的文章中介绍.
MeteoAI
2019/07/24
5.4K0
Metropolis Hastings采样和贝叶斯泊松回归Poisson模型
在本文中,我想向你展示如何使用R的Metropolis采样从贝叶斯Poisson回归模型中采样。
拓端
2021/12/21
7620
Metropolis Hastings采样和贝叶斯泊松回归Poisson模型
python lasso回归分析_解析python实现Lasso回归「建议收藏」
X = np.random.randn(n_samples, n_features) # randn(…)产生的是正态分布的数据
全栈程序员站长
2022/09/05
1.5K0
# LightGBM大战XGBoost,谁将夺得桂冠?
  如果你是一个机器学习社区的活跃成员,你一定知道 提升机器(Boosting Machine)以及它们的能力。提升机器从AdaBoost发展到目前最流行的XGBoost。XGBoost实际上已经成为赢得在Kaggle比赛中公认的算法。这很简单,因为他极其强大。但是,如果数据量极其的大,XGBoost也需要花费很长的时间去训练。
机器学习算法工程师
2019/01/31
9670
# LightGBM大战XGBoost,谁将夺得桂冠?
广义线性模型应用举例之泊松回归及R计算
在前文“广义线性模型”中,提到广义线性模型(GLM)可概括为服务于一组来自指数分布族的响应变量的模型框架,正态分布、指数分布、伽马分布、卡方分布、贝塔分布、伯努利分布、二项分布、负二项分布、多项分布、泊松分布、集合分布等都属于指数分布族,并通过极大似然估计获得模型参数。
王诗翔呀
2020/07/21
9.2K0
广义线性模型应用举例之泊松回归及R计算
十大机器学习库,超全面总结!!
也毫无疑问,Pytorch的同学以压倒性的话语权霸屏。其实无论使用哪种框架,适合自己适合项目是最合适的。
Python编程爱好者
2024/05/22
7060
十大机器学习库,超全面总结!!
【机器学习实战】kaggle背包价格预测(堆叠的实战用法)
该竞赛的数据集是学生背包价格预测数据集中训练的深度学习模型生成的。特征分布与原始分布接近但不完全相同。请随意将原始数据集用作本竞赛的一部分,既要探索差异,又要查看将原始内容纳入训练是否可以改善模型性能。
机器学习司猫白
2025/02/15
1380
基于LightGBM的UCI信贷数据建模(完整Python代码)
本文是UCI数据集建模的第3篇,第一篇是数据的探索性分析EDA部分,第二篇是基于LightGBM模型的baseline。
皮大大
2024/06/30
4150
数据挖掘实践(金融风控):金融风控之贷款违约预测挑战赛(下篇)xgboots/lightgbm/Catboost等模型--模型融合:stacking、blend
通过组合多个学习器来完成学习任务,通过集成方法,可以将多个弱学习器组合成一个强分类器,因此集成学习的泛化能力一般比单一分类器要好。
汀丶人工智能
2023/05/17
4.1K0
数据挖掘实践(金融风控):金融风控之贷款违约预测挑战赛(下篇)xgboots/lightgbm/Catboost等模型--模型融合:stacking、blend
用于时间序列数据的泊松回归模型
泊松和类泊松回归模型常用于基于计数的数据集,即包含整数计数的数据。例如,每小时走进医院急诊室的人数就是一个这样的数据集。
deephub
2021/04/30
2.2K0
Python业务分析实战|共享单车数据挖掘
自行车共享系统是传统自行车租赁的新一代,从注册会员、租赁到归还的整个过程都是自动化的。通过这些系统,用户可以很容易地从一个特定的位置租用自行车,并在另一个位置归还。目前,全球大约有500多个共享单车项目,这些项目由50多万辆自行车组成。今天,由于它们在交通、环境和健康问题上的重要作用,人们对这些系统产生了极大的兴趣。
数据STUDIO
2021/11/17
1.7K0
Python业务分析实战|共享单车数据挖掘
机器学习实战 | LightGBM建模应用详解
教程地址:http://www.showmeai.tech/tutorials/41
ShowMeAI
2022/03/21
3.1K0
机器学习实战 | LightGBM建模应用详解
推荐阅读
相关推荐
机器学习第5天:逻辑回归
更多 >
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档