Python 普通最小二乘法(OLS)进行多项式拟合的方法
作者:薛定谔的DBA 发布时间:2021-03-09 23:06:58
标签:Python,最小二乘法,多项式,拟合
多元函数拟合。如 电视机和收音机价格多销售额的影响,此时自变量有两个。
python 解法:
import numpy as np
import pandas as pd
#import statsmodels.api as sm #方法一
import statsmodels.formula.api as smf #方法二
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
df = pd.read_csv('http://www-bcf.usc.edu/~gareth/ISL/Advertising.csv', index_col=0)
X = df[['TV', 'radio']]
y = df['sales']
#est = sm.OLS(y, sm.add_constant(X)).fit() #方法一
est = smf.ols(formula='sales ~ TV + radio', data=df).fit() #方法二
y_pred = est.predict(X)
df['sales_pred'] = y_pred
print(df)
print(est.summary()) #回归结果
print(est.params) #系数
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d') #ax = Axes3D(fig)
ax.scatter(X['TV'], X['radio'], y, c='b', marker='o')
ax.scatter(X['TV'], X['radio'], y_pred, c='r', marker='+')
ax.set_xlabel('X Label')
ax.set_ylabel('Y Label')
ax.set_zlabel('Z Label')
plt.show()
拟合的各项评估结果和参数都打印出来了,其中结果函数为:
f(sales) = β0 + β1*[TV] + β2*[radio]
f(sales) = 2.9211 + 0.0458 * [TV] + 0.188 * [radio]
图中,sales 方向上,蓝色点为原 sales 实际值,红色点为拟合函数计算出来的值。其实误差并不大,部分数据如下。
同样可拟合一元函数;
import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
df = pd.read_csv('http://www-bcf.usc.edu/~gareth/ISL/Advertising.csv', index_col=0)
X = df['TV']
y = df['sales']
est = smf.ols(formula='sales ~ TV ', data=df).fit()
y_pred = est.predict(X)
print(est.summary())
fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(X, y, c='b')
ax.plot(X, y_pred, c='r')
plt.show()
Ridge Regression:(岭回归交叉验证)
岭回归(ridge regression, Tikhonov regularization)是一种专用于共线性数据分析的有偏估计回归方法,实质上是一种改良的最小二乘估计法,通过放弃最小二乘法的无偏性,以损失部分信息、降低精度为代价获得回归系数更为符合实际、更可靠的回归方法,对病态数据的拟合要强于最小二乘法。通常岭回归方程的R平方值会稍低于普通回归分析,但回归系数的显著性往往明显高于普通回归,在存在共线性问题和病态数据偏多的研究中有较大的实用价值。
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn import linear_model
from mpl_toolkits.mplot3d import Axes3D
df = pd.read_csv('http://www-bcf.usc.edu/~gareth/ISL/Advertising.csv', index_col=0)
X = np.asarray(df[['TV', 'radio']])
y = np.asarray(df['sales'])
clf = linear_model.RidgeCV(alphas=[i+1 for i in np.arange(200.0)]).fit(X, y)
y_pred = clf.predict(X)
df['sales_pred'] = y_pred
print(df)
print("alpha=%s, 常数=%.2f, 系数=%s" % (clf.alpha_ ,clf.intercept_,clf.coef_))
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(df['TV'], df['radio'], y, c='b', marker='o')
ax.scatter(df['TV'], df['radio'], y_pred, c='r', marker='+')
ax.set_xlabel('TV')
ax.set_ylabel('radio')
ax.set_zlabel('sales')
plt.show()
输出结果:alpha=150.0, 常数=2.94, 系数=[ 0.04575621 0.18735312]
来源:https://blog.csdn.net/kk185800961/article/details/79220724
0
投稿
猜你喜欢
- 现在假如要写一个按照"标题",'内容','作者'等等进行针对性的选择,这时需要涉及到使用
- 有时需要获取远程网站的某些信息,而服务器又限制了GET方式,只能通过POST数据提交,这个时候我们可以通过asp来实现模拟提交post数据,
- 问题描述使用pandas库的read_excel()方法读取外部excel文件报错, 截图如下好像是缺少了什么方法的样子问题分析分析个啥,
- 这篇论坛文章着重介绍了SQL Server数据库简体繁体数据混用的问题,详细内容请参考下文:我现在要说的是一个在简体繁体数据混用的时候很容易
- 关于作者 王丹丹 , IBM 中国系统与技术中心软件工程师,自从 2006 年加入 IBM,一直从事 Web 系统设计和开发工作,有五年 P
- 防止Application对象在多线程访问中出现错误asp代码处理代码如下(VB):<%Application.Lock()Appli
- Python来进行查询和替换一个文本字符串?可以使用sub()方法来进行查询和替换,sub方法的格式为:sub(replacement, s
- 今天终于抽出时间瞟了一眼JavaScript的新标准ECMA262v5,让我欣喜的是,不像是因为过于激进而被废除的v4,这个新标准对Java
- 前言通常我们在一个站站点进行采集的时候,如果是小站的话 我们使用scrapy本身就可以满足。但是如果在面对一些比较大型的站点的时候,单个sc
- 分享一下 IntelliJ IDEA 2021.1 的激活破解教程,相当于永久激活了,亲测有效,下面是详细文档哦~申明:本教程 Intell
- tensor计算三通道均值今天用pytorch处理图像时,涉及到了计算均值的问题,整理一下解决思路。第一种思路tensor转换为numpy再
- 前言词云的使用相信大家已经不陌生了,使用很简单,直接调用wordcloud包就可以了。它的主要功能是根据文本词汇和词汇频率生成图片,从中可以
- WinHttp; // Microsoft WinHTTP Services, version 5.1Alias HTTPREQUEST_P
- 获取评论贴的请求头与表单数据下一篇在这里这里,我们随便选取一个网站,获取该贴评论后的请求头,表单数据以及评论贴链接。(因为涉及敏感信息,自己
- python结构体数组在C语言中我们可以通过struct关键字定义结构类型,结构中的字段占据连续的内存空间,每个结构体占用的内存大小都相同,
- 前言看到这篇文章我就默认你已经在你的电脑上使用 pipenv搭建好了虚拟环境并且设置好了开发环境(pycharm)。如果没有,请参照这篇文章
- 简单说明:思路:从数据岛menuXML中读取数据,从树的根节点开始分析树,利用 hasChildNodes() [方法:是否含有子节点 ]
- <?php// 使用Memache 作为进程锁 class lock_processlock{// key 的前缀protected
- 1 The syntax of the SQL statement is verified.SQL的语法检查2 The data dicti
- 如何修改NT的登录密码? 代码见下:<%Sub ChangeUserPassword(C