optimization

optimazation is usually used by the \(optim()\) function to perform the optimization.
we could write the ordinary least square estimation and also the maximum likelihood and maximum a postires.
See the following examples.

1.Linear Regression OLS (e.g \(y=\theta_0+\theta_1 \times x_1 + \theta_2 \times x_2\))

by Default is to minimize the function.

L2.norm<-function(data,coef){
  with(data,sum((y-coef[1]-coef[2]*x1-coef[3]*x2)^2))
}

test=data.frame(x1=c(1,2,3,4,5,6,9), 
               x2=c(3,1,8,2,2,7,11), 
               y=c(1,3,5,6,8,12,6))

opt.model<-optim(par = c(0,1,1),L2.norm,data=test)
lm.model<-lm(formula = y~x1+x2,data = test)

opt.model$par
## [1]  2.4322797  1.0591250 -0.2294115
lm.model$coefficients
## (Intercept)          x1          x2 
##   2.4331797   1.0587558  -0.2292627
opt.model$convergence
## [1] 0
#0 indicates good convergence

nls()

xdata = c(-2,-1.64,-1.33,-0.7,0,0.45,1.2,1.64,2.32,2.9)
ydata = c(0.699369,0.700462,0.695354,1.03905,1.97389,2.41143,1.91091,0.919576,-0.730975,-1.42001)
plot(xdata,ydata)

\(y=p_1 \times cos(p_2 \times x) +p_2 \times sin(p_1 \times x)\)

#try some initial value
fit = nls(ydata ~ p1*cos(p2*xdata) + p2*sin(p1*xdata), start=list(p1=1,p2=0.2))
summary(fit)
## 
## Formula: ydata ~ p1 * cos(p2 * xdata) + p2 * sin(p1 * xdata)
## 
## Parameters:
##    Estimate Std. Error t value Pr(>|t|)    
## p1 1.881851   0.027430   68.61 2.27e-12 ***
## p2 0.700230   0.009153   76.51 9.50e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.08202 on 8 degrees of freedom
## 
## Number of iterations to convergence: 7 
## Achieved convergence tolerance: 2.189e-06
plot(xdata,ydata)
new = data.frame(xdata = seq(min(xdata),max(xdata),len=200))
lines(xdata,predict(fit,newdata=new$xdata))

len <- 24
x = runif(len)
y = x^3 + runif(len, min = -0.1, max = 0.1)
plot(x, y)
s <- seq(from = 0, to = 1, length = 50)
lines(s, s^3, lty = 2)

df <- data.frame(x, y)
m <- nls(y ~ I(x^power), data = df, start = list(power = 1), trace = T)
## 2.153021 :  1
## 0.3108513 :  2.122466
## 0.06309841 :  2.961643
## 0.05751992 :  3.144276
## 0.05751981 :  3.143411
## 0.05751981 :  3.143448
summary(m)
## 
## Formula: y ~ I(x^power)
## 
## Parameters:
##       Estimate Std. Error t value Pr(>|t|)    
## power   3.1434     0.1288   24.41   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.05001 on 23 degrees of freedom
## 
## Number of iterations to convergence: 5 
## Achieved convergence tolerance: 2.514e-06
lines(s, predict(m, list(x = s)), col = "green")

Details

Function I has two main uses.

In function data.frame. Protecting an object by enclosing it in I() in a call to data.frame inhibits the conversion of character vectors to factors and the dropping of names, and ensures that matrices are inserted as single columns. I can also be used to protect objects which are to be added to a data frame, or converted to a data frame via as.data.frame.

It achieves this by prepending the class “AsIs” to the object’s classes. Class “AsIs” has a few of its own methods, including for [, as.data.frame, print and format.

In function formula. There it is used to inhibit the interpretation of operators such as “+”, “-”, “*" and “^” as formula operators, so they are used as arithmetical operators. This is interpreted as a symbol by terms.formula.

(RSS.p <- sum(residuals(m)^2)) 
## [1] 0.05751981
(TSS <- sum((y - mean(y))^2)) 
## [1] 1.700365
1 - (RSS.p/TSS)
## [1] 0.9661721
m.exp<-nls(y ~ I(a * exp(b * x)), data = df, start = list(a = 1, b = 0), trace = T)
## 15.03929 :  1 0
## 6.246962 :  -0.1091652  0.7264928
## 3.408804 :  -0.09159764 -3.54043390
## 3.291543 :   -0.09570517 -37.46973193
## 3.247068 :    0.05818516 -69.42244658
## 3.239722 :   0.05891477 -8.70498552
## 2.670362 :  0.04957494 0.09054465
## 1.392467 :  0.01028939 3.56552552
## 0.1028067 :  0.00937403 4.74555979
## 0.07244489 :  0.01358137 4.37445195
## 0.05789463 :  0.01368032 4.44997640
## 0.05787189 :  0.01374887 4.44123798
## 0.05787188 :  0.01374357 4.44169841
## 0.05787188 :  0.01374386 4.44167334
plot(x, y)
lines(s, s^3, lty = 2)
lines(s, predict(m, list(x = s)), col = "green")
lines(s, predict(m.exp, list(x = s)), col = "red")

exp.eq <- function(x, a, b) {
    exp(1)^(a + b * sin(x^4))
}
exp.eq(2, 1, 3)  # testing the equation works: 3^2 + 1 = 10
## [1] 1.146014
m.sinexp <- nls(y ~ exp.eq(x, a, b), data = df, start = list(a = 1, b = 1), 
    trace = T)
## 240.302 :  1 1
## 28.77648 :  0.0260631 1.2064371
## 2.966551 :  -0.8970697  1.6580154
## 0.3549156 :  -1.664932  2.360828
## 0.2006881 :  -2.067160  2.860641
## 0.1993663 :  -2.110173  2.905472
## 0.1993615 :  -2.10701  2.89939
## 0.1993614 :  -2.107554  2.900366
## 0.1993614 :  -2.107467  2.900210
## 0.1993614 :  -2.107481  2.900235
plot(x, y)
lines(s, s^3, lty = 2)
lines(s, predict(m, list(x = s)), col = "green")
lines(s, predict(m.exp, list(x = s)), col = "red")
lines(s, predict(m.sinexp, list(x = s)), col = "blue")