如何解决是否有R函数可以从“ beta + beta * t”之类的字符串创建公式?
我正在运行一种交叉验证算法,以找到适用于每天变化的数据的最佳多项式。我想找到一种不麻烦的方法来在简单的图中显示拟合,而不必每次都为该图手动编写整个回归公式和beta系数。对于回归公式,求解很容易,我使用sprintf创建一个字符串,并在字符串上使用as.formula()。
问题在于绘制线条。我以相同的方式创建了一个字符串,但是as.formula()函数似乎仅适用于回归公式,不适用于“ beta + beta * t”形式的公式。我还尝试过使用eval(parse()),如下所示,但这只会创建NA的向量。
#Create strings
poly_form = "y ~ t"
beta_form = "beta[1]"
for (i in 1:pmin) { #pmin is the best polynomial fit,e.g. 4 or 9.
poly_form <- sprintf("%s + I(t^%s)",poly_form,i)
beta_form <- sprintf("%s + beta[%s]*t^%s",beta_form,i+1,i)
}
#Regression
poly.mod = lm(as.formula(poly_form))
beta = coef(poly.mod)
#Plot
plot(t,y,type = 'h')
lines(t,eval(parse(text = beta_form))) #This doesn't work.
因此,从本质上讲,我如何以自动生成与以下内容相同的输出的方式,将作为输入的一部分创建的字符串使用到lines函数中:
lines(t,beta2[1] + beta2[2]*t + beta2[3]*t^2 + beta2[4]*t^3 + beta2[5]*t^4 + beta2[6]*t^5 + beta2[7]*t^6)
解决方法
这不是您的操作方式。
首先,使用poly
函数。其次,使用predict
。
set.seed(42)
y <- rnorm(10)
t <- 1:10
DF <- data.frame(y,t) #important!
pmin <- 3
poly.mod <- lm(y ~ poly(t,degree = pmin,raw = TRUE),data = DF)
plot(t,y,type = 'h')
curve(predict(poly.mod,newdata = data.frame(t = x)),add = TRUE)
curve
计算传递给其第一个参数的表达式。 x
表示图的x值。它总是必须为x
!
我认为罗兰(Roland)的方法在这里更好,但是获得解释为什么您自己的代码不起作用总是很高兴的。
让我们使用一些虚拟数据来具体说明一下,以便我们了解问题出在哪里:
set.seed(69)
t <- 1:100
y <- 3 + 0.3 * t + 0.01*t^2 + 0.0002*t^3 + 4e-6*t^4 +
3e-10*t^5 + 4e-16*t^6 + rnorm(100,50)
plot(t,y)
现在让我们想象一下,我们已经决定适合六阶多项式回归:
pmin <- 6
poly_form = "y ~ t"
beta_form = "beta[1]"
for (i in 1:pmin) { #pmin is the best polynomial fit,e.g. 4 or 9.
poly_form <- sprintf("%s + I(t^%s)",poly_form,i)
beta_form <- sprintf("%s + beta[%s]*t^%s",beta_form,i+1,i)
}
到目前为止,太好了。现在,让我们看看我们的多边形形式和Beta形式:
poly_form
#> [1] "y ~ t + I(t^1) + I(t^2) + I(t^3) + I(t^4) + I(t^5) + I(t^6)"
beta_form
# > [1] "beta[1] + beta[2]*t^1 + beta[3]*t^2 + beta[4]*t^3 + beta[5]*t^4 +
beta[6]*t^5 + beta[7]*t^6"
这里有些问题。我们在回归分析中包括了t
和{strong>的术语t^1
。这些当然是同一回事。因此,如果我们创建poly_mod
,我们将得到:
poly.mod = lm(as.formula(poly_form))
poly.mod
#> Call:
#> lm(formula = as.formula(poly_form))
#>
#> Coefficients:
#> (Intercept) t I(t^1) I(t^2) I(t^3) I(t^4)
#> -1.910e+00 -2.444e-01 NA -4.095e-02 5.933e-03 -1.499e-04
#> I(t^5) I(t^6)
#> 1.611e-06 -5.903e-09
您可以看到我们为NA
获得了I(t^1)
。但是,这意味着coef(poly.mod)
现在将包含一个NA
:
beta = coef(poly.mod)
beta
#> (Intercept) t I(t^1) I(t^2) I(t^3) I(t^4)
#> 8.139958e+01 -1.494928e+01 NA 1.037905e+00 -3.454374e-02 6.267641e-04
#> I(t^5) I(t^6)
#> -5.534399e-06 1.904566e-08
这意味着当我们解析beta_form
时,总和中总是有一个NA
,因此它只会产生一个NA
的向量:
eval(parse(text = beta_form))
#> [1] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
#> [28] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
#> [55] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
#> [82] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
那有什么解决方案?
只需将原来的poly_form = "y ~ t"
更改为poly_form = "y ~ "
。
现在,您按原样运行其余代码,并获得所需的结果:
plot(t,type = 'h')
lines(t,eval(parse(text = beta_form)))
,
使用poly()
:
model = lm(y ~ poly(t,4,raw = TRUE,data = df)
beta = coef(model)
t = t0 ^ (0:4)
sum(beta * t)
# or
predict(model,newdata) # dataframe of t
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。