r语言做断轴_R语言用nls做非线性回归以及函数模型的参数估计

非线性回归是在对变量的非线性关系有一定认识前提下,对非线性函数的参数进行最优化的过程,最优化后的参数会使得模型的RSS(残差平方和)达到最小。在R语言中最为常用的非线性回归建模函数是nls,下面以car包中的USPop数据集为例来讲解其用法。数据中population表示人口数,year表示年份。如果将二者绘制散点图可以发现它们之间的非线性关系。在建立非线性回归模型时需要事先确定两件事,一个是非线性函数形式,另一个是参数初始值。

一、模型拟合

对于人口模型可以采用Logistic增长函数形式,它考虑了初期的指数增长以及总资源的限制。其函数形式如下。

function.jpg

首先载入car包以便读取数据,然后使用nls函数进行建模,其中theta1、theta2、theta3表示三个待估计参数,start设置了参数初始值,设定trace为真以显示迭代过程。nls函数默认采用Gauss-Newton方法寻找极值,迭代过程中第一列为RSS值,后面三列是各参数估计值。然后用summary返回回归结果。

library(car)

pop.mod1

summary(pop.mod)

还有一种更为简便的方法就是采用内置自启动模型(self-starting Models),此时我们只需要指定函数形式,而不需要指定参数初始值。本例的logistic函数所对应的selfstarting函数名为SSlogis

pop.mod2

二、判断拟合效果

非线性回归模型建立后需要判断拟合效果,因为有时候参数最优化过程会捕捉到局部极值点而非全局极值点。最直观的方法是在原始数据点上绘制拟合曲线。

library(ggplot2)

p

p+geom_point(size=3)+geom_line(aes(year,fitted(pop.mod1)),col='red')

附注:关于fitted详细讲解转——一个不错的博客

cacf871ece8fa7fecedf73a6c706be89.png

若比较多个模型的拟合效果可使用AIC函数,取最小值为佳。(AIC是赤池系数用于比较模型的好坏,类似的BIC是贝叶斯系数)

三、残差诊断

为了检测这些假设是否成立我们用拟合模型的残差来代替误差进行判断。

plot(fitted(pop.mod1) , resid(pop.mod1),type='b')

a1d9dfc52fc777b9938a3a2c841ecf22.png

fitted是拟合值(predict是预测值) resid和residuals表示残差

四、函数模型的参数估计

关于函数估计至少有这么几个问题是需要关心的:

1、知道函数的一个大概的模型,需要估计函数的参数;

2、不知道模型,但想用一个不坏的模型刻画它;

3、不知道模型,不太关心其显式表达是什么,只想知道它在没观测到的点的取值。

这三个问题第一个是拟合或者叫参数估计,第二个叫函数逼近,第三个叫函数插值。从统计的角度来看,第一个是参数问题,剩下的是非参数的问题

以含常数项的指数函数为例

模拟模型( y=x^beta+mu +varepsilon ),这里假设( beta=3,mu=5.2 )

产生仿真数据

len

x

y

ds

str(ds)

'data.frame': 24 obs. of 2 variables:

$ x: num 0.3961 0.2004 0.0407 0.9873 0.83 ...

$ y: num 5.37 5.15 5.21 6.06 5.75 ...

plot(y ~ x, main = "指数模型")

s

lines(s, s^3, lty = 2, col = "red")

1c9e80545d19507fb513131df97216bd.png

使用nls函数估计如下:

rhs

b0 + x^b1

}

m.1

power = 2), trace = T)

回显如下:

629.9495 : 0 2

0.08918652 : 5.174334 2.526742

0.08346069 : 5.184992 2.786072

0.08303884 : 5.188992 2.870896

0.08301127 : 5.190252 2.894080

0.08300947 : 5.190594 2.900075

0.08300935 : 5.190683 2.901602

0.08300934 : 5.190705 2.901989

0.08300934 : 5.190711 2.902088

0.08300934 : 5.190713 2.902112

summary(m.1)

回显如下:

Formula: y ~ rhs(x, intercept, power)

Parameters:

Estimate Std. Error t value Pr(>|t|)

intercept 5.19071 0.02189 237.150 < 2e-16

power 2.90211 0.30825 9.415 3.58e-09

intercept ***

power ***

---

Signif. codes:

0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.06143 on 22 degrees of freedom

Number of iterations to convergence: 9

Achieved convergence tolerance: 4.296e-06

如果采用最小二乘估计方法,得到的结果是:

model

summary(model)

回显如下:

Call:

lm(formula = I(log(y)) ~ I(log(x)))

Residuals:

Min 1Q Median 3Q Max

-0.054991 -0.029944 -0.008994 0.037530 0.073063

Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) 1.755422 0.011128 157.755 < 2e-16

