Matlab解决线性回归问题论文Word格式.docx

上传人:b****8 文档编号:22626610 上传时间:2023-02-04 格式:DOCX 页数:27 大小:345.34KB
下载 相关 举报
Matlab解决线性回归问题论文Word格式.docx_第1页
第1页 / 共27页
Matlab解决线性回归问题论文Word格式.docx_第2页
第2页 / 共27页
Matlab解决线性回归问题论文Word格式.docx_第3页
第3页 / 共27页
Matlab解决线性回归问题论文Word格式.docx_第4页
第4页 / 共27页
Matlab解决线性回归问题论文Word格式.docx_第5页
第5页 / 共27页
点击查看更多>>
下载资源
资源描述

Matlab解决线性回归问题论文Word格式.docx

《Matlab解决线性回归问题论文Word格式.docx》由会员分享,可在线阅读,更多相关《Matlab解决线性回归问题论文Word格式.docx(27页珍藏版)》请在冰豆网上搜索。

Matlab解决线性回归问题论文Word格式.docx

r残差

rint残差置信区间

stats用于检验回归模型的统计量,有三个数值:

相关系数R2、F值、与F对应的概率p,相关系数R2越接近1,说明回归方程越显著;

F>

F1-α(k,n-k-1)时拒绝H0,F越大,说明回归方程越显著;

与F对应的概率p时拒绝H0,回归模型成立。

Y为n*1的矩阵;

X为(ones(n,1),x1,…,xm)的矩阵;

alpha显著性水平(缺省时为0.05)。

(二)多元线性回归

1.命令regress

2.命令rstool多元二项式回归

命令:

rstool(x,y,’model’,alpha)

x为n*m矩阵

y为n维列向量

model由下列4个模型中选择1个(用字符串输入,缺省时为线性模型):

linear(线性):

purequadratic(纯二次):

interaction(交叉):

quadratic(完全二次):

alpha显著性水平(缺省时为0.05)

返回值beta系数

返回值rmse剩余标准差

返回值residuals残差

非线性回归

1.命令nlinfit

[beta,R,J]=nlinfit(X,Y,’’model’,beta0)

X为n*m矩阵

Y为n维列向量

model为自定义函数

beta0为估计的模型系数

beta为回归系数

R为残差

2.命令nlintool

nlintool(X,Y,’model’,beta0,alpha)

alpha显著性水平(缺省时为0.05)

3.命令nlparci

betaci=nlparci(beta,R,J)

返回值为回归系数beta的置信区间

4.命令nlpredci

[Y,DELTA]=nlpredci(‘model’,X,beta,R,J)

Y为预测值

DELTA为预测值的显著性为1-alpha的置信区间;

三、多元线性回归问题

在实际经济问题中,一个变量往往受到多个变量的影响。

例如,家庭消费支出,除了受家庭可支配收入的影响外,还受诸如家庭所有的财富、物价水平、金融机构存款利息等多种因素的影响,表现在线性回归模型中的解释变量有多个。

这样的模型被称为多元线性回归模型。

(multivariablelinearregressionmodel)

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

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

Y=0+1X1i+2X2i++kXki+i,i=1,2,,n

(1)

其中k为解释变量的数目,j(j1,2,k)称为回归系(regressioncoefficient)。

式也被称为总体回归函数的随机表达式。

它的非随机表达式为:

Y=0+1X1i+2X2i++kXki,i=1,2,n

(2)

®

j也被称为偏回归系数(partialregressioncoefficient)。

2、多元线性回归计算模型

Y=0+1x1+2x2+kxk+∼N(0,2)(3)

多元性回归模型的参数估计,同一元线性回归方程一样,也是在要求误差平方和(Σe)为最小的前提下,用最小二乘法或最大似然估计法求解参数。

设(x11,x12,x1p,y1),,(xn1,xn2,xnp,yn)是一个样本,用最大似然估计法估计参数:

取b0,b1bp,当b0b0,b1b1,b−b达到最小。

例1

输入以下程序:

输出结果如图1所示:

在Matlab命令窗口输入

z=b

(1)+b

(2)*x%运算

