怎么用python自制计算公式_如何使用Python和Numpy计算r平方?

我最初发布下面的基准是为了推荐numpy.corrcoef,愚蠢地没有意识到原来的问题已经使用了corrcoef,实际上是在询问高阶多项式拟合。我已经使用statsmodels为多项式r-squared问题添加了一个实际的解决方案,并且我已经离开了原始的基准测试,虽然偏离主题,但对某些人来说可能是有用的。

statsmodels能够直接计算多项式拟合的r^2,这里有2种方法......

import statsmodels.api as sm

import stasmodels.formula.api as smf

# Construct the columns for the different powers of x

def get_r2_statsmodels(x, y, k=1):

xpoly = np.column_stack([x**i for i in range(k+1)])

return sm.OLS(y, xpoly).fit().rsquared

# Use the formula API and construct a formula describing the polynomial

def get_r2_statsmodels_formula(x, y, k=1):

formula = 'y ~ 1 + ' + ' + '.join('I(x**{})'.format(i) for i in range(1, k+1))

data = {'x': x, 'y': y}

return smf.ols(formula, data).fit().rsquared

为了进一步利用statsmodels,还应该查看拟合的模型摘要,可以在Jupyter / IPython笔记本中打印或显示为丰富的HTML表。除了rsquared之外,结果对象还提供对许多有用的统计指标的访问。

model = sm.OLS(y, xpoly)

results = model.fit()

results.summary()

以下是我原来的答案,我对各种线性回归r ^ 2方法进行了基准测试...

问题中使用的corrcoef函数仅计算单个线性回归的相关系数r,因此它不能解决高阶多项式拟合的问题r^2。然而,对于它的价值,我发现对于线性回归,它确实是计算r的最快和最直接的方法。

def get_r2_numpy_corrcoef(x, y):

return np.corrcoef(x, y)[0, 1]**2

这些是我通过比较1000个随机(x,y)点的方法得到的时间结果:

纯Python(直接r计算)

1000循环,最佳3:每循环1.59毫秒

Numpy polyfit(适用于n次多项式拟合)

1000个循环,每循环最佳3:326μs

Numpy手册(直接r计算)

10000循环,最佳3:每循环62.1μs

Numpy corrcoef(直接r计算)

10000循环,最佳3:每循环56.6μs

Scipy(线性回归以r为输出)

1000个循环,最佳3:676μs/循环

Statsmodels(可以做n次多项式和许多其他拟合)

1000个循环,最佳3:每循环422μs

corrcoef方法使用numpy方法以“手动”方式勉强计算r ^ 2。它比polyfit方法快5倍,比scipy.linregress快12倍。只是为了加强numpy为你做的事情,它比纯蟒蛇快28倍。我并不精通像numba和pypy这样的东西,所以其他人不得不填补这些空白,但我认为这很有说服力,corrcoef是计算简单线性回归的最佳工具。

这是我的基准测试代码。我从Jupyter笔记本上复制粘贴(很难不称它为IPython笔记本......),所以如果有什么事情发生,我道歉。 %timeit magic命令需要IPython。

import numpy as np

from scipy import stats

import statsmodels.api as sm

import math

n=1000

x = np.random.rand(1000)*10

x.sort()

y = 10 * x + (5+np.random.randn(1000)*10-5)

x_list = list(x)

y_list = list(y)

def get_r2_numpy(x, y):

slope, intercept = np.polyfit(x, y, 1)

r_squared = 1 - (sum((y - (slope * x + intercept))**2) / ((len(y) - 1) * np.var(y, ddof=1)))

return r_squared

def get_r2_scipy(x, y):

_, _, r_value, _, _ = stats.linregress(x, y)

return r_value**2

def get_r2_statsmodels(x, y):

return sm.OLS(y, sm.add_constant(x)).fit().rsquared

def get_r2_python(x_list, y_list):

n = len(x)

x_bar = sum(x_list)/n

y_bar = sum(y_list)/n

x_std = math.sqrt(sum([(xi-x_bar)**2 for xi in x_list])/(n-1))

y_std = math.sqrt(sum([(yi-y_bar)**2 for yi in y_list])/(n-1))

zx = [(xi-x_bar)/x_std for xi in x_list]

zy = [(yi-y_bar)/y_std for yi in y_list]