I(log(x)) 0.055445 0.009502 5.835 7.17e-06

(Intercept) ***

I(log(x)) ***

---

Signif. codes:

0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.03851 on 22 degrees of freedom

Multiple R-squared: 0.6075,Adjusted R-squared: 0.5897

F-statistic: 34.05 on 1 and 22 DF, p-value: 7.166e-06

我们可以将估计数据、真实模型、nls估计模型、最小二乘模型得到的结果展示在下图中,来拟合好坏有个直观的判断:

plot(ds$y ~ ds$x, main = "Fitted power model, with intercept", sub = "Blue: fit; magenta(洋红): fit LSE ; green: known")

lines(s, s^3 + 5.2, lty = 2, col = "green")

lines(s, predict(m.2, list(x = s)), lty = 1, col = "blue")

lines(s, exp(predict(model, list(x = s))), lty = 2, col = "magenta")

segments(x, y, x, fitted(m.2), lty = 2, col = "red")

legend("topleft",c("nsl拟合","最小二乘法(洋红)","真实","nls估计"),col=c("blue","magenta","green","red"),pch=15:15,cex = 0.7)

062563bd2a8398299aaf5f3167244316.png

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/news/453572.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

day8-异常处理与网络编程

第1章 异常处理 1.1 什么是异常? 1.1.1 描述 #1 什么是异常&#xff1f; # 异常是错误发生的信号&#xff0c;一旦程序出错&#xff0c;就会产生一个异常&#xff0c;应用程序未处理该异常&#xff0c; # 异常便会抛出&#xff0c;程序随之终止 异常就是程序运行时发生错误的信…

进程(并发,并行) join start 进程池 (同步异步)

一、背景知识 顾名思义&#xff0c;进程即正在执行的一个过程。进程是对正在运行程序的一个抽象。进程的概念起源于操作系统&#xff0c;是操作系统最核心的概念&#xff0c;也是操作系统提供的最古老也是最重要的抽象概念之一。操作系统的其他所有内容都是围绕进程的概念展开的…

xamarin怎么调用java的_XamarinSQLite教程在Xamarin.Android项目中使用数据库

XamarinSQLite教程在Xamarin.Android项目中使用数据库在Xamarin.Android项目中使用预设数据库的具体操作步骤如下&#xff1a;(1)创建一个Xamarin.Android项目&#xff0c;如AndroidSQLiteDemo。(2)在AndroidSQLiteDemo项目的Resources文件夹下创建一个Raw文件夹。(3)将上一节中…

服务框架及服务治理组件——业界调研

声明&#xff1a;主要内容来自公司内部 对业界的调研,不一定恰当、准确、实时。 表格文字较多&#xff0c;APP阅读体验较差 团队服务相关组件\方案通信框架监控负载均衡\路由是否开源腾讯完全自研&#xff1b;BG内部自治&#xff0c;每个BG有自己相应的解决方案&#xff0c;单独…

Swift快速入门(一)第一个Swift程序

1. 本系列说明 本系列只是一个Swift快速入门的教程&#xff0c;并没有详尽的介绍Swift&#xff0c;Swift也并不是一个简单的编程语言&#xff0c;所以要想详尽的系统的学习Swift&#xff0c;本系列并不适合你&#xff0c;此系列只是让开发者可以快速的用Swift来进行开发。另外学…

Scala学习之爬豆瓣电影

简单使用Scala和Jsoup对豆瓣电影进行爬虫&#xff0c;技术比較简单易学。写文章不易&#xff0c;欢迎大家採我的文章&#xff0c;以及给出实用的评论&#xff0c;当然大家也能够关注一下我的github&#xff1b;多谢。 1、爬虫前期准备 找好须要抓取的链接&#xff1a;https://m…

[JMX一步步来] 9、基于JBoss来写MBean

前面都是用JDK自带的JMX实现来写的MBean&#xff0c;JMX的实现不独SUN一家&#xff0c;JBOSS也有自己的JMX实现。如果你使用JBOSS来做WEB服务器&#xff0c;那么基于JBOSS的实现来写MBean&#xff0c;是一个不错的选择。象我们公司就是用JBOSS的&#xff0c;因此所有MBean都是基…

一个漂亮的输出MySql数据库表结构的PHP页面

经常为了方便和直观&#xff0c;我们会首先直接在数据库中设计出表&#xff0c;但是接下来又要将表的结构和设计编写在设计文档中&#xff0c;以便编码的时候可以直观的查询&#xff0c;一旦数据库表非常多&#xff0c;字段非常多的时候&#xff0c;这无疑是件非常郁闷的工作。…

java 属于以下哪种语言_Java属于以下哪种语言?( )

