草庐IT

利用python进行回归分析

first青年危机 2023-07-07 原文

通常大家会认为曲线拟合和回归分析类似,但其实回归分析中是包含曲线拟合的。拟合是研究因变量和自变量的函数关系的。而回归是研究随机变量间的相关关系的。拟合侧重于调整参数,使得与给出的数据相符合。而回归则是侧重于研究变量的关系,对拟合问题做统计分析。

一元线性回归

一元线性回归模型的一般形式

数据通常呈一条直线,则y和x之间的关系通常可以看做近似线性关系。但是一般来说这些数据点并不在一条直线上,这说明y和x的关系并没有确切到给定x就可以唯一确定y的程度。其实y还受到很多因素的影响。如果主要研究y和x的关系,可以假设有如下关系。

(1)为未知待定常数称为回归系数,是其他随机因素对y的影响,并且服从分布。

2)对于参数的估计有多重方法,一方面我们可以使用拟合的方法去估计参数,另一方面我们可以用最小二乘法进行参数的估计。

最小二乘法:使得回归模型和直线方程所在的所有数据点都比较接近。

,我们为了求得最小Q,对这个式子求导最终计算出

(3)相关性检验,拟合优度和剩余标准差:

,这是原始数据的总变异平方和。

,拟合直线可解释的变异平方和。

,残差平方和。

由公式推导可以证明有:

对于一个具体的样本来说,是一个定值。所以可解释变异越大,则必有残差平方和越小。这个式子可以从两个方面解释拟合程度的优良。

(1)越大,在回归方程中解释因变量得变异部分越大,回归方程对原始数据的解释越好。

(2)越小,观测值围绕回归直线越紧密,回归方程对源数据的拟合程度越好,

 因此定义一个测量标量来说明回归方程对原始数据的拟合程度,这个系数叫拟合优度。拟合优度是指可解释变异占总体变异的百分比,用表示。有:

有以下简单的性质:

(1)

(2)当=1时,有,此时原数据的变异程度完全可以由拟合值的变异来解释,且残差为0(),即原数据和拟合点完全吻合。

(3) 当=0时,此时的回归方程完全不能解释原数据的总变异,且残差为0(),即y得变异和x完全无关。

 ( 4)回归方程的显著检验

在以上的讨论中,我们假设y是关于x的回归方程具有形式。在实际中,需要检验是否为x的线性函数,若与x成线性函数为真,则 不为0,如果为0,则y和x就无线性关系,所以我们需要进行假设检验。 

提出假设::,

如果接受原假设,则说明回归结果不显著,否则回归结果显著。

求解可到的拒绝域为:

 如果满足拒绝域条件,则可以在一定的置信水平上认为线性回归结果显著。

可以通过查表或者利用软件查出,是检验水平。

当回归效果显著的时候:

可以对系数做置信区间估计,系数的置信区间为

 同理得到相应的置信区间,

点预测和预测区间:

假设的处的观察结果,是利用回归方程计算的观察值。在给定的置信区间中,

这意味着对于观察值,有的把握认为在区间中。

 例子:合金强度y和其中一个碳含量x有比较密切的关系,今从生产中收集了一批数据如表所示,试拟合一个函数,并且对回归模型进行检验。

x0.100.110.120.130.140.150.160.170.180.200.220.24
y42.042.545.045.545.047.549.051.050.055.057.559.5

首先我们将图像的散点图画出来,判断使用哪一类方程进行拟合,然后求出方程的未知系数。然后对方程进行拟合优度的检验,最后对系数进行拟合优度的判断。

import pylab as pl
import numpy as np
import statsmodels.api as sm
x = np.array([0.10, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.20, 0.22, 0.24])
y = np.array([42.0, 42.5, 45.0, 45.5, 45.0, 47.5, 49.0, 51.0, 50.0, 55.0, 57.5, 59.5])
def sp(a,b):
    pl.plot(a, b, 'o', color='black')
    pl.show()
def main(d):
    re=sm.formula.ols('y~x',d).fit()
    print(re.summary())
