n = 100
eps = rnorm(n, 0, 2)
alpha = 1
beta1 = 2
beta2 = 5
beta3 = 0.1
beta4 = 0.01
x = runif(n, 0, 2)
y = alpha * x^beta1 * exp(beta2 * x + beta3 * x^2 + beta3 * x^3) * exp(eps)
y_ln = log(y)
par(mfrow = c(1, 2))
plot(x, y, main = "y")
plot(x, y_ln, main = "log(y)")
# Строим матрцу X
X = cbind(rep(1, n), log(x), x, x^2, x^3)
XtX = t(X) %*% X
det(XtX)
## [1] 1230563
coef = solve(XtX) %*% t(X) %*% y_ln
print(coef)
## [,1]
## 0.5865
## 2.3796
## x 3.4244
## 2.5434
## -1.0031
fit = lm(y_ln ~ I(log(x)) + x + I(x^2) + I(x^3))
summary(fit)
##
## Call:
## lm(formula = y_ln ~ I(log(x)) + x + I(x^2) + I(x^3))
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.692 -1.159 0.244 1.085 4.620
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.587 2.179 0.27 0.79
## I(log(x)) 2.380 0.557 4.27 4.6e-05 ***
## x 3.424 5.838 0.59 0.56
## I(x^2) 2.543 5.286 0.48 0.63
## I(x^3) -1.003 1.564 -0.64 0.52
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.95 on 95 degrees of freedom
## Multiple R-squared: 0.894, Adjusted R-squared: 0.889
## F-statistic: 199 on 4 and 95 DF, p-value: <2e-16
error = predict(fit) - y_ln
sum(error)
## [1] 7.405e-14
plot(error, main = "error")
plot(x, y_ln)
points(x, predict(fit), col = "red")
