2016-05-10 20 views
1

ich eine nlsfitting Aufgabe, die ich mit R. Mein erster Versuch machen wollte diese here und als @Roland zu tun, wies darauf hin,Failing Fitting mit nicht-linearen Anpassungsverfahren zu tun (nlsLM, nlxb und wrapnls)

"Der Punkt ist, dass komplexe Modelle schwierig zu montieren sind. Umso weniger, unterstützen die Daten das Modell, bis es unmöglich wird. Sie könnten es vielleicht anpassen, wenn Sie extrem gute Startwerte hätten."

kann ich mit @Roland zustimmen, aber wenn excel kann dies passend warum nicht R nicht tun kann?

Grundsätzlich kann diese Anpassung mit dem GRG Nonlinear Solver von Excel durchgeführt werden, aber der Prozess ist sehr zeitaufwendig und manchmal ist die Anpassung nicht gut. (Da es in der Realität viele Daten gibt).

Hier ist mein Beispiel data.frame. Ich mag jede set Gruppe mit dem Modell unterhalb,

set.seed(12345) 
    set =rep(rep(c("1","2","3","4"),each=21),times=1) 
    time=rep(c(10,seq(100,900,100),seq(1000,10000,1000),20000),times=1) 
    value <- replicate(1,c(replicate(4,sort(10^runif(21,-6,-3),decreasing=FALSE)))) 
    data_rep <- data.frame(time, value,set) 

    > head(data_rep) 
      # time  value set 
      #1  10 1.007882e-06 1 
      #2 100 1.269423e-06 1 
      #3 200 2.864973e-06 1 
      #4 300 3.155843e-06 1 
      #5 400 3.442633e-06 1 
      #6 500 9.446831e-06 1 
      *  *  *   * 

Versuch 1

ich schon eine Frage trouble-when-adding-3rd-fitting-parameter-in-nls

Grundsätzlich ist das Problem hier gepostet vorgesehen passen ist, dass ich wollte, passend zu tun in gruppierten Daten und macht eine Vorhersage basierend auf den Anpassungskoeffizienten.

verwendet I nlsLM von library(minpack.lm) I einen Fehler bekam

Fehler in nlsModel (Formel, mf, Start, WTS, upper): singulären Gradientenmatrix bei Anfangsparameterschätzer

das war Auf den ersten Blick waren der Modellfehler oder meine Startwerte laut @Roland nicht gut. Auf der anderen Seite könnte ich dieses Modell mit nur zwei Anpassungsparametern anpassen. Das Problem tritt auf, wenn ich der Anpassungsfunktion den Parameter third hinzufügen wollte.

Versuch 2

In diesem Beitrag trouble-when-adding-3rd-fitting-parameter-in-nls von @G folgen. Grothendieck Vorschlag, versuchte ich nlxb von nlmrt Paket und reparierte einen der Parameter d zu d=32 und die Anpassung wie folgt tun;

formula = value~Ps*(1-exp(-2*f*time*exp(-d)))*1/(sqrt(2*pi*sigma))*exp(-(d-d_ave)^2/(2*sigma))*d_step 

d_step <- 1 
f <- 1e9 
d <- 32  
library(plyr) 
library(nlmrt) 
get.coefs <- function(data_rep) { 
fit <- nlxb(formula , 
       data = data_rep, 
       start=c(d_ave=44,sigma=12,Ps=0.5), 
       lower=c(d_ave=25,sigma=2,Ps=0.5), 
       upper=c(d_ave=60,sigma=15,Ps=1), 
       trace=TRUE) 
} 

fit <- dlply(data_rep, c("set"), .fun = get.coefs) # Fit data grouped by "set" 

# > fit 
# $`1` 
# nlmrt class object: x 
# residual sumsquares = 1.474e-07 on 21 observations 
#  after 12 Jacobian and 13 function evaluations 
#  name   coeff   SE  tstat  pval  #gradient JSingval 
# d_ave   42.0126   NA   NA   NA #-7.082e-15 0.001733 
# sigma   12.8377   NA   NA   NA #2.408e-15 1.289e-19 
# Ps    0.973223   NA   NA   NA #9.33e-15 3.37e-20 
#  
# $`2` 
# nlmrt class object: x 
# residual sumsquares = 6.2664e-08 on 21 observations 
#  after 12 Jacobian and 13 function evaluations 
#  name   coeff   SE  tstat  pval  #gradient JSingval 
# d_ave    42.246   NA   NA   NA #-7.269e-15 0.001428 
# sigma   12.7429   NA   NA   NA #2.568e-15 3.098e-19 
# Ps    0.981517   NA   NA   NA #9.211e-15 2.746e-20 
#  
# $`3` 
# nlmrt class object: x 
# residual sumsquares = 1.773e-07 on 21 observations 
#  after 12 Jacobian and 13 function evaluations 
#  name   coeff   SE  tstat  pval  #gradient JSingval 
# d_ave    41.968   NA   NA   NA #-6.438e-15 0.001798 
# sigma   12.8561   NA   NA   NA #2.173e-15 2.414e-19 
# Ps    0.972988   NA   NA   NA #8.534e-15 5.922e-20 

