使用 Optim 拟合二次函数
Using Optmi to fit a quadrtic function
我需要使用 optim()
将二次函数 $f(x) = ax^2+bx+c$ 拟合到数据,并指定梯度函数。我编写了以下有效的代码:
{r linewidth=80}
load(url("http://www.stat.cmu.edu/~pfreeman/HW_07_Q7.Rdata"),verbose=TRUE)
set.seed(101)
my.fit.fun = function(my.par)
{
sum(sqrt(abs(my.par[1]*x^2+my.par[2]*x+my.par[3]-y^2)))
}
my.par = c(0.2,-4,-5)
gradient=function(my.par){
c(my.par[1]*2,my.par[2],0)
}
optim.out = optim(my.par,fn=my.fit.fun, gr=gradient, method = "BFGS")
round(optim.out$par,3)
round(optim.out$value,3)
plot(x,y)
lines(x,optim.out$par[1]*x^2+optim.out$par[2]*x+optim.out$par[3])
Loading objects:
x
y
[1] 0.192 0.038 -5.000
[1] 38.785
然而,指令说梯度是拟合度量的梯度,而不是$f(x)$的导数,并且梯度函数returns长度为三的向量:偏导数$a$ 的拟合度量,然后是 $b$,然后是 $c$。我很困惑这里的梯度函数到底是什么?我不认为我的 gr=gradient 是正确的(即使我的图表看起来不错)所以请有人就这部分提出建议。谢谢。
首先,我宁愿对函数使用平方和而不是绝对值。您可以执行以下操作:
x <- 1:10
y < c(-0.2499211,-4.6645685,-2.6280750,-2.0146818,1.5632500,0.2043376,2.9151158, 4.0967775,6.8184074,12.5449975)
d <- data.frame(x,y)
fun <- function(par, data){
y_hat <- data$x^2 * par[1] + data$x * par[2] + par[3]
sum((data$y - y_hat)^2)
}
optim(c(0.2,-4,-5), fun, data = d)
$par
[1] 0.2531111 -1.3135297 -0.6618520
$value
[1] 17.70251
$counts
function gradient
176 NA
$convergence
[1] 0
$message
NULL
我不会使用 optim
,而是使用 nls
。在这里,您只需提供公式。在这种情况下,我们将有:
nls(y~ a * x^2 + b * x + c, d, c(a=0.2, b=-4, c=-5))
Nonlinear regression model
model: y ~ a * x^2 + b * x + c
data: d
a b c
0.2532 -1.3147 -0.6579
residual sum-of-squares: 17.7
Number of iterations to convergence: 1
Achieved convergence tolerance: 2.816e-08
另外,为什么要从 0.2,-4, -5
任何先验知识开始?例如,如果您在优化中使用 0,0,0,您将得到 nls
结果
编辑:
既然你想要 BFGS 方法,你可以这样做:
fun <- function(par, data){
y_hat <- data$x^2 * par[1] + data$x * par[2] + par[3]
sum((y_hat - data$y)^2)
}
grad <- function(par, data){
y_hat <- data$x^2 * par[1] + data$x * par[2] + par[3]
err <- data$y - y_hat
-2 * c(sum(err * data$x^2), sum(err * data$x), sum(err))
}
optim(c(0.2,-4,-5), fun, grad,data = d, method = "BFGS")
$par
[1] 0.2531732 -1.3146636 -0.6579553
$value
[1] 17.70249
$counts
function gradient
38 7
$convergence
[1] 0
$message
NULL
我需要使用 optim()
将二次函数 $f(x) = ax^2+bx+c$ 拟合到数据,并指定梯度函数。我编写了以下有效的代码:
{r linewidth=80}
load(url("http://www.stat.cmu.edu/~pfreeman/HW_07_Q7.Rdata"),verbose=TRUE)
set.seed(101)
my.fit.fun = function(my.par)
{
sum(sqrt(abs(my.par[1]*x^2+my.par[2]*x+my.par[3]-y^2)))
}
my.par = c(0.2,-4,-5)
gradient=function(my.par){
c(my.par[1]*2,my.par[2],0)
}
optim.out = optim(my.par,fn=my.fit.fun, gr=gradient, method = "BFGS")
round(optim.out$par,3)
round(optim.out$value,3)
plot(x,y)
lines(x,optim.out$par[1]*x^2+optim.out$par[2]*x+optim.out$par[3])
Loading objects:
x
y
[1] 0.192 0.038 -5.000
[1] 38.785
然而,指令说梯度是拟合度量的梯度,而不是$f(x)$的导数,并且梯度函数returns长度为三的向量:偏导数$a$ 的拟合度量,然后是 $b$,然后是 $c$。我很困惑这里的梯度函数到底是什么?我不认为我的 gr=gradient 是正确的(即使我的图表看起来不错)所以请有人就这部分提出建议。谢谢。
首先,我宁愿对函数使用平方和而不是绝对值。您可以执行以下操作:
x <- 1:10
y < c(-0.2499211,-4.6645685,-2.6280750,-2.0146818,1.5632500,0.2043376,2.9151158, 4.0967775,6.8184074,12.5449975)
d <- data.frame(x,y)
fun <- function(par, data){
y_hat <- data$x^2 * par[1] + data$x * par[2] + par[3]
sum((data$y - y_hat)^2)
}
optim(c(0.2,-4,-5), fun, data = d)
$par
[1] 0.2531111 -1.3135297 -0.6618520
$value
[1] 17.70251
$counts
function gradient
176 NA
$convergence
[1] 0
$message
NULL
我不会使用 optim
,而是使用 nls
。在这里,您只需提供公式。在这种情况下,我们将有:
nls(y~ a * x^2 + b * x + c, d, c(a=0.2, b=-4, c=-5))
Nonlinear regression model
model: y ~ a * x^2 + b * x + c
data: d
a b c
0.2532 -1.3147 -0.6579
residual sum-of-squares: 17.7
Number of iterations to convergence: 1
Achieved convergence tolerance: 2.816e-08
另外,为什么要从 0.2,-4, -5
任何先验知识开始?例如,如果您在优化中使用 0,0,0,您将得到 nls
结果
编辑:
既然你想要 BFGS 方法,你可以这样做:
fun <- function(par, data){
y_hat <- data$x^2 * par[1] + data$x * par[2] + par[3]
sum((y_hat - data$y)^2)
}
grad <- function(par, data){
y_hat <- data$x^2 * par[1] + data$x * par[2] + par[3]
err <- data$y - y_hat
-2 * c(sum(err * data$x^2), sum(err * data$x), sum(err))
}
optim(c(0.2,-4,-5), fun, grad,data = d, method = "BFGS")
$par
[1] 0.2531732 -1.3146636 -0.6579553
$value
[1] 17.70249
$counts
function gradient
38 7
$convergence
[1] 0
$message
NULL