Fitting Fraktionen von multivariate Polynom

Question

Wie das multivariate fraktionierte Polynom der folgenden Form passen:

eine Funktion Gegeben: y = f(x,z), eine Funktion von zwei Variablen x und z. Genauer gesagt ist es der Form:

y = (x^2 + x^3)/(z^2 + z^3) 

Zähler ein Polynom eines dritten Grades eines Prädiktor x, und der Nenner ist auch ein Polynom eines dritten Grades einiger Prädiktor z.

Ich mag Polynome passen für jeden des Prädiktor x und z, also muß ich die Koeffizienten A, B, C finden, D:

y = (A*x^2 + B*x^3)/(C*z^2 + D*z^3) 

Grundsätzlich ist y ein Verhältnis von zwei Polynomen Grad 3. Wie passt man eine solche Funktion?

Beispiel eines Datenrahmens unten. Ich kann keinen vollständigen Datenrahmen posten, da er mehr als 1000 Zeilen enthält.

y = c(-4.10594369806545, -4.23691458506868, -4.24690667936422, -3.53677470254628, -4.30406509320417, -4.19442802077908, -4.66857169733859, -2.82271310942235, -4.19720194766181, 3.52164353473802, -4.3917019001973, -5.41654474791269, 2.87471821731616, -3.85922481986118, -4.25370811223789, -3.57887855889961, -5.33913936106829, -4.11775265312012, -2.89958841300109, -4.18661983833127) 

x = c(8.06526520889773, 9.39897529082673,9.07348918922699,7.5522372875608, 9.17294998275762,5.77455154554441, 9.2005930205213, 8.07309119969315, 7.42177579364465,8.18896686364888, 8.07868822922987, 8.50956416425175,9.71269017726113, 7.98378106897745, 7.69893619981345, 8.49576524400262, 8.02224091680654,8.25400859056484, 7.58171964012531, 8.35655484545343) 

z = c(2.56494935746154, 4.99043258677874, 4.43081679884331,3.66356164612965,4.53259949315326,1.79175946922805,4.23410650459726, 5.38449506278909,3.13549421592915,4.34380542185368, 3.43398720448515,2.77258872223978,6.94985645500077,3.97029191355212, 3.40119738166216,4.39444915467244,2.19722457733622,3.91202300542815,4.06044301054642, 3.87120101090789) 

dat = data.frame(cbind(y=y,x=x,z=z)) 

UPDATE:

Anruf zu nls:

nls(y~(a*(x**2) + b*(x**3))/(c*(z**2) + d*(z**3)), dat, start=list(a=1,b=1,c=1,d=1)) 

Answer 1

Dies ist ein schönes Problem, dass Sie hier sind. Wenigstens habe ich einiges daraus gelernt. Ich habe jedoch das Gefühl, dass sich diese Frage eher um die Lösung dieser spezifischen Aufgabe (Hochschulauftrag?) Dreht als um eine allgemeine Frage.

Aber lassen Sie mich zu der Lösung einen Ansatz teilen: Was wir hier haben

kann

Lösen von y^\ theta wird numerisch zu

vereinfachter werden überschaubar. Wie wir vielleicht sehen (und nach einem sehr harten Versuch, ein nicht lineares Problem nicht zu lösen), ist es tatsächlich eine Division von zwei linearen Modellen. Ein Ansatz besteht also darin, Koeffizienten der beiden Probleme getrennt zu schätzen. Das heißt, wir fixieren die Koeffizienten a und b, um c und d zu finden, und verwenden danach c und d, um a und b zu finden.

folgenden Code löst für Koeffizienten

Zunächst einige Daten

library(dplyr) 

sampleData <- data.frame(x = runif(100, -100, 100), z = runif(100, -100, 100)) %>% 
    mutate(y = ((-2 * x^2) + (5 * x^3))/(-4 * z^2 + 6 * z^3)) %>% 
    mutate(zxfactor = z^2/x^2, 
     yy = y * zxfactor) 

jetzt lösen wir für yy. Mit einigen zufälligen Startwerten ...

init2 <- structure(runif(4, -10, 10), names=c("A", "B", "C", "D")) 
coefab <- init2[c("A", "B")] 
coefcd <- init2[c("C", "D")] 

...wir brauchen ein lineares Modell für a und b von

