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


猜你喜欢
- 背景是:在实际开发中,可能会遇到网络问题或者查询量比较大的情况,上一个请求还没有完成,用户就发起了下一个请求。会造成什么情况呢?但是同一个请
- 概览(循环方式 - 常用)formapforEachfilter声明遍历的数组和异步方法声明一个数组:??const skills = [&
- <?php /********************************************** *&n
- 因为在做一个项目需要筛选掉一部分产品列表中的产品,使其在列表显示时排在最后,但是所有产品都要按照更新时间排序。研究了一下系统的数据库结构后,
- 进入PyCharm后,点击File→Open,然后在弹窗中选择需要导入项目的文件夹;打开了python项目后,需要配置该项目对应的pytho
- 一、概述:Oracle的注册就是将数据库作为一个服务注册到监听程序。客户端不需要知道数据库名和实例名,只需要知道该数据库对外提供的服务名就可
- 本文实例讲述了Python素数检测的方法。分享给大家供大家参考。具体如下:因子检测:检测因子,时间复杂度O(n^(1/2))def is_p
- webpack代码拆分webpack有两种组织模块依赖的方式,同步和异步。异步依赖作为分割点,形成一个新的块,在优化了依赖树后,每一个异步区
- demo: <div v-for="item in temps" :key="
- 目录利用python反转图片/视频安装库反转效果实现代码项目地址利用python反转图片/视频准备:一张图片/一段视频python库:Pil
- QQ和微信这两款都是非常受人喜欢的聊天交友软件!可能大家平时没有留意到,也或者是大家可能很少用微信,或者很少用QQ吧!所以可能没有留意这些小
- 日一二三四五六'.split('') ['日','一','二
- 首先要用designer设计ui界面打开后就和c#一样拖动控件做ui界面保存后是xxx.ui文件再添加个工具Arguments:-m PyQ
- 1.分析 我们在用 php 制作网站时,分类是很重要的,在分类下面又再分类这第二个分类称为次分类,而现在大多
- 在网上看到一个小需求,需要用正则表达式来处理。原需求如下:找出文本中包含”因为……所以”的句子,并以两个词为中心对齐输出前后3个字,中间全输
- 在做视觉设计时,如何高效地使用图标是一门学问:该使用什么样的图标?图标该放在哪里?大小如何?图标的使用是否帮助用户更好更快的理解内容,亦或是
- 演示:<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//
- 今天星期天,因数据库太慢,最后决定将数据库进行重新整理. (假定数据库名称为:DB_ste) 1、根据现在的数据库的脚本创建一个脚本文件(F
- <!DOCTYPE html><html lang="en"><head> &nbs
- 前言在尝试将结构体序列化为 Json 时,你可能会遇到 “omitempty” 标记,本小记就来浅看