r = sum(zxi*zyi for zxi, zyi in zip(zx, zy))/(n-1)

return r**2

def get_r2_numpy_manual(x, y):

zx = (x-np.mean(x))/np.std(x, ddof=1)

zy = (y-np.mean(y))/np.std(y, ddof=1)

r = np.sum(zx*zy)/(len(x)-1)

return r**2

def get_r2_numpy_corrcoef(x, y):

return np.corrcoef(x, y)[0, 1]**2

print('Python')

%timeit get_r2_python(x_list, y_list)

print('Numpy polyfit')

%timeit get_r2_numpy(x, y)

print('Numpy Manual')

%timeit get_r2_numpy_manual(x, y)

print('Numpy corrcoef')

%timeit get_r2_numpy_corrcoef(x, y)

print('Scipy')

%timeit get_r2_scipy(x, y)

print('Statsmodels')

%timeit get_r2_statsmodels(x, y)

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/news/370726.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

ASP .NET SVN emmet 插件

学习 ASP .NET 时间的第三周: 来讲讲如何在 visual studio 2013...上搭载 SVN吧: 废话不多说: One Step: 电脑上已安装 visual studio 2013 等版本(未安装时 红色区域是不存在的) Two Step: 从官网上下载对…

Python之路3【知识点】白话Python编码和文件操作(截载)

无意发现这篇文章讲的比较好,存下来供参考: http://www.cnblogs.com/luotianshuai/p/5735051.html转载于:https://www.cnblogs.com/shikaihong/p/7778880.html

Http协议入门

[在此处输入文章标题] 1 web web入门 1)web服务软件作用: 把本地资源共享给外部访问 2)tomcat服务器基本操作 : 启动: %tomcat%/bin/startup.bat 关闭: %tomcat%/bin/shutdown.bat 访问tomcat主页: http://loca…

计算机硬件系统都是看得见的,计算机组成硬件系统).doc

计算机组成硬件系统)各位计算机协会的成员大家好,很高兴大家能陪我们走过这段难忘的时光。为了让大家更好的学到东西,我们特地将计算机方面的东西整理成技术文档,共大家使用,祝大家学得愉快!湘信院计算机协会一&#x…

Google Guava –期货

这篇文章是我在Google Guava上的系列文章的延续,这次涵盖了Future。 Futures类是用于使用Future / ListenableFuture接口的静态实用程序方法的集合。 Future是提交给ExecutorService的异步任务(可运行或可调用)的句柄。 Future界面提供以下方…

iptables 配置后连接不上数据库_Linux服务器配置-VSFTP服务配置(三)

上文:Linux服务器配置-VSFTP服务配置(二)一、vsftpd服务防火墙配置1、主动(POST)模式 FTP 防火墙配置CentOS6 系统 iptables 的配置iptables -t filter --line-number -nL INPUT#显示现有防火墙规则,查看是否开启20、21号端口。iptables -t filter -I IN…

下标索引必须为正整数类型或逻辑类型_Python3 基本数据类型

Python中的变量不需要声明。每个变量在使用前都必须赋值,变量赋值以后该变量才会被创建。在Python中,变量就是变量,它没有类型,我们所说的"类型"是变量所指的内存中对象的类型。Python 3中有六个标准的数据类型&#xf…

noip模拟赛 写代码

分析:这其实就是括号匹配题,一眼贪心题,不过一开始贪错了,以为([)]是合法的......其实括号之间不能嵌套. 一开始的想法是尽量往左边填左括号,因为每种括号的数量都确定了,那么左括号和右括号的数量也就确定…

[CF797C] Minimal string(贪心,栈)

题目链接:http://codeforces.com/contest/797/problem/C 题意:给个字符串,求字典序最小的出栈顺序。 一开始想用优先队列记录全局最小的字符,然后每次入栈的时候检查当前字符是不是最小的,如果是那么同时pop。这样做的…

让我们紧缩大数据

作为开发人员,我们的重点是简单,有效的解决方案,因此,最有价值的原则之一就是“保持简单和愚蠢”。 但是使用Hadoop map-reduce很难坚持这一点。 如果我们要评估多个Map Reduce作业中的数据,那么最终将得到与业务无关但…

行内元素和块元素以及行内块元素的特点