und ein lineares Modell für c und d von

# don't use for loop but determine a terminal condition... but i'm too lazy :-) 
for(i in 1:100) { 
    # make linear prediction using coeff. c and d 
    sampleData <- sampleData %>% 
    mutate(yab = yy * (coefcd[1] + coefcd[2] * z)) 
    # and fit a model for a and b 
    coefab <- coef(lm(yab ~ x, sampleData)) 
    # then make a linear prediction using coeff. a and b 
    sampleData <- sampleData %>% 
    mutate(ycd = (coefab[1] + coefab[2] * x)/yy) 
    # and fit a model for c and d 
    coefcd <- coef(lm(ycd ~ z, sampleData)) 
} # repeat until satisfied 

coefab 
coefcd 

Sind zu passen glücklich, dass wir mit den Koeffizienten haben wir gefunden? Lets überprüfen:

optimFun <- function(params, out, x, z) { 
    res <- (params[1] + params[2]*x)/(params[3] + params[4]*z) 
    return(sqrt(sum((out - res)^2))) 
} 

optimFun(c(coefab, coefcd), x = sampleData$x, z = sampleData$z, out = sampleData$yy) 

> optimFun(c(coefab, coefcd), x = sampleData$x, z = sampleData$z, out = sampleData$yy) 
[1] 1.951043e-12 

Tatsächlich wir da die Differenz zwischen den modellierten Funktion Schätzungen und Daten yy sind (die skalierte eins) nahe bei Null liegt. Verschiedene Iterationen führen zu unterschiedlichen Parameterschätzungen, da das Problem überbestimmt ist. (Vielleicht hat jemand kann es näher erläutern)

---

Kommentare:

• Die Schätzungen sind näher an null als bestimmt Schätzungen nls
• Wenn Sie optim verwenden Sie die optimFun verwenden können . Und für schnellere Konvergenz können Sie sogar definieren, Derivat-Funktion
• Sein sehr schönes Problem zu zeigen, wo allgemeine Optimierung fehlschlagen kann und dass es immer wert ist, über das vorliegende Problem nachzudenken.
• Versuchen Sie lm(yyz ~ x, data = sampleData %>% mutate(yyz = yy*(-4 + 6*z))) , die genaue Werte für a = -2 und b = 5 zurückgibt. Bei zwei beliebigen Parametern können Sie ein passendes Paar finden, das die Funktion minimiert.

Answer 2

Ich wurde in der gleichen Richtung mit Drey Denken, durch Datentransformation: u = (x/z)^2, so dass y ~ (a * x^2 + b * x^3)/(c * z^2 + d * z^3)y ~ u * (a + b * x)/(c + d * z) wird.

Dann machte Drey einen Schritt vorwärts, durch weitere Transformation: v = y/u, so dass das Problem v ~ (a + b * x)/(c + d * z) wird.

Hinweis, wie ich bereits sagte, ist es nicht legitim alle a, b, c und d in Ihrem Modell zu haben, da dann das Modell nicht identifizierbar wird. Um dieses Problem zu verstehen, beachten Sie, dass das Ergebnis gleich ist, wenn wir alle Parameter um einen Faktor von k aufblasen. Mit anderen Worten, die Lösung ist nicht einzigartig. Ein üblicher Weg, um damit fertig zu werden, besteht darin, einen dieser Parameter zu verwerfen. Ich entscheide mich, die d aufzugeben, also sind wir Ziel, ein Modell zu passen: v ~ (a + b * x)/(c + z).

Keine Sorge! Wir verlieren hier keine Informationen. In der Tat, wenn wir Zähler und Nenner durch d in der ursprünglichen Form teilen, erhalten wir die letztere Form.

Nun lassen Sie uns nls nennen:

dat$u <- with(dat, (x/z)^2) 
dat$v <- with(dat, y/u) 

fit <- nls(v ~ (a + b * x)/(c + z), dat, start = list(a = 1, b = 1, c = 1)) 

#Nonlinear regression model 
# model: v ~ (a + b * x)/(c + z) 
# data: dat 
#  a  b  c 
#-3.8316 -0.2824 5.8768 
# residual sum-of-squares: 8.692 

#Number of iterations to convergence: 20 
#Achieved convergence tolerance: 6.138e-06