Identifizierbarkeit in einem nichtlinearen Regressionsproblem

8

Angenommen, ich arbeite mit dem folgenden Modell

yi=α(1exp(βti))+γ(1exp(δti))+εi.

Das εi sind iid Gauß mit dem Mittelwert Null und ich versuche, die besten Anpassungswerte von zu finden α,β,γ,δ.

Nehmen wir zur Verdeutlichung an, dies ist ein Modell für die Gesamtmenge einiger Bakterienarten mit zwei Unterarten, die gemäß dem ersten und zweiten Term auf der RHS zeitlich wachsen, aber wir messen nur die Gesamtpopulation. Hinweis: Dies ist nicht die tatsächliche Einstellung, reicht jedoch für die Frage aus.

Das Modell ist im üblichen Sinne nicht identifizierbar, da ich immer nur tauschen kann α und γ zum Beispiel und erhalten genau die gleiche Dichte / Wahrscheinlichkeit.

Wie zu erwarten ist, habe ich, wenn ich ein MCMC auf diesem Gebiet betreibe, schrecklich breite Posterioren, und jeder nichtlineare Ansatz der kleinsten Quadrate reagiert unglaublich empfindlich auf anfängliche Vermutungen - wir haben ein großes Plateau in der Wahrscheinlichkeitsfunktion.

Ein besseres experimentelles Design ist derzeit keine Option - eine getrennte Messung der Unterarten wäre natürlich die beste Option.

Kann ich mit diesem Problem etwas anfangen oder ist ein besseres experimentelles Design die einzige Option?

Wille
quelle
1
Warum haben Sie zwei Begriffe mit derselben Struktur? Wie Sie sagten, sind sie falsch identifiziert.
Anna SdTC
1
Sie repräsentieren zwei "Teilgrößen" der gemessenen Größe. Ich habe dieses Modell nicht formuliert, ich habe es erhalten, um zu versuchen, es an einige Daten anzupassen, und die Leute, die es mir gegeben haben, sagten, sie wollen es nicht ändern. Ich wollte sicherstellen, dass ich nichts tun konnte, bevor ich ihnen sagte, dass es nicht möglich sei.
Will
1
Ich glaube nicht, dass Sie auf diesem Modell basieren können. Können Sie eine zweistufige Schätzung durchführen, bei der es im ersten Schritt nur einen dieser Exponentialterme gibt und im zweiten Schritt die "Subquantitäten" berechnen?
Anna SdTC
1
Das könnte funktionieren. Ich muss definitiv einige Parameter korrigieren, wenn ich überhaupt irgendwohin kommen will, und das könnte ein guter Weg sein, dies zu tun. Was ist, wenn ich entweder einen Exponentialprior (mit 'großem' Ratenparameter) vorlege?α oder γ(Sie sind nicht negativ)? Denken nach dem Vorbild, irgendwie zu regulieren ...
Will
1
Zur allgemeinen Information, als nichtlineare Methode der kleinsten Quadrate gibt es "Levenberg-Marquardt" (LM) und "Trust-Region Reflective" (TRR). Haben Sie eine davon ausprobiert? Ich weiß, dass Sie mit dem TRR Grenzen für Parameter setzen können, während LM dies nicht tut. Für Ihr Problem denke ich, dass Sie etwas mehr haben müssen, mit dem Sie beide Begriffe unterscheiden können. Entweder können Sie Messungen mit nur einer Population identifizieren und eine zweistufige Schätzung durchführen, wie Anna feststellte, oder Sie können Ihr Modell ändern, um die Spezifitäten jeder Unterart zu berücksichtigen. Ansonsten sehe ich die Nützlichkeit dieses Modells nicht.
Learn_and_Share

Antworten:

5

Es gibt kein Identifizierbarkeitsproblem, außer in dem trivialen Sinne, dass ein bestimmtes Modell zwei Beschreibungen haben kann. Das eigentliche Problem scheint die Schwierigkeit bei der Anpassung des Modells zu sein - dies liegt jedoch eher an der Parametrisierung der Modelle als an der mangelnden Identifizierbarkeit.

Dieses Problem hat eine ebenso triviale Lösung: Erklären Sie dies ohne Verlust der Allgemeinheit βδ. Wenn Sie wirklich pingelig sein wollen, bestehen Sie auch darauf, dass wennβ=δ, dann αγ.

Leider erfordert dies ein Verfahren, um das Modell anzupassen und diese Einschränkungen zu berücksichtigen. Das Einführen einer Einschränkung ist hier jedoch nicht so schlimm, da die Anwendung so ist, dass offensichtlich alle Parameter ohnehin nicht negativ sind: Der Parameterraum hat bereits scharfe Grenzen. Das Einfügen einer weiteren Einschränkung erzwingt keine Änderungen bei der Anpassung des Modells.


