matlab回归分析结果输出,科学网—回归分析的MATLAB和R程序实现 - 王福昌的博文...
前面博客中已经讲过MATLAB中常用的命令拟合polyfit() , lsqcurvefit() ,nlinfit() 和 cftool等,这里简单介绍简单的回归分析的MATLAB和R语言实现。
例15.1
测得某种物质在不同温度下吸附另一种物质的重量如下表
$\begin{array}{c|ccccccccc}
x_i(^oC)&1.5& 1.8&2.4&3.0 &3.5 &3.9 &4.4 &4.8 &5.0 \\\hline
y_i(mg)&4.8& 5.7&7.0&8.3&10.9&12.4&13.1&13.6&15.1
\end{array}$
由所给定的样本观测值,如果画出散点图,可以看出 $9$ 个点近乎在一条直线上。因此,我们可以假设(严格地讲,应该先检验这假设)吸附量 Y 与温度 $x$ 具有线性关系:$Y = a+b,x+varepsilon$
,求 $Y$ 对 $x$ 的线性回归。
按照书上公式,在MATLAB中输入:
x = [1.5 , 1.8, 2.4, 3.0, 3.5, 3.9, 4.4, 4.8,5.0];
y = [4.8, 5.7, 7.0, 8.3, 10.9, 12.4, 13.1, 13.6,15.1];
sx = sum(x); sy = sum(y); sxy = sum(x.*y); sxx = sum(x.^2);
fprintf('sx = %6.2f, sy = %6.2f, sxy = %6.2f, sxx = %6.2fn', sx,sy,sxy,sxx)
mx =1/9* sum(x); my = 1/9*sum(y);
fprintf('mx = %6.4f, my = %6.4fn', mx, my)
lxy = sxy - 9*mx*my; lxx = sxx - 9*mx^2;
bhat = lxy/lxx; ahat = my - bhat*mx;
fprintf('ahat = %6.4f, bhat = %6.4fn', ahat, bhat)
即可得到
sx = 30.30, sy = 90.90, sxy = 344.09, sxx = 115.11
mx = 3.3667, my = 10.1000
ahat = 0.3187, bhat = 2.9053
R语言的程序代码为:
> x
> y
> slt
> summary(slt)
运行后得到
Call:
lm(formula = y ~ x)
Residuals:
Min 1Q Median 3Q Max
-0.7347 -0.2915 0.1233 0.2546 0.7505
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.3187 0.5151 0.619 0.556
x 2.9053 0.1440 20.170 1.84e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.5213 on 7 degrees of freedom
Multiple R-squared: 0.9831, Adjusted R-squared: 0.9807
F-statistic: 406.8 on 1 and 7 DF, p-value: 1.844e-07
显示参数的置信区间
> confint(slt)
2.5 % 97.5 %
(Intercept) -0.8994414 1.536795
x 2.5647360 3.245951
> x0
> lm.pred
> lm.pred
fit lwr upr
1 93.28967 83.45075 103.1286
利用xtable包可以输出为 LaTeX 代码(见:使用 xtable 将 R 表格输出 LaTeX url{http://elegantlatex.org/2015/03/31/xtable-export-r-to-latex/}),
源代码为
% latex table generated in R 3.2.3 by xtable 1.8-2 package
% Tue Mar 08 14:47:45 2016
begin{table}[ht]
centering
caption{linear model regression}
begin{tabular}{rrrrr}
toprule
& Estimate & Std. Error & t value & Pr($>$$|$t$|$) \
midrule
(Intercept) & 0.3187 & 0.5151 & 0.62 & 0.5558 \
x & 2.9053 & 0.1440 & 20.17 & 0.0000 \
bottomrule
end{tabular}
end{table}
编译结果如下
注意:
lm( )为线性模型程序包,它可以完成参数估计、假设检验。本程序中将 lm(y~x) 结果存储于slt, 再用summary(), confint( ) 和 ~predict( ) 等函数显示或调用它。
例15.5
在某种产品表面进行腐蚀刻线试验,得到腐蚀深度~$Y$ 与腐蚀时间~$x$对应的一组数据:
$\begin{array}{c|ccccccccccc}
x_i(s)& 5& 10&15&20 &30 &40 &50 &60 &70&90&120 \\\hline
y_i(\mu g)&6&10&10&13&16&17&19&23&25&29&46
\end{array}$
(1)预测腐蚀时间为 $75s$ 时,腐蚀深度的范围($1-alpha = 95%$ );
(2)若要求腐蚀深度在 $10sim20mu m$ 之间,问腐蚀时间应如何控制?
x = [5 10 15 20 30 40 50 60 70 90 120];
y = [6 10 10 13 16 17 19 23 25 29 46];
n = length(x); lxx = sum(x.^2) - n*mean(x)^2; lxy = sum(x.*y)-n*mean(x)*mean(y);
lyy = sum(y.^2) - n*mean(y)^2; bhat = lxy/lxx; ahat = mean(y) - bhat*mean(x);
fprintf('ahat = %6.4f, bhat = %6.4fn', ahat, bhat)
disp(strcat('回归方程 y = ',num2str(ahat),' + ',' ',num2str(bhat),' x'))
%(1) x0 = 75
x0 = 75; y0hat = ahat + bhat*x0;sig2hat = 1/(n-2)*(lyy-lxy^2/lxx);sighat = sqrt(sig2hat)
deltax0 = sighat*tinv(0.975,n-2)*sqrt(1+1/n+(x0-mean(x))^2/lxx)
得到
ahat = 5.3444, bhat = 0.3043
回归方程 y =5.3444 +0.30434 x
sighat =
2.2356
deltax0 =
5.4315
x0 = 75; y0hat = ahat + bhat*x0;sig2hat = 1/(n-2)*(lyy-lxy^2/lxx);sighat = sqrt(sig2hat)
deltax0 = sighat*tinv(0.975,n-2)*sqrt(1+1/n+(x0-mean(x))^2/lxx)
fprintf('腐蚀时间为75s时,腐蚀深度Y0的预测区间为[%6.4f, %6.4f]n',y0hat-deltax0,y0hat+deltax0)
%(2) 当要求腐蚀深度在10~20 mu m 之间时,近似地有
x1 = 1/bhat*(10+sighat*norminv(0.975)-ahat); x2 = 1/bhat*(20-sighat*norminv(0.975)-ahat);
fprintf('即腐蚀时间应控制在%6.4f s到%6.4f s之间,就能得到10~20 之间的腐蚀深度。n',x1,x2)
得到
sighat =
2.2356
deltax0 =
5.4315
腐蚀时间为75s时,腐蚀深度Y0的预测区间为[22.7381, 33.6011]
即腐蚀时间应控制在29.6950 s到33.7584 s之间,就能得到10~20 之间的腐蚀深度。
例15.6
出钢时所用盛钢水的钢包,由于钢水对耐火材料的侵蚀,容积不断扩大,我们希望找出使用次数 $x$ 与增大容积 $Y$之间的关系。试验数据列于下表:
$\begin{array}{c|cccccccc}
\hline
x_i&2& 3&4&5&6 &7 &8 &9\\
y_i&6.42& 8.20&9.58&9.50&9.70&10.00&9.93&9.99\\ \hline
\hline
x_i&10&11&12&13&14&15&16& \\
y_i&10.49&10.59&10.60&10.80&10.60&10.90&11.76&\\
\hline
\end{array}$
x = 2:16; y = [6.42,8.20,9.58,9.50,9.70,10.00,9.93,...
9.99,10.49,10.59,10.60,10.80,10.60,10.90,11.76];
z = 1./y'; t = 1./x'; n = length(x);
[abhat] = [ones(n,1) t]z;
fprintf('ahat = %6.4f, bhat = %6.4fn', abhat(1), abhat(2))
disp(strcat('回归方程 z = ',num2str(abhat(1)),' + ',' ',num2str(abhat(2)),' t'))
得到
ahat = 0.0812, bhat = 0.1349
回归方程 z =0.081193 +0.13491 t
R 语言程序如下
> x
> x1
> y1
> slt
> summary(slt)
得到
Call:
lm(formula = y1 ~ x1)
Residuals:
Min 1Q Median 3Q Max
-0.0105351 -0.0017476 0.0009717 0.0022768 0.0071176
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.081193 0.001915 42.40 2.52e-15 ***
x1 0.134905 0.009701 13.91 3.50e-09 ***
---
Signif. codes:0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.004408 on 13 degrees of freedom
Multiple R-squared: 0.937, Adjusted R-squared: 0.9322
F-statistic: 193.4 on 1 and 13 DF, p-value: 3.5e-09
> confint(slt)
2.5 % 97.5 %
(Intercept) 0.07705633 0.08532936
x1 0.11394772 0.15586327
非线性回归nlinfit()
nlinfit( ) 函数可以做一般的非线性回归,其调用格式为
[beta,r,J]=nlinfit(x,y,’model’, beta0)
其中beta表示估计出的回归系数,r表示残差,J表示Jacobian矩阵,x,y表示输入数据,x、y分别为矩阵和n维列向量,对一元非线性回归,x为n维列向量,model表示是事先定义的非线性函数,beta0表示回归系数的初值.
相关命令还有产生交互界面的nlintool(x,y,’model’, beta0,alpha)、非线性模型置信区间预测的[Y,DELTA]=nlpredci(’model’, x,beta,r,J),表示nlinfit 或nlintool所得的回归函数在x处的预测值Y及预测值的显著性为1-alpha的置信区间Y±DELTA.
在 例15.6中,也可以直接用MATLAB解算过程如下,输入
>>x=2:16;
y=[6.42 8.20 9.58 9.5 9.7 10 9.93 9.99 10.49 10.59 10.60 10.80 10.60 10.90 10.76];
Fit1=@(beta,x)(x./(beta(1)*x+beta(2)));% 非线性模型1
[beta,r ,J]=nlinfit(x,y,Fit1,[0.1,0.1])
得到
beta =
0.084463 0.11524
r =
-0.61815
0.061728
…………
-0.14924
J = -49.535 -24.768
-66.231 -22.077
…………
-119.01 -7.4382
即回归模型为$y = dfrac{x}{0.084463x+0.11524}$,下面预测并作图,输入
>> [yh,delta]=nlpredci(Fit1,x,beta,r ,J); plot(x,y,'k+',x,yh,'r')
得到
图15.13 散点图及非线性拟合曲线
若用例15.7模型,则继续输入
>> Fit2=@(beta,x) (beta(1)*exp(beta(2)./x)); [beta,r ,J]=nlinfit(x,y,Fit2,[11,8])
得到
beta = 11.604 -1.0641
r = -0.39589
0.061455
…………
-0.097033
J =
0.58739 3.4079
0.70138 2.7128
…………
0.93566 0.67856
即拟合模型为 $y = 11.604 e^{-1.641/x}$.
>> [yhh,delta]=nlpredci(Fit2,x,beta,r ,J); plot(x,y,'k+',x,yhh,'r')
可得到类似于图15.13 的散点图和新的非线性拟合曲线。
转载本文请联系原作者获取授权,同时请注明本文来自王福昌科学网博客。
链接地址:http://blog.sciencenet.cn/blog-292361-1019140.html
上一篇:工程硕士《高等工程数学》期末考试近三年试卷
下一篇:方差分析的MATLAB和R程序实现
matlab回归分析结果输出,科学网—回归分析的MATLAB和R程序实现 - 王福昌的博文...相关推荐
- 单纯性搜索算法 matlab函数,科学网—一种有效的最优化方法——Nelder-Mead单纯形直接搜索算法 - 王福昌的博文...
虽然MATLAB本身自带了fminsearch()函数,可以求解目标函数无梯度的最优化问题,但是感觉下面的程序在很多时候更好用,特别是自变量有边界和非线性约束的时候. 这里是John D'Errico ...
- python牛顿法解非线性方程组_科学网—求解多元非线性方程组F(x)=0的Newton-Raphson方法及其MATLAB实现 - 王福昌的博文...
科学网对公式支持不太好,在博客园有相同博文 牛顿迭代法可以推广到多元非线性方程组 $boldsymbol{F}(boldsymbol{x})=boldsymbol{0}$的情况,称为牛顿-- 拉夫逊方 ...
- matlab trapz二重积分函数_科学网—MATLAB中的数值积分方法 - 王福昌的博文
实际应用中在MATLAB里面都有开发好的命令可以使用,如 quad(), quadl(),quad2d(),triplequad() .需要掌握这些命令的用法. 1. 定积分 trapz(),qua ...
- matlab对数收益直方图,科学网—MATLAB中绘制数据直方图的新函数histogram2 - 王福昌的博文...
MATLAB中有命令hist3() 可以绘制直方图,竖坐标是频数,这与一些教科书中用纵轴表示频率的做法不一致,有些时候不便于使用.当然,使用者可以自己编写定制能够在纵轴绘出频率的直方图.在MATLAB ...
- matlab统计水文参数,科学网—[转载]利用MATLAB计算水文极值 - 刘朋的博文
利用MATLAB计算水文极值(年最大值,年连续5日最大,连续干/湿日,连续极端径流低值日数) [filename,filepath]=uigetfile('*.*','请选择文件'); %计算水文极端 ...
- hdc mfc 画扇形图_科学网—画扇形图(idl程序) - 张国印的博文
IDL画扇形图还是有些麻烦的,今天中午没午休,以红移和RA为例写了程序,希望以后能用上 pro sector set_plot,'ps' device,file='F:Aprilmap.ps' REA ...
- ipadpython代码_科学网—如何用iPad运行Python代码? - 王树义的博文
其实,不只是iPad,手机也可以. 痛点 我组织过几次线下编程工作坊,带着同学们用Python处理数据科学问题. 其中最让人头疼的,就是运行环境的安装. 实事求是地讲,参加工作坊之前,我已经做了认真准 ...
- matlab print 白边,科学网-[原创] matlab输出图片无白边-杨光的博文
今天要做一个gif动画,可惜GIF Movie Gear不认eps文件,无奈只好输出png格式的文件,麻烦来了,输出的图像有白边!之前挥之不去的问题再一次来了.在网上搜索一个多小时,都是说什么先ims ...
- matlab 画qq图,科学网—[转载]R语言绘制QQ图 - 刘朋的博文
R语言绘制QQ图 实例1: #############加载数据 data R R=apply(R,2,as.numeric) #R语言将字符串矩阵转化为数值型矩阵,apply()函数里面的第2个值,如 ...
最新文章
- HAProxy+Keepalived高可用负载均衡配置
- 山东省百万奖金赛事来了!
- php 的包管理,php composer包管理器
- PCL: 根据几何规则的曲面剖分-贪婪法表面重建三角网格
- mysql大数据量分页的一些做法
- Apache ab并发负载压力测试
- [导入]WCF后传系列(8):深度通道编程模型Part 1—设计篇
- webmin升级php,Centos linux下webmin安装及配置
- 关于给电鼓音源增加鼓盘或者DIY鼓盘(DIY镲片)的方法
- 上海大华条码称代码_上海大华条码秤的调试方法
- win10默认壁纸_Win10系统待机锁频壁纸怎么提取?
- 云计算是一种商业模式
- Kotlin入门-数据类与密封类 的解脱,由繁至简
- arm linux not syncing,Linux系统启动中途停止,提示Kernel panic - not syncing: Attempted to kill init!...
- Python经典实验4-字典和集合的应用
- border:none以及border:0的区别
- “75后”大学副校长,当选院士!
- 展锐芯片之GPU频率(一百一十四)
- 项目开发合作协议合同范本
- 《金瓶梅》是何等小说
热门文章
- Sniffer的讨论
- SAP ABAP 业务对象 BUS6038 AssetDownPayment 资产:预付款 BAPI 清单和相关 TCODE
- 什么是java类,类怎么理解,类的含义
- 一键ghost提示“Cannot open image file'1.4 I:/~1/c_pan.gho'”解决方法
- Hadoop国内镜像下载地址
- HTML+css3个人博客html源码
- 状态管理的概念,都是纸老虎
- 黑苹果N卡NVIDIA Web驱动合集,随时更新!
- C语言程序设计——实验教学大纲
- Linux离线安装Kafka(超级精简亲测安装)