sp(x,y)
x1=x;x1=np.delete(x1,[8,4]);y1=y;y1=np.delete(y1,[8,4])
d={'x':x1,'y':y1}
main(d)

(1)这里首先介绍一下基于公式构建并拟合模型的调用公式,调用格式为:

import statsmodels.api as sm
sm.formula.ols(formula,data=df).fit()

这当中formula是引号括起来的公式,df为数据框或字典格式的数据。使用方法summary获取一个关于此回归的统计检验。

(2)我们在进行数据拟合时通常需要对异常值进行剔除,这里是利用t检验将样本中的异常值进行了剔除:

print(re.outlier_test())

注意:通常使用这个函数,样本数量最好大于20个。

多元线性回归

当自变量的个数多余一个时,回归模型是关于自变量的线性表达形式,则称此模型为多元线性回归,数学模型可以表达为:

 式中:为回归系数;为随机误差,服从正态分布,未知。

回归分析的主要步骤是:

(1)由观测值确定参数(回归系数)的估计值

(2)对线性关系,自变量的显著性进行统计检验

(3)利用回归方程进行预测。

(1) 确定参数

对于模型式中的参数用最小二乘法估计,即选取,使得误差平方和最小。原理和上面一元线性回归模型类似,在这里不过多介绍。

(2) 统计检验

前面是在假设随机变量y和变量j具有线性关系的条件下建立线性回归方程,但是变量y与变量是否为线性关系?所有的变量对变量y是否有影响?都需要做统计检验。

对于总平方,有,是定义的残差平方和,反映随机误差对y的影响;称回归平方和,反映自变量对y的影响。这里是因变量的均值。,其中SST的自由度,SSE的自由度,SSR的自由度

因变量y与自变量之间是否存在线性规关系是需要检验的,显然如果所有都很小,y和的线性关系就不明显,所以我们可令原假设为

:

成立时有分解式定义的SSE,SSR满足

~

在显著性水平上,有上分位数,若回归方程效果显著,若,回归效果不显著。

注:

1.当y和的关系不明显时,可能存在非线性关系,例如平方关系。

2.这里的F分布检验直接针对这个模型的所有系数,也就说当这个模型的回归效果显著时,依旧有可能某个自变量与y是不相关的,因此我们还需要对各个回归系数做假设检验,这里使用的是t检验。

上面在一元回归中已经有过定义在这里就不在详述。

(3)回归方程的预测

对于需要预测的自变量,我们直接将其代入我们建立的模型当中,然后就得到了预测值。

我们也可以进行区间估计。记

,则的置信度为的预测区间可以近似为这是样本量较大时的预测区间。

statsmodels可以使用两种模式求解回归分析模型:一种是基于公式的模型,一种是基于数组的模型。

基于公式的模型:我们在上面的一元线性回归中已经介绍过。

基于数组的模型:

import statsmodels.api as sm
sm.OLS(y,X).fit()

其中y是因变量的观察值向量,而X为自变量观测值矩阵在添加第一列全部为1的增广矩阵。

例题:某企业为了研究车工的平均年龄,平均文化程度和平均产量之间的变化关系,随机抽取了8个班组,测得的数据如下表所示,试求平均产量对平均工龄,平均文化程度的回归方程。并求时,y的预测值。

12345678
x17.16.89.211.48.76.610.310.6
x211.110.812.410.99.69.010.512.4
x315.41522.827.819.513.124.926.2
import pylab as pl
import numpy as np
import statsmodels.api as sm
a=np.loadtxt('data10_2.txt')
pl.subplot(121)
pl.plot(a[0],a[2],'o',label='x1',color='black')
pl.plot(a[1],a[2],'*',label='x2',color='black')
pl.legend()
d={'x1':a[0], 'x2':a[1],'y':a[2]}
re=sm.formula.ols('y~x1+x2',d).fit()
print(re.summary())
print(re.mse_resid)
yh=re.predict({'x1':[9,10],'x2':[10,9]})
print(yh)
pre=re.get_prediction(d)
df=pre.summary_frame(alpha=0.05)
dfv=df.values;low,upp=dfv[:,4:].T
r=(upp-low)/2
num=np.arange(1,len(a[0])+1)
pl.subplot(122)
pl.errorbar(num,re.resid,r,fmt='o')
pl.show()