Eine bekannte Methode, um eine eingeschränkte Optimierung in eine nicht eingeschränkte umzuwandeln, besteht darin, das Problem neu zu parametrisieren, so dass im neuen Parameterraum die Grenzen ins Unendliche verschoben werden. Hier gibt es viele Möglichkeiten, dies zu erreichen. Eine Überlegung, was die Parameter bedeuten, wird uns leiten. Bestimmtes,ν=α+γ ist das von der Funktion erreichte Maximum

tg(t;α,β,γ,δ)=α(1eβt)+γ(1eδt)
zum t0. Gegebenνdann unbedingt 0αν und γ=να. Wenn sich nicht negative Werte zu einem festen Ganzen summieren, funktioniert es oft, ihre Proportionen des Ganzen in Bezug auf Winkel zu parametrisieren: Sei ein Anteil der quadratische Kosinus und der andere der quadratische Sinus. Darüber hinaus ein einfacher Weg, um sicherzustellenν, β, und δPositiv ist, sie exponentiell zu machen - das heißt, ihre Logarithmen als Parameter zu verwenden. Schließlich durchzusetzenδβ, einstellen δ der quadratische Kosinus einiger Winkelzeiten zu sein β. Daher können wir das Problem durch Anpassen der Funktion neu parametrisieren

tf(t;n,a,b,d)=en(1cos(a)2exp(ebt)sin(a)2exp(ebcos(d)2t)).

Aus Schätzungen dieser Parameter (die übrigens aufgrund der Mehrdeutigkeit der Winkel nicht "identifizierbar" sind a und d) können Sie die ursprünglichen als wiederherstellen

α=encos(a)2β=ebγ=ensin(a)2δ=ebcos(d)2.

Die Eigenschaften der Exponential- und Triggerfunktionen stellen sicher, dass alle Einschränkungen gelten: α>0, βδ>0, und γ>0. (Da Schwimmer mit doppelter Genauigkeit astronomisch klein werden können, gibt es keinen praktischen Unterschied zwischen> und in diesen Einschränkungen.)

In diesem genau definierten Sinne ist das Modell identifizierbar, obwohl die zur Anpassung verwendeten Parameter nicht identifizierbar sind.


Obwohl man MCMC verwenden könnte, ist es einfacher, einen numerischen Löser wie Newton-Raphson zu verwenden, wenn der Zweck nur darin besteht, die Kurve anzupassen. Der Trick besteht darin , einen guten Startwert zu finden . Das Maximum deryi wäre eine leichte Überschätzung von en;; Beginnen Sie also vielleicht mitn=log(max(yi)/2). Sie könnten mit beginnena=π/4Angenommen, jede Komponente leistet einen wesentlichen Beitrag zum Ganzen. Machen Sie einige vernünftige Vermutungen übereb und edbasierend auf erwarteten Zerfallsraten. Zum Beispiel, wenn der Bereich vont ist vernünftig, dann nimm b ein Bruchteil der größten sein t und vielleicht willkürlich auswählen d=π/4;; Verwenden Sie möglicherweise einen kleineren Startwert. ( Abhängig von diesen Auswahlmöglichkeiten erhalten Sie häufig unterschiedliche Werte für die Parameterschätzungen, die sich jedoch in der Regel nicht wesentlich auf die Funktion auswirkenfselbst .)

In vielen Fällen funktioniert dieser Ansatz auffallend gut. Außer wenn die Varianz der Fehler gleich groß ist wiemaxyi oder größer (wo es ohne eine große Datenmenge schwierig ist, ein Signal überhaupt zu erkennen), funktioniert die Anpassung auch mit winzigen Datenmengen: Alles, was benötigt wird, sind vier.

Beachten Sie, dass unabhängig von der Anpassung des Modells normalerweise große Unsicherheiten bei den Parametern bestehen: Diese Kurvenfamilie ist im Wesentlichen eine winzige Störung der Exponentialfamilie mit zwei Parametern tAeBt. In vielen Fällen also zwei der Parameter (entsprechend der AmplitudeA und längste Zerfallsrate B) können mit angemessener Genauigkeit identifiziert werden, aber die beiden anderen Parameter, die kleine Abweichungen von dieser Exponentialform widerspiegeln, sind normalerweise sehr unsicher.

Zahl

