2013-01-15 40 views
6

Ich habe eine Frage, die vielleicht eher eine statistische Abfrage als eine direkt mit r ist, aber es kann sein, dass ich gerade ein r-Paket falsch aufrufen, so werde ich die Frage hier. Ich habe folgende Daten-Set:Stückweise lineare und nichtlineare Regression in R

x<-c(1e-08, 1.1e-08, 1.2e-08, 1.3e-08, 1.4e-08, 1.6e-08, 1.7e-08, 
1.9e-08, 2.1e-08, 2.3e-08, 2.6e-08, 2.8e-08, 3.1e-08, 3.5e-08, 
4.2e-08, 4.7e-08, 5.2e-08, 5.8e-08, 6.4e-08, 7.1e-08, 7.9e-08, 
8.8e-08, 9.8e-08, 1.1e-07, 1.23e-07, 1.38e-07, 1.55e-07, 1.76e-07, 
1.98e-07, 2.26e-07, 2.58e-07, 2.95e-07, 3.25e-07, 3.75e-07, 4.25e-07, 
4.75e-07, 5.4e-07, 6.15e-07, 6.75e-07, 7.5e-07, 9e-07, 1.15e-06, 
1.45e-06, 1.8e-06, 2.25e-06, 2.75e-06, 3.25e-06, 3.75e-06, 4.5e-06, 
5.75e-06, 7e-06, 8e-06, 9.25e-06, 1.125e-05, 1.375e-05, 1.625e-05, 
1.875e-05, 2.25e-05, 2.75e-05, 3.1e-05) 

y2<-c(-0.169718017273307, 7.28508517630734, 71.6802510299446, 164.637259265704, 
322.02901173786, 522.719633360006, 631.977073772459, 792.321270345847, 
971.810607095548, 1132.27551798986, 1321.01923840546, 1445.33152600664, 
1568.14204073109, 1724.30089942149, 1866.79717333592, 1960.12465709003, 
2028.46548012508, 2103.16027631327, 2184.10965255236, 2297.53360080873, 
2406.98288043262, 2502.95194879366, 2565.31085776325, 2542.7485752473, 
2499.42610084412, 2257.31567571328, 2150.92120390084, 1998.13356362596, 
1990.25434682546, 2101.21333152526, 2211.08405955931, 1335.27559108724, 
381.326449703455, 430.9020598199, 291.370887491989, 219.580548355043, 
238.708972427248, 175.583544448326, 106.057481792519, 59.8876372379487, 
26.965143266819, 10.2965349811467, 5.07812046132922, 3.19125838983254, 
0.788251933518549, 1.67980552001939, 1.97695007279929, 0.770663673279958, 
0.209216903989619, 0.0117903221723813, 0.000974437796492681, 
0.000668823762763647, 0.000545308757270207, 0.000490042305650751, 
0.000468780182460397, 0.000322977916070751, 0.000195423690538495, 
0.000175847622407421, 0.000135771259866332, 9.15607623591363e-05) 

, die, wenn Handlung wie folgt aussieht: Segmentation test http://i48.tinypic.com/25pltoy.png

Ich habe dann die Segmentierung Paket zu verwenden versuchte, in drei Regionen, die drei lineare Regressionen (durchgezogene schwarze Linie) zu erzeugen (10^⁻8-10^77,10^⁻7--10^⁻6 und> 10^-6), da ich eine theoretische Basis habe, um verschiedene Beziehungen in diesen verschiedenen Regionen zu finden. Offensichtlich aber mein Versuch, mit dem folgenden Code war nicht erfolgreich:

lin.mod <- lm(y2~x) 
segmented.mod <- segmented(lin.mod, seg.Z = ~x, psi=c(0.0000001,0.000001)) 

So meine erste Frage- gibt es weitere Parameter der Segmentierung I andere zwicken können als die Stützpunkte? Soweit ich es verstehe, habe ich hier Iterationen als Voreinstellung eingestellt.