# $`4` 
# nlmrt class object: x 
# residual sumsquares = 2.5219e-07 on 21 observations 
#  after 12 Jacobian and 13 function evaluations 
#  name   coeff   SE  tstat  pval  #gradient JSingval 
# d_ave   41.8532   NA   NA   NA #-4.454e-15 0.001976 
# sigma   12.9045   NA   NA   NA #1.474e-15 3.443e-19 
# Ps    0.974319   NA   NA   NA #5.987e-15 3.124e-20 

# attr(,"split_type") 
# [1] "data.frame" 
# attr(,"split_labels") 
#  set 
# 1 1 
# 2 2 
# 3 3 
# 4 4 

Fitting-Koeffizienten sind sinnvoll wolaa! Aber dieses Mal habe ich gemerkt, dass (@ Gr. Grothendieck auch später darauf hingewiesen hat) es unmöglich ist, neue Werte nach nlxb vorherzusagen (warum =? Ich weiß es nicht!))

predvals <- ldply(fit, .fun=predictvals, xvar="time", yvar="value",xrange=range(range)) # predict values 

:: Sie predictvals Funktion von here

Fehler in UseMethod ("vorhersagen") finden: keine anwendbare Methode für 'vorhersagen' angewendet auf ein Objekt der Klasse "nlmrt"

Es gibt keine! coef oder predict methods für "nlmrt" Klassenobjekte.

Versuch 3

Nach @G folgen. Grothendieck ein weiterer Vorschlag Als nächstes versuchte ich wrapnls von nlmrt.

weil in diesem Beitrag erklärte er, can-we-make-prediction-with-nlxb-from-nlmrt-package

„, weil nlmrt Paket wrapnls bietet, die nlmrt ausgeführt wird und dann nls so dass ein "nls" Objekt Ergebnisse und dann das Objekt mit allen "nls" Klasse verwendet werden Methoden.

Aus dem gleichen nlmrt Paket immer noch Probleme wie bei unter

Ich gebe auf plyr nach meinem ersten post zu verwenden, weil das Laden von plyr und dplyr mein Problem komplexer macht. Also bleibe ich bei dplyr und verwende stattdessen do Funktion.

library(dplyr) 
library(nlmrt) 
formula = value~Ps*(1-exp(-2*f*time*exp(-d)))*1/(sqrt(2*pi*sigma))*exp(-(d-d_ave)^2/(2*sigma))*d_step 

d_step <- 1 
f <- 1e9 
d <- 32 

dffit = data_rep %>% group_by(set) %>% 
    do(fit = wrapnls(formula , 
       data = ., 
       start=c(d_ave=44,sigma=12,Ps=0.5), 
       lower=c(d_ave=25,sigma=2,Ps=0.5), 
       upper=c(d_ave=60,sigma=15,Ps=1), 
       trace=TRUE)) 

Fehler in nlsModel (Formel, mf, Start, WTS, upper): singuläre Gradientenmatrix bei Anfangsparameterschätzer

I kehrte zurück, wo ich mit diesem Fehler gestartet. Ich denke, ich habe alles versucht, was ich tun kann, suche nach relevanten Beispielen (obwohl nur 3), lese das Buch und folge den Vorschlägen.

+0

Ich sagte, dass das 'nlmrt'-Paket keine 'predicate'-Methode für Objekte hat, die von' nlxb' erzeugt werden, sondern dass das Paket 'wrapnls' liefert. Das ist nicht dasselbe wie zu sagen, dass Voraussagen unmöglich sind. –

+0

@ G.Grothiedieke Es tut mir leid, wenn ich dich falsch verstehe. Ich meinte das mit 'nlxb' ist es unmöglich.Übrigens habe ich herausgefunden, dass ich die passenden Koeffizienten erhalten kann, nachdem ich 'coeff <- ldply (fit, funktion (mod) c (coef (mod))) gemacht habe' ist das sinnvoll für die Vorhersage nach 'nlxb'? – Alexander

Antwort

2

Verwenden nls2 vom NLS2 Paket nach nlxb wie folgt aus (unter der Annahme, data_rep, formula, d_step, f und d von der Frage). Um das Beispiel minimal wir dplyr haben zu machen, eliminiert und zeigen nur die Berechnung für Satz == 2.

library(nlmrt) 
library(nls2) 

data_rep2 <- subset(data_rep, set == 2) 

fit.nlxb <- nlxb(formula , data = data_rep2, 
       start = c(d_ave = 44, sigma = 12, Ps = 0.5), 
       lower = c(d_ave = 25, sigma = 2, Ps = 0.5), 
       upper = c(d_ave = 60, sigma = 15, Ps = 1)) 

fit.nls <- nls2(formula, data = data_rep2, start = fit.nlxb$coefficients, 
    algorithm = "brute-force") 

identical(fit.nlxb$coefficients, coef(fit.nls)) 
## [1] TRUE 

fit.nls ist ein "nls" Klasse-Objekt mit den gleichen Koeffizienten wie fit.nlxb und wir können fitted() und predict() und alle verwenden die anderen "nls" Methoden darauf.

+0

Das ist großartig. Ich verstehe, dass Sie nur das Ergebnis für 'set == 2' anzeigen, um alle passenden Coefs zu erhalten. Was soll ich tun? – Alexander

+1

Verwenden Sie 'lapply (1: 4, f)', wobei f output fit.nls für die entsprechende Teilmenge von data_rep ausgibt. –

+0

Ich werde es schätzen, wenn Sie dies zur Antwort implementieren können, da meine Hauptfrage in Gruppen passt. Ich konnte diese kurze Antwort nicht verstehen. – Alexander