这里是基于公式的代码,不做详细介绍。下面附上基于数组的代码,我个人比较推荐基于公式的算法。

import pylab as pl
import numpy as np
import statsmodels.api as sm
a=np.loadtxt('data10_2.txt')
pl.plot(a[0],a[2],'o',label='x1',color='black')
pl.plot(a[1],a[2],'*',label='x2',color='black')
pl.legend()
X=sm.add_constant(a[:2].T)
re=sm.OLS(a[2],X).fit()
print(re.summary())
print(re.mse_resid)
yh=re.predict(np.array([[1,9,10],[1,10,9]]))
print(yh)
pl.show()

(1)这里的函数add_comstant(x)作用是为矩阵x添加一列全为一的作用。

(2)对于基于数组的函数预测和基于公式的有所不同,它的格式和填入OLS函数的格式相同,第一列全为一,第二列为x1,一次类推。

多项式回归

如果我们做出散点图发现y和x呈较明显的二次(或者高次)函数关系,或者用线性关系模拟的效果不太好,就可以使用多项式回归。随机意义下,一元n次多项式回归的数学模型可以表达为

式中:为随机误差,满足,

例子:某厂生产一种电器的销量y和竞争对手的价格和本场的价格有关,我们给出下表,是此商品在10个城市的销售情况,试根据这些数据建立y和的关系式,并对得到的模型和系数进行检验。并预测若本厂商品在本市售价为160元,那么竞争对手售价为170元。预测这个商品在本市销量。

120140190130155175125145180150
10011090150210150250270300250
10210012077469326696585

我们画出散点图,从中可以看出,y和有线性关系,而y和之间的关系难以确定,所以我们多拟合几种模型,去取残差方差最小的这个模型为最佳模型。

线性模型:

纯二次模型:

交叉二次:

纯二次:

import pylab as pl
import numpy as np
import statsmodels.api as sm
a=np.loadtxt('data10_4.txt')
x1=a[0];x2=a[1];y=a[2];
pl.plot(x1,y,'o',label='x1',color='black')
pl.plot(x2,y,'*',label='x2',)
pl.legend()
pl.show()
d={'x1':x1, 'y1':y, 'x2':x2}
resid=[]
idea=['y~x1+x2','y~x1+x2+I(x1**2)','y~x1*x2','y~x1*x2+I(x1**2)+I(x2**2)']
for i in range(0,4):
    re=sm.formula.ols(idea[i],d).fit()
    resid.append(re.mse_resid)
i=resid.index(max(resid))
re=re=sm.formula.ols(idea[i],d).fit()
print('最佳预测值:',re.predict({'x1':170,'x2':160}))

通过选择残差方差这个准则,确定纯二次多项式模型式是拟合的最好的。

逐步回归

实际中影响因变量的因素可能很多,有的可能关联性很强,而有的可能弱一些,人们总希望从中挑选对因变量影响显著的自变量来建立回归模型,下面来讨论多元线性回归模型的情形。

基本思想:为候选自变量的集合,是集合S的一个子集。设中有l个自变量,由和因变量y构造的回归模型的残差平方和为,则模型的残差方差,n为数据样本容量。所选子集应该使尽可能的小一般来说回归模型所含的自变量越多,那么残差平方和就越小,如果变量中存在对因变量影响很小的量,那么并不会影响,但是会增大导致分母变大,使得残差方差增大。

当自变量较多时,我们有自主选择自变量的方法:

通常我们最常用的是逐步回归法。

前进法:逐步添加变量,每步添加一个,直到没有变量可以引入为止。

第一步:对所有的变量和因变量做一元线性回归,然后得到每个回归模型的F统计量。如果,则将对应变量引入方程中。

第二步:假设我们第一步共引入了m个变量,然后把这m个变量再逐一和剩下的变量建立回归方程,并将其F值和做对比。

第三步:知道所有引入后所有变量的F值均小于。(P为引入变量的个数) 