Meine zweite Frage ist: Könnte ich vielleicht eine Segmentierung mit dem Paket nls versuchen? Es sieht so aus, als ob die ersten zwei Regionen auf der Zeichnung (10^& supmin; & sup8; - 10^& sup7; und 10^-7-10 & sup6; -6) weiter von dem linearen als dem letzten Abschnitt entfernt sind, so daß vielleicht eine Polynomfunktion besser wäre Hier?

Als ein Beispiel für ein Ergebnis finde ich akzeptabel Ich habe die ursprüngliche Handlung von Hand notiert: Annotated segmentation example http://i45.tinypic.com/zjb439.jpg.

Bearbeiten: Der Grund für die Verwendung von linearen Anpassungen ist die Einfachheit, die sie bieten, für mein untrainiertes Auge würde es eine ziemlich komplexe nichtlineare Funktion erfordern, um den Datensatz als eine Einheit zurückzufahren. Ein Gedanke, der mir in den Sinn kam, war ein lognormales Modell an die Daten anzupassen, da dies angesichts der Schräglage entlang einer logarithmischen x-Achse funktionieren könnte. Ich habe nicht genug Kompetenz in R, um dies zu tun, da mein Wissen sich nur auf fitdistr erstreckt, was, soweit ich es verstehe, hier nicht funktionieren würde.

Jede Hilfe oder Anleitung in einer relevanten Richtung würde am meisten geschätzt werden.

Antwort

4

Wenn Sie nicht mit dem segmented Paket zufrieden sind, können Sie versuchen, das earth Paket mit dem Mars Algorithmus. Aber hier finde ich, dass das Ergebnis des segmentierten Modells sehr akzeptabel ist. siehe unten den R-Squared.

lin.mod <- lm(y2~x) 
segmented.mod <- segmented(lin.mod, seg.Z = ~x, psi=c(0.0000001,0.000001)) 
summary(segmented.mod) 

Meaningful coefficients of the linear terms: 
       Estimate Std. Error t value Pr(>|t|)  
(Intercept) -2.163e+02 1.143e+02 -1.893 0.0637 . 
x   4.743e+10 3.799e+09 12.485 <2e-16 *** 
U1.x  -5.360e+10 3.824e+09 -14.017  NA  
U2.x   6.175e+09 4.414e+08 13.990  NA  

Residual standard error: 232.9 on 54 degrees of freedom 
Multiple R-Squared: 0.9468, Adjusted R-squared: 0.9419 

Convergence attained in 5 iterations with relative change 3.593324e-14 

Sie können das Ergebnis überprüfen, indem das Modell Plotten:

plot(segmented.mod) 

enter image description here

Um den Koeffizienten der Grundstücke zu erhalten, können Sie dies tun:

 intercept(segmented.mod) 
$x 
              Est. 
intercept1 -216.30 
intercept2 3061.00 
intercept3   46.93 

> slope(segmented.mod) 
$x 
             Est.   St.Err.  t value  CI(95%).l  CI(95%).u 
slope1  4.743e+10 3.799e+09  12.4800  3.981e+10  5.504e+10 
slope2 -6.177e+09 4.414e+08 -14.0000 -7.062e+09 -5.293e+09 
slope3 -2.534e+06 5.396e+06  -0.4695 -1.335e+07  8.285e+06 
+0

Ah seltsam dass ich es zum ersten Mal nicht zur Arbeit bringen konnte. Danke für die Tipps und die zusätzlichen Informationen zum Erhalten der Koeffizienten usw. Aus Interesse wissen Sie, wie ich zu etwas wie Lognormal zu solchen Daten passen könnte? Dies ist wahrscheinlich eine separate Frage, aber für den Fall, dass es einen schnellen Workaround gibt, dachte ich, ich würde es hier als Kommentar hinzufügen. Danke noch einmal. – user1912925

Verwandte Themen