一、背景 初学html&#xff0c;接触很多标签 <h1>、<p>、<span>、<ul>、<em>等&#xff0c;当写出简单的小页面的时候&#xff0c;例如仅仅是一篇带有标题的文章&#xff0c;标题 <h1>标签单独一行&#xff0c;不管后面有多大的空间&…

软件测试的功能测试和性能测试,大型软件的功能测试流程及性能测试流程

大型软件具有涉及子模块繁多、建设过程复杂、功能全面、性能具有较高要求的特点。依据ISO/IEC 9126软件产品评估标准&#xff0c;需要对软件的功能性、可靠性、可用性、效率、可维护性、可移植性等方面进行评估。因此&#xff0c;需要有一种方法能够对大型软件进行测试&#xf…

vue 分模块打包 脚手架_Vue面试官最爱的底层源码问题,你可以这样回答!

最近看到身边很多人都在投简历&#xff0c;有因为企业裁员的&#xff0c;有因为自己想跳槽的&#xff0c;原因不一&#xff0c;但是最终大家都会需要接触到面试这个事情。但是很多人对待面试不够认真&#xff0c;只会等待结果&#xff0c;不去努力。所以这边想整理一些懒人面试…

re.containerbase.startinternal 子容器启动失败_Python项目容器化实践(二) Docker Machine和Docker Swarm...

前言这篇文章介绍Docker生态中的常被提到的Engine、Machine和Swarm&#xff0c;大家以了解为主&#xff0c;工作需要再深入。EngineDocker Engine其实就是我们常说的「Docker」&#xff0c;它是一个C/S模型(Client/Server)的应用&#xff0c;包含如下组件:Daemon。守护进程&…

基于设备树的TQ2440的中断(2)

下面以按键中断为例看看基于设备数的中断的用法&#xff1a; 设备树&#xff1a; tq2440_key {compatible "tq2440,key";interrupt-parent <&gpf>;interrupts <0 IRQ_TYPE_EDGE_FALLING>, <1 IRQ_TYPE_EDGE_FALLING>;key_3 <&gpf 2…

计算机里有个不能进入的磁盘分区,新电脑只有一个分区怎么办? 教你们如何不进pe给硬盘创建新分区!...

很多朋友新电脑刚买回来打开发现明明自己机械硬盘1T或者1T机械加128G固态&#xff0c;但是却只有一个或者两个分区&#xff0c;但是又不会分区现在教大家如何不用老毛桃大白菜之类的进pe系统里面就能直接创建新分区1 WinR输入diskmgmt.msc2进入磁盘管理可以查看本机的C盘与E盘的…

OSGi中的权限

在上一篇文章中 &#xff0c;我们介绍了为Java应用程序实现沙箱的方法&#xff0c;在其中我们可以安全地运行移动代码 。 这篇文章探讨了如何在OSGi环境中执行相同的操作。 OSGi OSGi规范 为Java定义了一个动态模块系统 。 因此&#xff0c;它是实施那种可以使您的应用程序动…

HTTP简单教程

目录 HTTP简介 HTTP工作原理 HTTP消息结构 客户端请求消息服务器响应消息实例 HTTP请求方法HTTP响应头信息HTTP状态码 HTTP状态码分类HTTP状态码列表 HTTP content-type对照表 HTTP简介 HTTP协议是Hyper Text Transfer Protocol&#xff08;超文本传输协议&#xff09;的缩写&…

Reversed-Z详解

在3D渲染管线中&#xff0c;Z这个家伙几乎无处不在&#xff0c;如Z-Buffer&#xff0c;Early-Z&#xff0c;Z-Cull&#xff0c;Z-Test&#xff0c;Z-Write等等&#xff0c;稍有接触图形学的人都会对这些术语有所耳闻。 那么Z到底是什么呢&#xff1f;首先Z当然可以是任意坐标系…

pyqt开发的程序模板_小程序定制开发和模板开发要多少钱?有什么区别?

到现在&#xff0c;小程序开发已经有了1年多的历史&#xff0c;已经达到百万数量级。无论是小程序商城还是小程序游戏&#xff0c;其开发方式不外乎两种&#xff0c;一种是定制开发&#xff0c;另一种是模板开发。对于很多初次接触小程序的客户来说&#xff0c;还不知道小程序的…