Rapporto tra polinomi approssimazione

Question

Sto cercando di adattare un polinomio al mio set di dati, che assomiglia a quella (set di dati completo è alla fine del post):

La teoria predice che la formulazione della curva è :

che assomiglia a questo (per x tra 0 e 1):

Quando provo a fare un modello lineare in R facendo:

mod <- lm(y ~ poly(x, 2, raw=TRUE)/poly(x, 2)) 

ottengo il seguente curva:

che è molto diverso da quello che mi sarei aspettato. Hai idea di come adattare una nuova curva a questi dati in modo che sia simile a quella che predice la teoria? Inoltre, dovrebbe avere solo un minimo.

completi set di dati:

---

vettore di valori x:

x <- c(0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.10, 0.11, 0.12, 
0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.20, 0.21, 0.22, 0.23, 0.24, 0.25, 
0.26, 0.27, 0.28, 0.29, 0.30, 0.31, 0.32, 0.33, 0.34, 0.35, 0.36, 0.37, 0.38, 
0.39, 0.40, 0.41, 0.42, 0.43, 0.44, 0.45, 0.46, 0.47, 0.48, 0.49, 0.50, 0.51, 
0.52, 0.53, 0.54, 0.55, 0.56, 0.57, 0.58, 0.59, 0.60, 0.61, 0.62, 0.63, 0.64, 
0.65, 0.66, 0.67, 0.68, 0.69, 0.70, 0.71, 0.72, 0.73, 0.74, 0.75, 0.76, 0.77, 
0.78, 0.79, 0.80, 0.81, 0.82, 0.83, 0.84, 0.85, 0.86, 0.87, 0.88, 0.89, 0.90, 
0.91, 0.92, 0.93, 0.94, 0.95) 

vettore di valori y:

y <- c(4.104, 4.444, 4.432, 4.334, 4.285, 4.058, 3.901, 4.382, 
    4.258, 4.158, 3.688, 3.826, 3.724, 3.867, 3.811, 3.550, 3.736, 3.591, 
    3.566, 3.566, 3.518, 3.581, 3.505, 3.454, 3.529, 3.444, 3.501, 3.493, 
    3.362, 3.504, 3.365, 3.348, 3.371, 3.389, 3.506, 3.310, 3.578, 3.497, 
    3.302, 3.530, 3.593, 3.630, 3.420, 3.467, 3.656, 3.644, 3.715, 3.698, 
    3.807, 3.836, 3.826, 4.017, 3.942, 4.208, 3.959, 3.856, 4.157, 4.312, 
    4.349, 4.286, 4.483, 4.599, 4.395, 4.811, 4.887, 4.885, 5.286, 5.422, 
    5.527, 5.467, 5.749, 5.980, 6.242, 6.314, 6.587, 6.790, 7.183, 7.450, 
    7.487, 8.566, 7.946, 9.078, 9.308, 10.267, 10.738, 11.922, 12.178, 13.243, 
    15.627, 16.308, 19.246, 22.022, 25.223, 29.752) 

Answer 1

Utilizzare nls per adattarsi a un modello non lineare. Si noti che la formula del modello non è definita in modo univoco come visualizzata nella domanda poiché se moltiplichiamo tutti i coefficienti per qualsiasi numero, il risultato darà comunque le stesse previsioni. Per evitare ciò, dobbiamo fissare un coefficiente. Una prima prova utilizzava i coefficienti mostrati nella domanda come valori iniziali (tranne il fixing di uno) ma che non riusciva a far cadere C così come i coefficienti risultanti alimentati in un secondo adattamento con C = 1.

st <- list(a = 43, b = -14, c = 25, B = 18) 
fm <- nls(y ~ (a + b * x + c * x^2)/(9 + B * x), start = st) 
fm2 <- nls(y ~ (a + b * x + c * x^2)/(9 + B * x + C * x^2), start = c(coef(fm), C = 1)) 

plot(y ~ x) 
lines(fitted(fm2) ~ x, col = "red") 

(continua dopo la tabella)

Nota: Ecco un esempio di utilizzo nls2 per ottenere i valori iniziano con ricerca casuale. Supponiamo che i coefficienti siano compresi tra -50 e 50.

library(nls2) 

set.seed(123) # for reproducibility 
v <- c(a = 50, b = 50, c = 50, B = 50, C = 50) 
st0 <- as.data.frame(rbind(-v, v)) 
fm0 <- nls2(y ~ (a + b * x + c * x^2)/(9 + B * x + C * x^2), start = st0, 
    alg = "random", control = list(maxiter = 1000)) 

fm3 <- nls(y ~ (a + b * x + c * x^2)/(9 + B * x + C * x^2), st = coef(fm0)) 

Answer 2

Dal momento che si dispone già di una previsione teorica, non fare sembra aver bisogno di un nuovo modello, ed è davvero solo un compito di plotting:

png(); plot(y~x) 
lines(x,mod,col="blue") 
dev.off() 

Non si può pretendere lm per produrre una buona approssimazione di un problema non lineare. Il denominatore che coinvolge x in quell'espressione teoretica rende questo intrinsecamente non lineare.