重复上述步骤,知道添加完所有变量。

后退法:使用与F统计量等价的T检验量

第一步:用所有自变量和因变量做拟合。

第二步:每一步都把未通过t检验的自变量中选择最小一个删除,如果均通过t检验,可以为了使得最小而删除,最小的变量。

第三步:直到所有的自变量均通过t检验。

逐步回归法:

把变量一个一个的引入,每次引入一个变量就对已选入的变量进行逐个检验,当原引入的变量因为后引入的变量而不再显著时,需要进行剔除。引入一个变量或者删除一个变量都要保证模型中的所有变量均为显著的。直到既无显著变量引入,也无不显著的变量从模型里删除。

要求:引入自变量的显著性水平要高于剔除变量的显著性水平。

Logistic回归模型

通常我们会面临二分类问题,比方说我们去医院抽血化验,我们会得到一张化验单上面会标注血液中各成分的含量,医生需要判断我们的身体到底有病或者没病。

这个时候就是一个二分类问题,而Logistic通常给出的是一个概率值,如果有病1,为没病为0的话,我们对一个样本进行预测会得到一个在0~1之间的预测值,如果这个值大于0.5我们就认为他有病,反之则无。

调用格式为:

import statsmodels.api as sm
md=sm.formula.logit('y~x1+x2+x3',d).fit()

调用格式基本和前面的线性回归类似。