Die Abbildung zeigt ein Beispiel für eine herausfordernde Passform. Die zugrunde liegende Kurve ist schwarz dargestellt. Letztendlich erreicht es ein Maximum von4/3, sehr langsam. Nur24Datenpunkte sind verfügbar und als graue Punkte dargestellt. Die Standardabweichung der Zufallsfehler beträgt1/2ein beträchtlicher Anteil dieses Maximums. Viele der Fehler waren positiv, was dazu führte, dass die angepasste Kurve in Rot etwas höher war. Die beiden Exponentialkomponenten der angepassten Kurve sind als gestrichelte und gepunktete graue Linien dargestellt. Man zeigt einen raschen Anstieg auf eine Schwelle von1/3 Zu der Zeit t=1;; das andere spiegelt das andere Exponential wider, das bis zu seiner Schwelle von ansteigt1. (Sie werden wenig Hoffnung haben, diese scharfe "Schulter" in der Nähe zu reproduzierent=1 bis du eine hast 1000Datenpunkte oder mehr: Probieren Sie es aus, indem Sie nden folgenden Code variieren .)

Ihr Erfolg bei einem bestimmten Problem hängt von der Größe der Fehler ab. der Wertebereich vontdie abgetastet werden; wie diese Werte beabstandet sind; wie viele Werte sind verfügbar; und Wahl der Startwerte. Trotzdem scheint dies im Allgemeinen ein nachvollziehbares Problem zu sein, mit Lösungen, die schnell erhalten werden können. Darüber hinaus wird jeder Monteur mit maximaler Wahrscheinlichkeit ähnlich vorgehen, um die Summe der Quadrate der Residuen zu minimieren - und zusätzlich Konfidenzbereiche für die Parameter bereitstellen.

Dies ist der RCode, mit dem ich diesen Vorschlag getestet habe. Es gibt die Abbildung wieder und kann leicht geändert werden - ändern Sie die Werte der Variablen am Anfang -, um Daten zu untersuchen, die wie die von Ihnen möglicherweise vorhandenen aussehen.

#
# Describe the underlying model
#
set.seed(17)
alpha <- 1
beta <- 2
gamma <- 1/3
delta <- 1/10
sigma <- 1/2  # Error SD.
n <- 24
x.max <- 20   # Largest value of t.
#
# The original parameterization.
#
g <- function(x, alpha, beta, gamma, delta) {
  alpha * (1 - exp(-beta * x)) + gamma * (1 - exp(-delta * x))
}
#
# The re-parameterization.  `f.1` and `f.2` are the two exponential components.
#
f <- function(x, nu, t.a, log.b, t.d) {
  n <- exp(nu)
  a <- cos(t.a)^2
  alpha <- n*a
  gamma <- n*(1-a)
  beta <- exp(log.b)
  delta <- cos(t.d)^2 * beta
  n - alpha * exp(-beta * x) - gamma * exp(-delta * x)
}
f.1 <- function(x, nu, t.a, log.b, t.d) {
  n <- exp(nu)
  a <- cos(t.a)^2
  alpha <- n*a
  beta <- exp(log.b)
  alpha * (1 - exp(-beta * x))
}
f.2 <- function(x, nu, t.a, log.b, t.d) {
  n <- exp(nu)
  a <- cos(t.a)^2
  gamma <- n*(1-a)
  beta <- exp(log.b)
  delta <- cos(t.d)^2 * beta
  gamma * (1 - exp(-delta * x))
}
#
# The objective to minimize is the mean squared residual.
# This is equivalent to finding the MLE for Gaussian errors.
#
obj <- function(theta, x, y) {
  crossprod(y - f(x, theta[1], theta[2], theta[3], theta[4])) / length(x)
}
#
# Create data and plot them.
#
x <- seq(0, x.max, length.out=n)
y <- g(x, alpha, beta, gamma, delta) + rnorm(length(x), 0, sigma)
plot(x,y, pch=16, col="#00000040", xlab="t")
#
# Fit the curve.
#
theta <- c(nu=log(max(y)/2), t.a=pi/4, log.b=log(max(x)/10), t.d=pi/4)
fit <- nlm(obj, theta, x=x, y=y, gradtol=1e-14)
theta.hat <- fit$estimate
#
# Plot relevant curves.
#
curve(g(x, alpha, beta, gamma, delta), add=TRUE, lwd=2)
curve(f(x, theta.hat[1], theta.hat[2], theta.hat[3], theta.hat[4]),
      add=TRUE, col="Red", lwd=2)
curve(f.1(x, theta.hat[1], theta.hat[2], theta.hat[3], theta.hat[4]),
      add=TRUE, col="Gray", lty=2, lwd=2)
curve(f.2(x, theta.hat[1], theta.hat[2], theta.hat[3], theta.hat[4]),
      add=TRUE, col="Gray", lty=3, lwd=2)
whuber
quelle