plot(x,Y,‘k+'

x,z,'

r'

)%输出运算

得到预测比较图如图2、3所示:

图2

例2:

(多元线性回归)

水泥凝固时放出的热量y与水泥中的四种化学成分x1,x2,x3,x4有关,今测得一组数据如下,试确定多元线性模型

%数据

%运行方程

输出结果如图4所示:

图4

结果如图5所示:

图5

在matlab命令框中输入

得到如图6所示:

图6

四、线性回归拟合

对于多元线性回归模型:

设变量

的n组观测值为

,则

的估计值为

(11.2)

在Matlab中,用regress函数进行多元线性回归分析,应用方法如下:

语法:

b=regress(y,x)

[b,bint,r,rint,stats]=regress(y,x)

[b,bint,r,rint,stats]=regress(y,x,alpha)

b=regress(y,x)得到的

维列向量b即为(11.2)式给出的回归系数

的估计值.

[b,bint,r,rint,stats]=regress(y,x)给出回归系数

的估计值b,

的95%置信区间(

向量)bint,残差r以及每个残差的95%置信区间(

向量)rint;

向量stats给出回归的R2统计量和F以及临界概率p的值.

如果

的置信区间(bint的第

行)不包含0,则在显著水平为

时拒绝

的假设,认为变量

是显著的.

[b,bint,r,rint,stats]=regress(y,x,alpha)给出了bint和rint的100(1-alpha)%的置信区间.

三次样条插值函数的MATLAB程序

matlab的spline

x=0:

10;

y=sin(x);

%插值点

xx=0:

.25:

%绘图点

yy=spline(x,y,xx);

plot(x,y,'

o'

xx,yy)

结果如图7所示

图7

五、非线性拟合

非线性拟合可以用以下命令(同样适用于线形回归分析):

1.beta=nlinfit(X,y,fun,beta0)

X给定的自变量数据,Y给定的因变量数据,fun要拟合的函数模型(句柄函数或者内联函数形式),

beta0函数模型中系数估计初值,beta返回拟合后的系数

2.x=lsqcurvefit(fun,x0,xdata,ydata)

fun要拟合的目标函数,x0目标函数中的系数估计初值,xdata自变量数据,ydata函数值数据

X拟合返回的系数(拟合结果)

nlinfit

格式:

[beta,r,J]=nlinfit(x,y,’model’,beta0)

Beta估计出的回归系数

r残差

JJacobian矩阵

x,y输入数据x、y分别为n*m矩阵和n维列向量,对一元非线性回归,x为n维列向量。

model是事先用m-文件定义的非线性函数

beta0回归系数的初值

例1已知数据:

x1=[0.5,0.4,0.3,0.2,0.1];

x2=[0.3,0.5,0.2,0.4,0.6];

x3=[1.8,1.4,1.0,1.4,1.8];

y=[0.785,0.703,0.583,0.571,0.126]’;

且y与x1,x2,x3关系为多元非线性关系(只与x2,x3相关)为:

 

y=a+b*x2+c*x3+d*(x2.^2)+e*(x3.^2)

求非线性回归系数a,b,c,d,e。

(1)对回归模型建立M文件model.m如下:

functionyy=myfun(beta,x)

x1=x(:

1);

%数据

x2=x(:

2);

x3=x(:

3);

yy=beta

(1)+beta

(2)*x2+beta(3)*x3+beta(4)*(x2.^2)+beta(5)*(x3.^2);

%运算出图

(2)主程序如下:

x=[0.5,0.4,0.3,0.2,0.1;

0.3,0.5,0.2,0.4,0.6;

1.8,1.4,1.0,1.4,1.8];

y=[0.785,0.703,0.583,0.571,0.126];

beta0=[1,1,1,1,1];

[beta,r,j]=nlinfit(x,y,@myfun,beta0)%运算出图

输出结果如图8所示

图8

例题2:

混凝土的抗压强度随养护时间的延长而增加,现将一批混凝土作成12个试块,记录了养护日期(日)及抗压强度y(kg/cm2)的数据:

养护时间:

x=[234579121417212856]

抗压强度:

y=[35+r42+r47+r53+r59+r65+r68+r73+r76+r82+r86+r99+r]

建立非线性回归模型,对得到的模型和系数进行检验。

注明:

此题中的+r代表加上一个[-0.5,0.5]之间的随机数

模型为:

y=a+k1*exp(m*x)+k2*exp(-m*x);

Matlab程序:

x=[234579121417212856];

r=rand(1,12)-0.5;

y1=[354247535965687376828699];

y=y1+r;

%运算路径

myfunc=inline('

beta

(1)+beta

(2)*exp(beta(4)*x)+beta(3)*exp(-beta(4)*x)'

'

beta'

x'

);

%运算路径

beta=nlinfit(x,y,myfunc,[0.50.50.50.5]);

%运算方法

a=beta

(1),k1=beta

(2),k2=beta(3),m=beta(4)%数据

%testthemodel%运算方法

xx=min(x):

max(x);

%运算方法

yy=a+k1*exp(m*xx)+k2*exp(-m*xx);

xx,yy,'

结果如图9所示:

图9

lsqnonlin

非线性最小二乘(非线性数据拟合)的标准形式为

其中:

L为常数

在MATLAB5.x中,用函数leastsq解决这类问题,在6.0版中使用函数lsqnonlin。

则目标函数可表达为

x为向量,F(x)为函数向量。

函数lsqnonlin

格式x=lsqnonlin(fun,x0)%x0为初始解向量;

fun为

,i=1,2,…,m,fun返回向量值F,而不是平方和值,平方和隐含在算法中,fun的定义与前面相同。

x=lsqnonlin(fun,x0,lb,ub)%lb、ub定义x的下界和上界:

x=lsqnonlin(fun,x0,lb,ub,options)%options为指定优化参数,若x没有界,则lb=[],ub=[]。

[x,resnorm]=lsqnonlin(…)%resnorm=sum(fun(x).^2),即解x处目标函数值。

[x,resnorm,residual]=lsqnonlin(…)%residual=fun(x),即解x处fun的值。

[x,resnorm,residual,exitflag]=lsqnonlin(…)%exitflag为终止迭代条件。

[x,resnorm,residual,exitflag,output]=lsqnonlin(…)%output输出优化信息。

[x,resnorm,residual,exitflag,output,lambda]=lsqnonlin(…)%lambda为Lagrage乘子。

[x,resnorm,residual,exitflag,output,lambda,jacobian]=lsqnonlin(…)%fun在解x处的Jacobian矩阵。

例求下面非线性最小二乘问题

初始解向量为x0=[0.3,0.4]。

解:

先建立函数文件,并保存为myfun.m,由于lsqnonlin中的fun为向量形式而不是平方和形式,因此,myfun函数应由

建立:

k=1,2,…,10

functionF=myfun(x)

k=1:

F=2+2*k-exp(k*x

(1))-exp(k*x

(2));

然后调用优化程序:

x0=[0.30.4];

[x,resnorm]=lsqnonlin(@myfun,x0)%主运算方法

结果为:

Optimizationterminatedsuccessfully:

NormofthecurrentstepislessthanOPTIONS.TolX

x=

0.25780.2578

resnorm=%求目标函数值

lsqcurvefit

非线性曲线拟合是已知输入向量xdata和输出向量ydata,并且知道输入与输出的函数关系为ydata=F(x,xdata),但不知道系数向量x。

今进行曲线拟合,求x使得下式成立:

在MATLAB5.x中,使用函数curvefit解决这类问题。

函数lsqcurvefit

格式x=lsqcurvefit(fun,x0,xdata,ydata)%运算方法

x=lsqcurvefit(fun,x0,xdata,ydata,lb,ub)

x=lsqcurvefit(fun,x0,xdata,ydata,lb,ub,options)

[x,resnorm]=lsqcurvefit(…)%运算方法

[x,resnorm,residual]=lsqcurvefit(…)%运算方法

[x,resnorm,residual,exitflag]=lsqcurvefit(…)%运算方法

[x,resnorm,residual,exitflag,output]=lsqcurvefit(…)

