注意：关闭R之前务必保存工作空间，保证学习的连续性。这样以前数据的控制台命令执行的效果以及相关变量仍然保存在内存中。
1 访问数据框变量
建议：在read.table命令执行names查看要处理的变量
& names(Squid) [1] "Sample" & "Year" & & "Month" & &"Location" "Sex" & & &"GSI" & && &&
1.1 str函数
& & str函数可以查看数据框中每个变量的属性：
& str(Squid) 'data.frame': & 2644 obs. of &6 variables: &$ Sample &: int &1 2 3 4 5 6 7 8 9 10 ... &$ Year & &: int &1 1 1 1 1 1 1 1 1 1 ... &$ Month & : int &1 1 1 1 1 1 1 1 1 2 ... &$ Location: int &1 3 1 1 1 1 1 3 3 1 ... &$ Sex & & : int &2 2 2 2 2 2 2 2 2 2 ... &$ GSI & & : num &10.44 9.83 9.74 9.31 8.99 ... &&
Sample &,Yead,Month,Location,Sex这几个变量是整型
GSI这个变量是数值型
GSI这个变量是存在于数据框Squid中的，不能通过在R控制台中输入GSI查看
& GSI 错误: 找不到对象'GSI' &&
1.2 函数中的数据参数--访问数据框中的变量的最佳方式
& M1 &- lm(GSI ~ factor(Location)+factor(Year),data = Squid) & M1
Call: lm(formula = GSI ~ factor(Location) + factor(Year), data = Squid)
Coefficients: & & & (Intercept) &factor(Location)2 &factor(Location)3 &factor(Location)4 & & & & & & &1.3939 & & & & & &-2.2178 & & & & & &-0.1417 & & & & & & 0.3138 & & & factor(Year)2 & & &factor(Year)3 & & &factor(Year)4 & & & & & & &1.3548 & & & & & & 0.9564 & & & & & & 1.2270 &
lm 是做线性回归的函数，data = Squid表示从数据框Squid中取变量
data = 并不是适用于任何函数，eg：
& mean(GSI,data = Squid)& 错误于mean(GSI, data = Squid) : 找不到对象'GSI' && 1.3 $ 符号 访问变量的另外一种方法
Squid$GSI&
& Squid$GSI & &[1] 10.1 &9.7 &8.7 &8.5 & &[9] &7.2 &6.2 &6.6 &5.0 & [17] &5.7 &5.1 &5.0 &5.0 & [25] &5.7 &5.2 &5.2 &5.5
& Squid[,6] & &[1] 10.1 &9.7 &8.7 &8.5 & &[9] &7.2 &6.2 &6.6 &5.0 & [17] &5.7 &5.1 &5.0 &5.0 & [25] &5.7 &5.2 &5.2 &5.5
此时可以通过mean求平均值
& mean(Squid$GSI)
[1] 2.187034
1.4 attach 函数
attach函数将数据框添加到R的搜索路径中，此时就可以通过GSI命令直接查看GSI数据
& attach(Squid) & GSI & &[1] 10.1 &9.7 &8.7 &8.5 & &[9] &7.2 &6.2 &6.6 &5.0 & [17] &5.7 &5.1 &5.0 &5.0 & [25] &5.7 &5.2 &5.2 &5.5
此时就可以直接使用相关函数了。
& boxplot(GSI) &&
(额、、看不懂这个图)
使用attach函数显然应该小心保证变量名字的唯一性，如果与R自带函数名字或者变量一样肯定会出问题。
attach使用总结：
(1)为了避免复制变量，避免输入Squid$GSI两次以上
(2)使用attach命令应该保证变量的唯一性
(3)如果要处理多个数据集,而且一次只处理一个数据集，使用detach函数将数据集从R搜索路径中删除
2 访问数据集
首先执行detach(Squid)命令！！！
查看Squid中Sex的值
& Squid$Sex & &[1] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2 2 2 2 2 2 2 2 2 & [36] 2 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2 2 2 2 2 2 2 1 2 1 1 1 1 2 1 1 1 1 1 1
显示位移值
& unique(Squid$Sex) [1] 2 1 &&
其中1表示雄性2表示雌性
& Sel &- Squid$Sex == 1 & SquidM &- Squid[Sel,] & SquidM & & &Sample Year Month Location Sex & &GSI 24 & & & 24 & &1 & & 5 & & & &1 & 1 5.2970 48 & & & 48 & &1 & & 5 & & & &3 & 1 4.2968 58 & & & 58 & &1 & & 6 & & & &1 & 1 3.5008 60 & & & 60 & &1 & & 6 & & & &1 & 1 3.2487 61 & & & 61 & &1 & & 6 & & & &1 & 1 3.2304
Sel &- Squid$Sex == 1这条命令生成一个向量与Sex具有相同的长度,如果Sex的值等于1则该变量的值为TRUE，否则为FALSE，这样一个变量可称为布尔变量，可以用来选择行。
SquidM &- Squid[Sel,]这条命令表示选择Squid中Sel等于TRUE的行，并将数据存储到SquidM中。因为是选择行，所以需要使用方阔号。
第三章未完待续...
获得雌性数据
& SquidF &- Squid[Squid$Sex == 2,] & SquidF & & &Sample Year Month Location Sex & & GSI 1 & & & & 1 & &1 & & 1 & & & &1 & 2 10.4432 2 & & & & 2 & &1 & & 1 & & & &3 & 2 &9.8331 3 & & & & 3 & &1 & & 1 & & & &1 & 2 &9.7356 4 & & & & 4 & &1 & & 1 & & & &1 & 2 &9.3107 5 & & & & 5 & &1 & & 1 & & & &1 & 2 &8.9926
下面几条命令不解释：
unique(Squid$Location) Squid123 &- Squid[Squid$Location == 1 | Squid$Location ==2 | Squid$Location == 3,] Squid123 &- Squid[Squid$Location != 4,] Squid123 &- Squid[Squid$Location & 4 ,] Squid123 &- Squid[Squid$Location &=3 ,] Squid123 &- Squid[Squid$Location &=1 &Squid$Location &=3 ,]
都是获得Location值为1,2,3的行
& unique(Squid$Location) [1] 1 3 4 2 & Squid123 &- Squid[Squid$Location == 1 | Squid$Location ==2 | Squid$Location == 3,] & Squid123 & & &Sample Year Month Location Sex & & GSI 1 & & & & 1 & &1 & & 1 & & & &1 & 2 10.4432 2 & & & & 2 & &1 & & 1 & & & &3 & 2 &9.8331 3 & & & & 3 & &1 & & 1 & & & &1 & 2 &9.7356 4 & & & & 4 & &1 & & 1 & & & &1 & 2 &9.3107 5 & & & & 5 & &1 & & 1 & & & &1 & 2 &8.9926 6 & & & & 6 & &1 & & 1 & & & &1 & 2 &8.7707
获得Location值为1的雄性数据行
& SquidM.1 &- Squid[Squid$Sex == 1 & Squid$Location == 1,] & SquidM.1 & & &Sample Year Month Location Sex & &GSI 24 & & & 24 & &1 & & 5 & & & &1 & 1 5.2970 58 & & & 58 & &1 & & 6 & & & &1 & 1 3.5008 60 & & & 60 & &1 & & 6 & & & &1 & 1 3.2487
获得位置为1或2的雄性数据
& SquidM.12 &- Squid[Squid$Sex == 1 &( Squid$Location == 1 | Squid$Location == 2),] & SquidM.12 & & &Sample Year Month Location Sex & &GSI 24 & & & 24 & &1 & & 5 & & & &1 & 1 5.2970 58 & & & 58 & &1 & & 6 & & & &1 & 1 3.5008 60 & & & 60 & &1 & & 6 & & & &1 & 1 3.2487
& SquidM1 &- SquidM[Squid$Location == 1,]
& & & & Sample Year Month Location Sex & &GSI
24 & & & & &24 & &1 & & 5 & & & &1 & 1 5.2970
58 & & & & &58 & &1 & & 6 & & & &1 & 1 3.5008
..........
..........
NA & & & & &NA & NA & &NA & & & NA &NA & & NA NA.1 & & & &NA & NA & &NA & & & NA &NA & & NA NA.2 & & & &NA & NA & &NA & & & NA &NA & & NA NA.3 & & & &NA & NA & &NA & & & NA &NA & & NA NA.4 & & & &NA & NA & &NA & & & NA &NA & & NA ..........
原因分析：
之前得到的SquidM表示雄性数据，显然SquidM的行数与Squid$Location == 1 布尔向量的长度不一致。因此导出出现上面的现象。
2.1 数据排序
& Ord1 &- order(Squid$Month) & Squid[Ord1,] & & &Sample Year Month Location Sex & & GSI 1 & & & & 1 & &1 & & 1 & & & &1 & 2 10.4432 2 & & & & 2 & &1 & & 1 & & & &3 & 2 &9.8331 3 & & & & 3 & &1 & & 1 & & & &1 & 2 &9.7356 4 & & & & 4 & &1 & & 1 & & & &1 & 2 &9.3107
根据月份排序
&也可以只对一个变量进行排序
& Squid$GSI[Ord1] & &[1] 10.1 &9.7 &8.7 &8.5 & &[9] &7.2 &6.7 &1.7 &0.6 & [17] &0.8 &0.8 &0.9 &0.6 & [25] &0.2 &0.8 &0.1 &0.9
3 使用相同的标识符组合两个数据集
& setwd("E:/R/R-beginer-guide/data/RBook") & Sql1 &- read.table(file = "squid1.txt",header = TRUE) & Sql2 &- read.table(file = "squid2.txt",header = TRUE) & SquidMerged &- merge(Sql1,Sql2,by = "Sample") & SquidMerged & & &Sample & & GSI YEAR MONTH Location Sex 1 & & & & 1 10.4432 & &1 & & 1 & & & &1 & 2 2 & & & & 2 &9.8331 & &1 & & 1 & & & &3 & 2 3 & & & & 3 &9.7356 & &1 & & 1 & & & &1 & 2 4 & & & & 5 &8.9926 & &1 & & 1 & & & &1 & 2 5 & & & & 6 &8.7707 & &1 & & 1 & & & &1 & 2 6 & & & & 7 &8.2576 & &1 & & 1 & & & &1 & 2
merge 命令采用两个数据框Sql1 ,Sql2作为参数并使用变量Sample作为形同的标识符合并两个数据。merger函数还有一个选项是all,缺省状态值是FALSE:即如果Sql1或Sql2中的值有缺失,则将被忽略。如果all的值设置为TRUE,可能会产生NA值
Sql11 &- read.table(file = "squid1.txt",header = TRUE) Sql21 &- read.table(file = "squid2.txt",header = TRUE) SquidMerged1 &- merge(Sql11,Sql21,by = "Sample") SquidMerged1
额、、这里好像没有出现NA，看来是数据没有丢失
4 输出数据
通过write.table将数据输出为ascii文件
write.table(SquidM,file = "MaleSquid_wujiahua.txt",sep = " ",quote = &FALSE,append = FALSE,na = "NA")
查看工作目录，生成了一个MaleSquid_wujiahua.txt文件，
Sample Year Month Location Sex GSI 24 24 1 5 1 1 5.297 48 48 1 5 3 1 4.2968 58 58 1 6 1 1 3.5008 60 60 1 6 1 1 3.2487 61 61 1 6 1 1 3.2304
write.table第一个参数表示要输出的数据，第二参数是数据保存的文件名，sep = " " 宝成数据通过空格隔开，qoute=FALSE消除字符串的引号标识，na="NA"表示缺失值通过NA替换。append=TRUE表示把数据添加到文件的尾部
5 重新编码分类变量
& str(Squid) 'data.frame': & 2644 obs. of &6 variables: &$ Sample &: int &1 2 3 4 5 6 7 8 9 10 ... &$ Year & &: int &1 1 1 1 1 1 1 1 1 1 ... &$ Month & : int &1 1 1 1 1 1 1 1 1 2 ... &$ Location: int &1 3 1 1 1 1 1 3 3 1 ... &$ Sex & & : int &2 2 2 2 2 2 2 2 2 2 ... &$ GSI & & : num &10.44 9.83 9.74 9.31 8.99 ...
其中Sex和locaton的值确定，属于分类变量。
在数据框中一般根据分类变量生成新的变量
& Squid$fLocation &- factor(Squid$Location) & Squid$fSex &- factor(Squid$Sex) & Squid$fLocation & &[1] 1 3 1 1 1 1 1 3 3 1 1 1 1 1 1 1 3 1 3 1 3 1 1 1 1 1 1 1 1 1 1 1 1 1 1 & [36] 1 1 1 1 3 1 1 1 1 3 1 1 3 1 1 1 1 1 1 1 3 1 1 1 1 1 3 1 1 1 1 1 1 1 1
& Squid$fSex
& &[1] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2 2 2 2 2 2 2 2 2
& [36] 2 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2 2 2 2 2 2 2 1 2 1 1 1 1 2 1 1 1 1 1 1
& [71] 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 2 1 1 1
[ 1 1 1 1 1 1 1 2 1 1 2 1 1 2 1 1 1 2 1 2 1 1 2 1 1 2 1 2 2 1 1 1 1 [ 1 1 1 2 1 1 1 2 1 2 1 2 1 2 1 1 1 Levels: 1 2
&fLocation和fSex只是名义变量，f表示他们是因子
levels：1，2可以对其修改
& Squid$fSex &- factor(Squid$Sex,levels = c(1,2),labels = c("M","F")) & Squid$fSex & &[1] F F F F F F F F F F F F F F F F F F F F F F F M F F F F F F F F F F F & [36] F F F F F F F F F F F F M F F F F F F F F F M F M M M M F M M M M M M
[2556] F M M M M F F M M M M M M M F M M M M M M F M M F M M M F M M F M M M
[2591] M M M M M M M M M F M M F M M F M M M F M F M M F M M F M F F M M M M
[2626] M M M M M F M M M F M F M F M F M M M
Levels: M F
这样每个1被M替换，2被F替换
使用重新分类的因子变量
boxplot(GSI ~ fSex,data = Squid)
& M1 &- lm(GSI ~ fSex+fLocation,data = Squid) & M1
Call: lm(formula = GSI ~ fSex + fLocation, data = Squid)
Coefficients: (Intercept) & & & &fSexF & fLocation2 & fLocation3 & fLocation4 & & & &1.3593 & & & 2.0248 & & &-1.8552 & & &-0.1425 & & & 0.5876 &
& summary(M1)
lm(formula = GSI ~ fSex + fLocation, data = Squid)
Residuals:
-3.5 -0.9 11.2159
Coefficients:
Estimate Std. Error t value Pr(&|t|)
(Intercept)
&2e-16 ***
&2e-16 ***
fLocation2
&2e-16 ***
fLocation3
fLocation4
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2.415 on 2639 degrees of freedom
Multiple R-squared: 0.1759,
Adjusted R-squared: 0.1746
F-statistic: 140.8 on 4 and 2639 DF,
p-value: & 2.2e-16
& (才发现有这么一个插入脚本功能)
& M2 &- lm(GSI ~ factor(Sex)+factor(Location),data = Squid)
& summary(M2)
lm(formula = GSI ~ factor(Sex) + factor(Location), data = Squid)
Residuals:
-3.5 -0.9 11.2159
Coefficients:
Estimate Std. Error t value Pr(&|t|)
(Intercept)
&2e-16 ***
factor(Sex)2
&2e-16 ***
factor(Location)2 -1.85525
&2e-16 ***
factor(Location)3 -0.14248
factor(Location)4
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2.415 on 2639 degrees of freedom
Multiple R-squared: 0.1759,
Adjusted R-squared: 0.1746
F-statistic: 140.8 on 4 and 2639 DF,
p-value: & 2.2e-16
& 估计的参数是一致的，但是第二种方式占用的屏幕空间更大，传说在二阶，三阶交互作用时将是一个严重的问题。
& Squid$fLocation
[1] 1 3 1 1 1 1 1 3 3 1 1 1 1 1 1 1 3 1 3 1 3 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[36] 1 1 1 1 3 1 1 1 1 3 1 1 3 1 1 1 1 1 1 1 3 1 1 1 1 1 3 1 1 1 1 1 1 1 1
[ 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
Levels: 1 2 3 4 Levels:的顺序可以更改
& Squid$fLocation &- factor(Squid$Location,levels= c(2,3,1,4))
& Squid$fLocation
[1] 1 3 1 1 1 1 1 3 3 1 1 1 1 1 1 1 3 1 3 1 3 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[36] 1 1 1 1 3 1 1 1 1 3 1 1 3 1 1 1 1 1 1 1 3 1 1 1 1 1 3 1 1 1 1 1 1 1 1
[71] 1 1 1 1 1 3 1 1 3 1 1 3 1 1 3 1 1 1 1 1 1 1 1 1 1 1 1 3 1 3 3 3 1 3 1
] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
Levels: 2 3 1 4
boxplot(GSI ~ fLocation,data = Squid)
SquidM &- Squid[Squid$Sex == 1,]
SquidM &- Squid[Squid$fSex == "1"] 在定义了fSex这个因子之后上面两种写法都是一样的效果。
但是1有双引号是必须的，因为fSex是因子
定义新的变量之后也可以通过str命令查看
& Squid$fSex &- factor(Squid$Sex,labels = c("M","F"))
& Squid$fLocation &- factor(Squid$Location)
& str(Squid)
'data.frame':
2644 obs. of
8 variables:
1 2 3 4 5 6 7 8 9 10 ...
1 1 1 1 1 1 1 1 1 1 ...
1 1 1 1 1 1 1 1 1 2 ...
$ Location : int
1 3 1 1 1 1 1 3 3 1 ...
2 2 2 2 2 2 2 2 2 2 ...
10.44 9.83 9.74 9.31 8.99 ...
$ fLocation: Factor w/ 4 levels "1","2","3","4": 1 3 1 1 1 1 1 3 3 1 ...
: Factor w/ 2 levels "M","F": 2 2 2 2 2 2 2 2 2 2 ...
第三章总结：
write.table & &把一个变量写入到ascii文件中 & & &write.table(Squid,file="test.txt")
order & & & & & 确定数据的排序 & & & & & & & & & & &order(x)
merge & & & & &合并两个数据框 & & & & & & & & & & &merege(a,b,by="ID")&
str & & & & & & & 显示一个对象的内部结构 & & & & &str(Squid)
factor & & & & & 定义变量作为因子 & & & & & & & & &factor(x)
& 开源中国(OSChina.NET) |
开源中国社区(OSChina.net)是工信部
指定的官方社区后使用快捷导航没有帐号？
R语言与函数估计学习笔记（函数模型的参数估计）
查看: 25194|
评论: 0|来自: CSDN博客
摘要: 毫无疑问，函数估计是一个比参数估计要复杂得多的问题，当然也是一个有趣的多的问题。这个问题在模型未知的实验设计的建模中十分的常见，也是我正在学习的内容的一部分。关于函数估计我想至少有这么几个问题是我们关 ...
毫无疑问，函数估计是一个比参数估计要复杂得多的问题，当然也是一个有趣的多的问题。这个问题在模型未知的实验设计的建模中十分的常见，也是我正在学习的内容的一部分。关于函数估计我想至少有这么几个问题是我们关心的：1、我知道函数的一个大概的模型，需要估计函数的参数；2、我不知道它是一个什么模型，但是我想用一个不坏的模型刻画它；3、我不知道它是一个什么模型，我也不太关心它的显式表达是什么，我只想知道它在没观测到的点的取值。这三个问题第一个是拟合或者叫参数估计，第二个叫函数逼近，第三个叫函数插值。从统计的角度来看，第一个是参数问题，剩下的是非参数的问题。函数模型的参数估计这类的问题有很多，一个比较典型的例子是柯布-道格拉斯函数( Y=L^alpha k^beta mu )。我们要估计参数常用的就是最小化残差平方和，如果是密度函数或者分布函数常用的办法在加上矩估计与似然估计（MLE）两种办法。我们在这里介绍一下R中的用于函数拟合的函数nls(),其调用格式如下：nls(formula, data, start, control, algorithm, trace, subset, weights, na.action, model, lower, upper, …)其用法与线性回归函数lm()用法类似，这里就不作过多介绍了，我们来看几个例子来说明函数的用法：情形一：指数模型模拟模型( y=x^beta+varepsilon ),这里假设( beta=3 )len &- 24
x &- runif(len, 0.1, 1)
y &- x^3 + rnorm(len, 0, 0.06)
ds &- data.frame(x = x, y = y)
## 'data.frame':
24 obs. of
2 variables:
0.238 0.482 0.787 0.145 0.232 ...
0.48 0.87 -0.00321 ...
plot(y ~ x, main = "Known cubic, with noise")
s &- seq(0, 1, length = 100)
lines(s, s^3, lty = 2, col = "green")
使用函数nls估计参数( beta )m &- nls(y ~ I(x^power), data = ds, start = list(power = 1), trace = T)
## 1.637 :
## 0.2674 :
## 0.07229 :
## 0.06273 :
## 0.06264 :
## 0.06264 :
## 0.06264 :
summary(m)
## Formula: y ~ I(x^power)
## Parameters:
Estimate Std. Error t value Pr(&|t|)
&2e-16 ***
## Signif. codes:
0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Residual standard error: 0.0522 on 23 degrees of freedom
## Number of iterations to convergence: 6
## Achieved convergence tolerance: 6.07e-06
当然，也可以两边取对数，通过最小二乘来处理这个问题。其R代码如下：model &- lm(I(log(y)) ~ I(log(x)))
summary(model)
## lm(formula = I(log(y)) ~ I(log(x)))
## Residuals:
## -1.7 -0.6
## Coefficients:
Estimate Std. Error t value Pr(&|t|)
## (Intercept)
## I(log(x))
1.3e-06 ***
## Signif. codes:
0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Residual standard error: 0.695 on 21 degrees of freedom
(1 observation deleted due to missingness)
## Multiple R-squared:
Adjusted R-squared:
## F-statistic: 44.8 on 1 and 21 DF,
p-value: 1.27e-06
如果这个模型还有常数项，两边取对数就不好使了，不过，我们的nls函数还是能解决的。情形二：含常数项的指数模型模拟模型( y=x^beta+mu +varepsilon ),这里假设( beta=3,mu=5.2 )len &- 24
x &- runif(len)
y &- x^3 + 5.2 + rnorm(len, 0, 0.06)
ds &- data.frame(x = x, y = y)
## 'data.frame':
24 obs. of
2 variables:
0.277 0.831 0.127 0.464 0.734 ...
5.17 5.79 5.22 5.37 5.64 ...
plot(y ~ x, main = "Known cubic, with noise")
s &- seq(0, 1, length = 100)
lines(s, s^3, lty = 2, col = "green")
使用nls函数估计如下：rhs &- function(x, b0, b1) {
m.2 &- nls(y ~ rhs(x, intercept, power), data = ds, start = list(intercept = 0,
power = 2), trace = T)
## 632.5 :
## 0.05006 :
5.171 2.331
## 0.04934 :
5.173 2.395
## 0.04934 :
5.174 2.404
## 0.04934 :
5.174 2.404
## 0.04934 :
5.174 2.405
## 0.04934 :
5.174 2.405
summary(m.2)
## Formula: y ~ rhs(x, intercept, power)
## Parameters:
Estimate Std. Error t value Pr(&|t|)
## intercept
& 2e-16 ***
3.7e-12 ***
## Signif. codes:
0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Residual standard error: 0.0474 on 22 degrees of freedom
## Number of iterations to convergence: 6
## Achieved convergence tolerance: 1.67e-06
如果这时我们还是采用最小二乘估计的办法处理，那么得到的结果是：model &- lm(I(log(y)) ~ I(log(x)))
summary(model)
## lm(formula = I(log(y)) ~ I(log(x)))
## Residuals:
## -0.083 -0.040
## Coefficients:
Estimate Std. Error t value Pr(&|t|)
## (Intercept)
& 2e-16 ***
## I(log(x))
6.3e-06 ***
## Signif. codes:
0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Residual standard error: 0.0287 on 22 degrees of freedom
## Multiple R-squared:
Adjusted R-squared:
## F-statistic: 34.7 on 1 and 22 DF,
p-value: 6.32e-06
我们可以将估计数据、真实模型、nls估计模型、最小二乘模型得到的结果展示在下图中，来拟合好坏有个直观的判断：plot(ds$y ~ ds$x, main = "Fitted power model, with intercept", sub = "Blue: 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")
从图就可以看出，化为最小二乘的办法不总是可行的。情形三：分段函数模型我们来看下面的模型：f.lrp &- function(x, a, b, t.x) {
ifelse(x & t.x, a + b * t.x, a + b * x)
f.lvls &- seq(0, 120, by = 10)
b.0 &- 0.05
t.x.0 &- 70
test &- data.frame(x = f.lvls, y = f.lrp(f.lvls, a.0, b.0, t.x.0))
test &- rbind(test, test, test)
test$y &- test$y + rnorm(length(test$y), 0, 0.2)
plot(test$y ~ test$x, main = "Linear response and plateau yield response", xlab = "Fertilizer added",
ylab = "Crop yield")
(max.yield &- a.0 + b.0 * t.x.0)
## [1] 5.5
lines(x = c(0, t.x.0, 120), y = c(a.0, max.yield, max.yield), lty = 2)
abline(v = t.x.0, lty = 3)
abline(h = max.yield, lty = 3)
显然用一个线性模型解决不了，二次模型解决不好，分段函数倒是一个很好的选择，那么在哪里划分比较合理呢？我们还是用nls函数来解决这个问题：m.lrp &- nls(y ~ f.lrp(x, a, b, t.x), data = test, start = list(a = 0, b = 0.1,
t.x = 50), trace = T, control = list(warnOnly = T, minFactor = 1/2048))
## 32.74 :
## 7.352 :
2.119 59.34899
2.119 70.24081
## 1.116 :
2.139 72.09071
## 1.116 :
2.139 72.08250
summary(m.lrp)
## Formula: y ~ f.lrp(x, a, b, t.x)
## Parameters:
Estimate Std. Error t value Pr(&|t|)
&2e-16 ***
&2e-16 ***
## t.x 72.08250
&2e-16 ***
## Signif. codes:
0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Residual standard error: 0.176 on 36 degrees of freedom
## Number of iterations to convergence: 4
## Achieved convergence tolerance: 3.63e-09
画图来看看拟合的可靠性：plot(test$y ~ test$x, main = "Linear response and plateau yield response", xlab = "Fertilizer added",
ylab = "Crop yield")
(max.yield &- a.0 + b.0 * t.x.0)
## [1] 5.5
lines(x = c(0, t.x.0, 120), y = c(a.0, max.yield, max.yield), lty = 2, col = "blue")
abline(v = t.x.0, lty = 3, col = "blue")
abline(h = max.yield, lty = 3, col = "blue")
(max.yield &- coefficients(m.lrp)["a"] + coefficients(m.lrp)["b"] * coefficients(m.lrp)["t.x"])
lines(x = c(0, coefficients(m.lrp)["t.x"], 120), y = c(coefficients(m.lrp)["a"],
max.yield, max.yield), lty = 1)
abline(v = coefficients(m.lrp)["t.x"], lty = 4)
abline(h = max.yield, lty = 4)
text(120, 4, "known true model", col = "blue", pos = 2)
text(120, 3.5, "fitted model", col = "black", pos = 2)
可以看到拟合的结果还是不错的。这也显示了nls函数的优秀之处，几乎可以拟合所有的连续函数，哪怕他们存在不可微的点。它的算法是怎么样的我没有深究，不过光是分段线性模型，CART算法可是一个不错的选择，模型树（model tree）就是拟合这种模型的极好的选择。最近在整理机器学习的笔记，model tree的R代码确实是写好了，不过由于人懒，敲字慢，最终也没形成文字发出来与大家分享。我们对参数估计大概就介绍这么多，关于矩估计，极大似然估计可以参见之前的博文.当然，如果一个函数分离掉已知部分是一个密度函数的话，矩估计与极大似然仍然是可用的，如你想估计函数( f(x)=e^{-(x-mu)^2} )中的参数( mu )。
刚表态过的朋友 ()【图文】R语言回归分析_百度文库
两大类热门资源免费畅读
续费一年阅读会员，立省24元！
评价文档：
R语言回归分析
上传于||文档简介
&&一​元​线​性​回​归​和​多​元​线​性​回​归​,​还​包​括​残​差​分​析​、​异​常​值​处​理​等
大小：3.63MB
登录百度文库，专享文档复制特权，财富值每天免费拿！
你可能喜欢}