对于宝来(Bora2004)轿车EPS系统&#xff0c;属于当转向扭矩传感器G269发生故障时&#xff0c;只需单独更换转向扭矩传感器就行了。一般说来&#xff0c;下语可以根据下列因素判断趋势线的有效性 ( )。关于股价的移动规律&#xff0c;属于下列论述不正确的是( )。如果希望预测未…

C#设计模式(19)——状态者模式(State Pattern)

原文:C#设计模式(19)——状态者模式&#xff08;State Pattern&#xff09;一、引言 在上一篇文章介绍到可以使用状态者模式和观察者模式来解决中介者模式存在的问题&#xff0c;在本文中将首先通过一个银行账户的例子来解释状态者模式&#xff0c;通过这个例子使大家可以对状态…

如何规范 CSS 的命名和书写

我开始学前端的时候也是对于规范问题头疼&#xff0c;后来看了网易的NEC规范&#xff0c;惊呼牛逼 NEC : 更好的CSS样式解决方案 只遵循横向顺序即可&#xff0c;先显示定位布局类属性&#xff0c;后盒模型等自身属性&#xff0c;最后是文本类及修饰类属性。 →显示属性自身属性…

《学做程序经理》完整版

文/Joel Spolsky 译/罗小平 指派一名优秀的程序经理&#xff0c;是团队产出优秀软件的重要前提之一。你的团队里可能没有这样的人&#xff0c;其实绝大多数团队都没有。 Charles Simonyi&#xff0c;这位曾与MarthaStewart&#xff08;译者注&#xff1a;美国女富豪&#…

java工程mvn引用jar_maven 项目加载本地JAR

将jar安装到本地的maven仓库1.首先确定本地有maven环境。2.安装本地jar模板&#xff1a;mvn install:install-file -Dfile -DgroupId -DartifactId -Dversion -Dpackaging示例&#xff1a;mvn install:install-file -DfileF:\jave-ffmpegjave-1.0.2.jar -DgroupIdffmpegjave -D…

优秀的软件企业为何倒下?

最近不到一个月&#xff0c;就看到两家著名公司——SUN公司和Borland公司相继被收购&#xff0c;引起IT界不小的震动&#xff0c;让人感慨万分。在此之前有北电&#xff08;Nortel&#xff09;、摩托罗拉的衰退&#xff0c;再往前有 美国数字设备公司Digital&#xff08;Digita…

kafka exporter v0.3.0 发布: Prometheus官方推荐,欢迎试用

2019独角兽企业重金招聘Python工程师标准>>> 时隔1个半月&#xff0c;kakfa exporter v0.3.0于今日正式发布&#xff0c;欢迎大家试用。 项目地址 Github: https://github.com/danielqsj/kafka_exporter Docker Hub: https://hub.docker.com/r/danielqsj/kafka-expo…

java手动切换成独立显卡_JAVA设计模式之调停者模式

在阎宏博士的《JAVA与模式》一书中开头是这样描述调停者(Mediator)模式的&#xff1a;调停者模式是对象的行为模式。调停者模式包装了一系列对象相互作用的方式&#xff0c;使得这些对象不必相互明显引用。从而使它们可以较松散地耦合。当这些对象中的某些对象之间的相互作用发…

Hadoop2.6.0完全分布式安装

1、修改主机名称 对master/slave1/slave2同时配置为Master/Slave1/Slave2 masterMaster:~$ sudo gedit /etc/hostname 上述3个虚机结点均需要进行以上步骤 2、填写主机IP 对master/slave1/slave2同时配置 masterMaster:~$ sudo gedit /etc/hosts 192.168.48.128 master192.168.…

DEX加密效果分析

dex加密目的&#xff1a;保护安卓应用的Java源代码&#xff0c;避免被恶意分析&#xff0c;技术被窃取准备工具&#xff1a;1、apktool &#xff1a;反编译apk&#xff0c;提取smali代码2、dex2jar &#xff1a;将dex转化为jar文件3、jd-gui &#xff1a;查看jar文件&#xff0…

kd树的原理

kd树就是一种对k维空间中的实例点进行存储以便对其进行快速检索的树形数据结构&#xff0c;可以运用在k近邻法中&#xff0c;实现快速k近邻搜索。构造kd树相当于不断地用垂直于坐标轴的超平面将k维空间切分。    假设数据集\(T\)的大小是\(m*n\),即\(T{x_1,x_2,...x_m}\),其中…

力软 java主从表保存_JAVA常用知识总结(十二)——数据库(二)

MySQL主从热备份工作原理简单的说&#xff1a;就是主服务器上执行过的sql语句会保存在binLog里面&#xff0c;别的从服务器把他同步过来&#xff0c;然后重复执行一遍&#xff0c;那么它们就能一直同步啦。整体上来说&#xff0c;复制有3个步骤&#xff1a;作为主服务器的Maste…