[x,resnorm,residual,exitflag,output,lambda]=lsqcurvefit(…)%运算方法

[x,resnorm,residual,exitflag,output,lambda,jacobian]=lsqcurvefit(…)%运算方法

参数说明:

x0为初始解向量;

xdata,ydata为满足关系ydata=F(x,xdata)的数据;

lb、ub为解向量的下界和上界

,若没有指定界,则lb=[],ub=[];

options为指定的优化参数;

fun为拟合函数,其定义方式为:

x=lsqcurvefit(@myfun,x0,xdata,ydata),

其中myfun已定义为functionF=myfun(x,xdata)

F=…%计算x处拟合函数值fun的用法与前面相同;

resnorm=sum((fun(x,xdata)-ydata).^2),即在x处残差的平方和;

residual=fun(x,xdata)-ydata,即在x处的残差;

exitflag为终止迭代的条件;

output为输出的优化信息;

lambda为解x处的Lagrange乘子;

jacobian为解x处拟合函数fun的jacobian矩阵。

例求解如下最小二乘非线性拟合问题

已知输入向量xdata和输出向量ydata,且长度都是n,拟合函数为

即目标函数为

初始解向量为x0=[0.3,0.4,0.1]。

先建立拟合函数文件,并保存为myfun.m

functionF=myfun(x,xdata)

F=x

(1)*xdata.^2+x

(2)*sin(xdata)+x(3)*xdata.^3;

然后给出数据xdata和ydata

>

xdata=[3.67.79.34.18.62.81.37.910.05.4];

ydata=[16.5150.6263.124.7208.59.92.7163.9325.054.3];

x0=[10,10,10];

%初始估计值

[x,resnorm]=lsqcurvefit(@myfun,x0,xdata,ydata)

RelativefunctionvaluechangingbylessthanOPTIONS.TolFun

x=0.22690.33850.3021

resnorm=6.2950

对于上述这些可化为线性模型的回归问题,一般先将其化为线性模型,然后再用最小二乘法求出参数的估计值,最后再经过适当的变换,得到所求回归曲线。

在熟练掌握最小二乘法的情况下,解决上述问题的关键是确定曲线类型和怎样将其转化为线性模型。

如图10所示

图10

确定曲线类型一般从两个方面考虑:

一是根据专业知识,从理论上推导或凭经验推测、二是在专业知识无能为力的情况下,通过绘制和观测散点图确定曲线大体类型。

六、Subplot问题

程序:

x1=[111112222333344445556666788881010101011111212131314151616161720]'

;

x2=[1010010000111010000101011011000111001010110000]'

x5=[1,3,3,2,3,2,2,1,3,2,1,2,3,1,3,3,2,2,3,1,1,3,2,2,1,2,1,3,1,1,2,3,2,2,1,2,3,1,2,2,3,2,2,1,2,1]'

x3=[1000000100100100000110001010110000100100000101]'

x4=[0001011001010000110000110100001011010011011010]'

y=[13876116081870111283117672087211772105351219512313149752137119800114172026313231128841324513677159651236621352138392288416978148031740422184135481446715942231742378025410148611688224170159902633017949256852783718838174831920719346];

size(x5)

x=[ones(46,1),x1,x2,x3,x4]

[b,bint,r,rint,stats]=regress(y,x,0.05)%运算方法

plot(x1,r,'

r+'

holdon

fori=1:

1:

46

if(x2(i,1)==0&

x5(i,1)==1)

k(i)=1;

%交替运算

elseif(x2(i,1)==1&

x5(i,1)==1)

k(i)=2;

elseif(x2(i,1)==0&

x5(i,1)==2)

k(i)=3;

k(i)=4;

%交替运算

x5(i,1)==3)

k(i)=5;

else

k(i)=6;

end

end

k;

plot(k,r,'

g+'

输出如图11所示

图11

第二个

x2=[101001000011101

展开阅读全文
相关资源
猜你喜欢
相关搜索

当前位置:首页 > 总结汇报 > 学习总结

copyright@ 2008-2022 冰豆网网站版权所有

经营许可证编号:鄂ICP备2022015515号-1