有关利用python进行回归分析的更多相关文章

  1. ruby-on-rails - 使用 Ruby on Rails 进行自动化测试 - 最佳实践 - 2

    很好奇,就使用ruby​​onrails自动化单元测试而言,你们正在做什么?您是否创建了一个脚本来在cron中运行rake作业并将结果邮寄给您?git中的预提交Hook?只是手动调用?我完全理解测试,但想知道在错误发生之前捕获错误的最佳实践是什么。让我们理所当然地认为测试本身是完美无缺的,并且可以正常工作。下一步是什么以确保他们在正确的时间将可能有害的结果传达给您? 最佳答案 不确定您到底想听什么,但是有几个级别的自动代码库控制:在处理某项功能时,您可以使用类似autotest的内容获得关于哪些有效,哪些无效的即时反馈。要确保您的提

  2. python - 如何使用 Ruby 或 Python 创建一系列高音调和低音调的蜂鸣声? - 2

    关闭。这个问题是opinion-based.它目前不接受答案。想要改进这个问题?更新问题,以便editingthispost可以用事实和引用来回答它.关闭4年前。Improvethisquestion我想在固定时间创建一系列低音和高音调的哔哔声。例如:在150毫秒时发出高音调的蜂鸣声在151毫秒时发出低音调的蜂鸣声200毫秒时发出低音调的蜂鸣声250毫秒的高音调蜂鸣声有没有办法在Ruby或Python中做到这一点?我真的不在乎输出编码是什么(.wav、.mp3、.ogg等等),但我确实想创建一个输出文件。

  3. ruby-on-rails - 按天对 Mongoid 对象进行分组 - 2

    在控制台中反复尝试之后,我想到了这种方法,可以按发生日期对类似activerecord的(Mongoid)对象进行分组。我不确定这是完成此任务的最佳方法,但它确实有效。有没有人有更好的建议,或者这是一个很好的方法?#eventsisanarrayofactiverecord-likeobjectsthatincludeatimeattributeevents.map{|event|#converteventsarrayintoanarrayofhasheswiththedayofthemonthandtheevent{:number=>event.time.day,:event=>ev

  4. ruby - 使用 C 扩展开发 ruby​​gem 时,如何使用 Rspec 在本地进行测试? - 2

    我正在编写一个包含C扩展的gem。通常当我写一个gem时,我会遵循TDD的过程,我会写一个失败的规范,然后处理代码直到它通过,等等......在“ext/mygem/mygem.c”中我的C扩展和在gemspec的“扩展”中配置的有效extconf.rb,如何运行我的规范并仍然加载我的C扩展?当我更改C代码时,我需要采取哪些步骤来重新编译代码?这可能是个愚蠢的问题,但是从我的gem的开发源代码树中输入“bundleinstall”不会构建任何native扩展。当我手动运行rubyext/mygem/extconf.rb时,我确实得到了一个Makefile(在整个项目的根目录中),然后当

  5. ruby - 如何进行排列以有效地定制输出 - 2

    这是一道面试题,我没有答对,但还是很好奇怎么解。你有N个人的大家庭,分别是1,2,3,...,N岁。你想给你的大家庭拍张照片。所有的家庭成员都排成一排。“我是家里的friend,建议家庭成员安排如下:”1岁的家庭成员坐在这一排的最左边。每两个坐在一起的家庭成员的年龄相差不得超过2岁。输入:整数N,1≤N≤55。输出:摄影师可以拍摄的照片数量。示例->输入:4,输出:4符合条件的数组:[1,2,3,4][1,2,4,3][1,3,2,4][1,3,4,2]另一个例子:输入:5输出:6符合条件的数组:[1,2,3,4,5][1,2,3,5,4][1,2,4,3,5][1,2,4,5,3][

  6. ruby - 即使失败也继续进行多主机测试 - 2

    我已经构建了一些serverspec代码来在多个主机上运行一组测试。问题是当任何测试失败时,测试会在当前主机停止。即使测试失败,我也希望它继续在所有主机上运行。Rakefile:namespace:specdotask:all=>hosts.map{|h|'spec:'+h.split('.')[0]}hosts.eachdo|host|begindesc"Runserverspecto#{host}"RSpec::Core::RakeTask.new(host)do|t|ENV['TARGET_HOST']=hostt.pattern="spec/cfengine3/*_spec.r

  7. ruby - 是否可以覆盖 gemfile 进行本地开发? - 2

    我们的git存储库中目前有一个Gemfile。但是,有一个gem我只在我的环境中本地使用(我的团队不使用它)。为了使用它,我必须将它添加到我们的Gemfile中,但每次我checkout到我们的master/dev主分支时,由于与跟踪的gemfile冲突,我必须删除它。我想要的是类似Gemfile.local的东西,它将继承从Gemfile导入的gems,但也允许在那里导入新的gems以供使用只有我的机器。此文件将在.gitignore中被忽略。这可能吗? 最佳答案 设置BUNDLE_GEMFILE环境变量:BUNDLE_GEMFI

  8. ruby - 在 Windows 机器上使用 Ruby 进行开发是否会适得其反? - 2

    这似乎非常适得其反,因为太多的gem会在window上破裂。我一直在处理很多mysql和ruby​​-mysqlgem问题(gem本身发生段错误,一个名为UnixSocket的类显然在Windows机器上不能正常工作,等等)。我只是在浪费时间吗?我应该转向不同的脚本语言吗? 最佳答案 我在Windows上使用Ruby的经验很少,但是当我开始使用Ruby时,我是在Windows上,我的总体印象是它不是Windows原生系统。因此,在主要使用Windows多年之后,开始使用Ruby促使我切换回原来的系统Unix,这次是Linux。Rub

  9. Python 相当于 Perl/Ruby ||= - 2

    这个问题在这里已经有了答案:关闭10年前。PossibleDuplicate:Pythonconditionalassignmentoperator对于这样一个简单的问题表示歉意,但是谷歌搜索||=并不是很有帮助;)Python中是否有与Ruby和Perl中的||=语句等效的语句?例如:foo="hey"foo||="what"#assignfooifit'sundefined#fooisstill"hey"bar||="yeah"#baris"yeah"另外,类似这样的东西的通用术语是什么?条件分配是我的第一个猜测,但Wikipediapage跟我想的不太一样。

  10. java - 什么相当于 ruby​​ 的 rack 或 python 的 Java wsgi? - 2

    什么是ruby​​的rack或python的Java的wsgi?还有一个路由库。 最佳答案 来自Python标准PEP333:Bycontrast,althoughJavahasjustasmanywebapplicationframeworksavailable,Java's"servlet"APImakesitpossibleforapplicationswrittenwithanyJavawebapplicationframeworktoruninanywebserverthatsupportstheservletAPI.ht

随机推荐