首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >numpy/pandas瞎搞系列(一):OLS,WLS的numpy实现

numpy/pandas瞎搞系列(一):OLS,WLS的numpy实现

作者头像
量化小白
发布2020-08-20 16:03:55
3.2K0
发布2020-08-20 16:03:55
举报

python里很多模块都有OLS的实现,之前总结过一次,详见《从零开始学量化(五):用Python做回归》。今天这个是自己用numpy实现OLS,WLS的一些内容。

自己写的好处是比内置的要快一些,倒也不是说内置的代码写的不好,而是内置的函数一般比较完善,会算出来很多东西,但有的时候就只需要个回归系数,需要个预测值,结果内置函数给你整一堆t值,r方,mse,summary拖慢你的速度,也没啥意义。

01

OLS的beta

先从OLS开始,OLS比较常用的是statsmodels里的OLS模块,除此之外,numpy里本身也有一个函数能实现。这里从定义出发直接算一个,另外做一个简单测试对比numpy和statsmodels里的速度差异。

OLS的beta定义:

公式推导就省略了,随便找概率书都有,直接代码。

def getOLSBeta(x,y):
    x = np.hstack([x.reshape(-1,1),np.ones([x.shape[0],1])])
    beta = np.linalg.inv(x.T.dot(x)).dot(x.T).dot(y)
    return beta

接下来分别用这个函数和statsmodel的OLS模块算beta,做一个简单测试。

首先随机生成两个(10000,1)的向量作为x,y:

np.random.seed(1)
x = np.random.uniform(0,1, (10000,))
y = np.random.uniform(0,1, (10000,))

然后分别用这两个函数算beta,重复100次,计算总的时间,结果对比如下

OLS需要86735microseconds,numpy函数的结果:

自定义的函数需要14958microseconds,自定义的要快6-7倍。

02

WLS的beta

同样的道理,定义WLS的beta函数,这个就不做测试了,不用想都知道肯定是比statsmodel里的WLS更快一些。WLS的beta表达式:

python代码如下

def getWLSBeta(x,y,w):
    
    x = np.hstack([x.reshape(-1,1),np.ones([x.shape[0],1])])
    
    x1 = np.sqrt(w.reshape(-1,1)) * x
    y1 = np.sqrt(w) * y
    
    beta = np.linalg.inv(x1.T.dot(x1)).dot(x1.T).dot(y1)
    
    return beta

03

OLS的预测值

OLS的预测值,有两种,一般大家只看点预测,也就是拟合出来的值,这个很简单,不管是新来的点还是回归数据里的点,直接和beta乘一下就出来了,没啥难度。

还有一个是区间预测,给定一个置信度alpha,预测值是会有一个范围的,点预测实际上是预测区间的中点,预测区间有多宽,取决于方程本身的拟合程度,拟合的越好,区间越窄,说明预测的结果准确性越高。所以这里用numpy实现OLS的区间预测,推导也不写了,直接写公式

numpy计算的函数如下

def getOLSfitted(x,y,x_new = None,alpha = 0.05):
    
    x = np.hstack([x.reshape(-1,1),np.ones([x.shape[0],1])])
    y_pred_se = np.linalg.inv(x.T.dot(x))
    beta = y_pred_se.dot(x.T).dot(y)

    y_pred = x.dot(beta)
    sigma2_est = np.sum((y - y_pred)**2) / (x.shape[0] - x.shape[1])
    
    if x_new == None:
        x_new = x
        y_pred_new = y_pred
    else:
        x_new = np.hstack([x_new.reshape(-1,1),np.ones([x_new.shape[0],1])])
        y_pred_new = x_new.dot(beta)

    y_pred_se = sigma2_est + (x_new * np.dot(y_pred_se * sigma2_est,x_new.T).T).sum(1)
    y_pred_se = np.sqrt(y_pred_se)

    y_pred_lower = y_pred_new - stats.t.ppf(q = 1- alpha/2 , df = x.shape[0] - x.shape[1])*y_pred_se
    y_pred_upper = y_pred_new + stats.t.ppf(q = 1- alpha/2 , df = x.shape[0] - x.shape[1])*y_pred_se
    
    return y_pred_se,y_pred,y_pred_lower,y_pred_upper

解释一下,函数输入x,y为拟合方程的,x_new为用来预测的点,如果不输入的话返回的是x的拟合值,alpha是置信度。返回值有四个,se,y的点预测,y的区间预测下限,y的区间预测上限。写的比较粗糙,自己用的话可以再根据需求改一改。

另外statsmodel里也可以直接求OLS的预测区间,需要用到wls_prediction_std函数,所以还是之前的那个例子,做一个测试。

wls_prediction_std的结果

需要130651microseconds。

自定义函数的结果:

需要81749microseconds,要快一些,但还有提升空间。

03

WLS的预测值

WLS同样的道理,可以直接用wls_prediction_std计算,另外也自己定义了getWLSfitted函数,对比结果如下,还是之前的例子,再随机生成一个权重w。

np.random.seed(1)
x = np.random.uniform(0,1, (10000,))
y = np.random.uniform(0,1, (10000,))
w = np.random.uniform(0,1, (10000,))

wls_prediction_std函数的结果

需要166555microseconds。

自定义函数的结果

需要98703microseconds,所以还是自定义的要快一些,但还有提升空间。

本文所有代码包括getWLSfitted的函数定义请在后台回复“WLS“获取。

当然除了以上这些, 还可以自己定义算OLS 的t值,OLS的控制等等,这里就略过啦。

本文参与 腾讯云自媒体分享计划,分享自微信公众号。
原始发表:2020-08-18,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 量化小白躺平记 微信公众号,前往查看

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

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

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档