numpy/pandas瞎搞系列(一):OLS,WLS的numpy实现
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的控制等等,这里就略过啦。
- HDU 1250 Hat's Fibonacci
- Scrapy在Ubuntu下的安装与配置
- Selenium2+python自动化20-引入unittest框架
- HDU 1002 A + B Problem II(高精度加法(C++/Java))
- POJ 1018 Communication System
- POJ 1017 Packets
- Codeforces 725B Food on the Plane
- Codefoces 723B Text Document Analysis
- Codefoces 723A The New Year: Meeting Friends
- ECJTUACM16 Winter vacation training #1 题解&源码
- 信息学奥赛一本通算法(C++版)基础算法:高精度计算
- 看破欧拉函数的奥秘
- 线段树入门总结
- 从零基础学三分查找
- JavaScript 教程
- JavaScript 编辑工具
- JavaScript 与HTML
- JavaScript 与Java
- JavaScript 数据结构
- JavaScript 基本数据类型
- JavaScript 特殊数据类型
- JavaScript 运算符
- JavaScript typeof 运算符
- JavaScript 表达式
- JavaScript 类型转换
- JavaScript 基本语法
- JavaScript 注释
- Javascript 基本处理流程
- Javascript 选择结构
- Javascript if 语句
- Javascript if 语句的嵌套
- Javascript switch 语句
- Javascript 循环结构
- Javascript 循环结构实例
- Javascript 跳转语句
- Javascript 控制语句总结
- Javascript 函数介绍
- Javascript 函数的定义
- Javascript 函数调用
- Javascript 几种特殊的函数
- JavaScript 内置函数简介
- Javascript eval() 函数
- Javascript isFinite() 函数
- Javascript isNaN() 函数
- parseInt() 与 parseFloat()
- escape() 与 unescape()
- Javascript 字符串介绍
- Javascript length属性
- javascript 字符串函数
- Javascript 日期对象简介
- Javascript 日期对象用途
- Date 对象属性和方法
- Javascript 数组是什么
- Javascript 创建数组
- Javascript 数组赋值与取值
- Javascript 数组属性和方法
- laravel实现按时间日期进行分组统计方法示例
- Pytorch 使用CNN图像分类的实现
- PHP实现无限极分类的两种方式示例【递归和引用方式】
- 记录模型训练时loss值的变化情况
- phpstorm 配置xdebug的示例代码
- 利用Python实现Excel的文件间的数据匹配功能
- PHP设计模式之简单工厂和工厂模式实例分析
- PHP实现数据四舍五入的方法小结【4种方法】
- 如何在Windows中安装多个python解释器
- PHP设计模式之抽象工厂模式实例分析
- 使用python matploblib库绘制准确率,损失率折线图
- Django REST Swagger实现指定api参数
- matplotlib.pyplot.matshow 矩阵可视化实例
- php+mysql开发的最简单在线题库(在线做题系统)完整案例
- python中元组的用法整理