Hallo allerseits,
Wie berechne ich die Wurzeln einer Funktion in JSL für Konfidenzintervalle?
Ich habe einen Sonderfall, in dem ich Konfidenzintervalle berechnen muss, die nicht in JMP implementiert sind.
Unter Verwendung des Wilks-Theorems möchte ich Konfidenzintervalle berechnen, indem ich die Wurzeln der folgenden Funktion berechne:
f(x)=-2log(Likelihood(x|Data))-ChiSquare Quantile( 1-alpha, 1).
Wobei x eine reelle Zahl ist und die Daten eine univariate Stichprobe sind
Gibt es eine eingebaute JSL-Funktion, die ich verwenden kann? Oder zumindest etwas ähnliches?
Ich konnte nichts im Web oder im Scripting-Index finden.
Über ein paar Anregungen würde ich mich sehr freuen.
Danke im Voraus
Hallo, @olithun -san, @rcast15 -san und @DunnTiger764 -san,
Ich bin Yusuke Ono, Senior Tester bei JMP Japan. Ich entschuldige mich im Voraus für mein schlechtes Englisch.
Wenn Sie das Profil-Likelihood-Konfidenzintervall in einem allgemeinen Modell berechnen möchten, das durch die Maximum-Likelihood-Schätzung geschätzt wird, möchte ich Ihnen die Verwendung der nichtlinearen Plattform empfehlen.
Um diese Funktion in der nichtlinearen Plattform zu nutzen, müssen Sie die negative Log-Likelihood selbst als Verlustfunktion angeben. Sobald Sie die Verlustfunktion in einer Datentabelle angeben, berechnet die nichtlineare Plattform automatisch die Profil-Likelihood-Vertrauensgrenzen, indem Sie einfach auf die Schaltfläche [Vertrauensgrenzen] klicken.
[Schritte]
(1) Geben Sie die Verlustformel in einer Spalte einer JMP-Datentabelle an. Diese Verlustfunktion ist eine negative Log-Likelihood (= -LogLikelihood(theta|data)), und die Modellparameter werden in der Formel als Parameter definiert.
(2) Wählen Sie [Analysieren] > [Spezialisierte Modellierung] > [Nichtlinear] .
(3) Wählen Sie die Spalte aus, die die negative Log-Likelihood-Funktion für [Verlust] enthält.
Klicken Sie auf [OK] .
(4) Klicken Sie im Bericht auf die Schaltfläche [Los] .
(5) Ändern Sie den Alpha-Wert, falls erforderlich.
Ändern Sie bei Bedarf das Konvergenzkriterium.
Click
[Vertrauensgrenzen] button.
Das folgende Beispiel dient zur Berechnung des Profil-Likelihood-Konfidenzintervalls für den Wahrscheinlichkeitsparameter in einer Einstichproben-Binomialverteilung. In der Formel lasse ich den konstanten Term (Log(_nC_k)) weg.
Beachten Sie, dass das Score-Intervall (die in der Distributionsplattform implementierte Methode) in diesem Beispiel statistisch besser sein muss als das Profil-Likelihood-Intervall. Dieses Beispiel zeigt nur, wie die nichtlineare Plattform zur Berechnung des Profil-Likelihood-Intervalls verwendet wird.
dt = New Table( "Binomial Data",
New Column( "k", Values( [10] )),
New Column( "n", Values( [500] )),
New Column( "Negative Loglikehood",
Formula( Parameter( {p = 0.5}, -:k * Log( p ) - (:n - :k) * Log( 1 - p ) ) ),
)
);
dt << Nonlinear(
Loss( :Negative Loglikehood ),
CL Limit( 1e-15 ),
Loss is Neg LogLikelihood( 1 ),
Newton,
Finish,
Plot( 0 ),
Confidence Limits
);
Der folgende Code überprüft, ob die vom obigen Skript berechnete Obergrenze tatsächlich unseren Erwartungen entspricht.
p0 = 0.0348142211885808;
p1 = 0.02;
Show(ChisquareQuantile(0.95,1));
Show(2 *(Log(Binomial Probability(p1,500,10)) - Log(Binomial Probability(p0,500,10))));
/*
ChiSquare Quantile(0.95, 1) = 3.84145882069412;
2 * (Log(Binomial Probability(p1, 500, 10)) - Log(Binomial Probability(p0, 500, 10))) = 3.84145882069407;
*/
Die nichtlineare Plattform kann auch die Profil-Likelihood-Konfidenzintervalle berechnen, wenn Ihr Modell über mehr als einen Modellparameter verfügt.
Beachten Sie auch, dass einige Plattformen in JMP und JMP Pro die Konfidenzintervalle mit der Profilwahrscheinlichkeitsmethode berechnen.
Wenn Sie die Wurzel einer allgemeinen Gleichung berechnen möchten, hängt die Lösung der Gleichung in JMP meiner Meinung nach von der Komplexität der allgemeinen Gleichung ab. JMP bietet keinen allgemeinen Löser für allgemeine Gleichungen. Der Profiler kann in manchen Situationen nützlich sein. Sie erhalten möglicherweise eine vernünftige Lösung, indem Sie (f(x))^2 als Zielfunktion angeben und die JSL-Funktion Minimize() oder die nichtlineare Plattform verwenden. Und ich denke, dass es nicht sehr schwer, wenn auch nicht sehr einfach ist, mit der JMP Scripting Language (JSL) einen spezifischen Code für die Regula-Falsi-Methode, die Ridders-Methode oder die Newton-Methode zu programmieren. Wenn Ihre allgemeine Gleichung sehr kompliziert ist, ist es möglicherweise besser, den Python-Code von JMP zu verwenden. In JMP 18 können Sie Python von JMP aus sehr reibungslos verwenden.
Das Folgende ist ein sehr grobes Beispiel zur Programmierung der Regula-Falsi-Methode und der Ridder-Methode. Sie können Ihre JSL-Funktion auch als Argument für Ihre andere JSL-Funktion verwenden. Der folgende Code berechnet das Profil-Likelihood-Konfidenzintervall für den Wahrscheinlichkeitsparameter in der Binomialverteilung mit k = 10 und n = 500.
Dieses Skript ist nur ein sehr grobes Beispiel, um zu zeigen, wie schwierig es ist, die Algorithmen mit JSL zu programmieren. Ich habe die Details nicht überprüft und im folgenden Code gibt es keine Fehlerbehandlung.
Names Default To Here(1);
/* Regula Falsi method */
/* https://en.wikipedia.org/wiki/Regula_falsi */
Regula Falsi = Function({f, data, x0, x2, tol=1e-8, niter = 1000}, {Default Local},
fx0 = f(x0, data);
fx2 = f(x2, data);
If(fx0 * fx2>0, Write("f(pl)*f(pu) should be less than zero"); Return(.));
For(i = 1, i <= niter, i++,
If(Abs(fx0)
<TOL>
0, fx0=fx1; x0=x1, fx1*fx2>0, fx2=fx1; x2=x1, Return(x1) ); ); Return(.); ); /* Ridders' method */ /* https://en.wikipedia.org/wiki/Ridders%27_method */ Ridders = Function({f, data, x0, x2, tol=1e-8, niter = 1000}, {Default Local}, fx0 = f(x0, data); fx2 = f(x2, data); If(fx0 * fx2>0, Write("f(pl)*f(pu) should be less than zero"); Return(.)); For(i = 1, i <= niter, i++, If(Abs(fx0)
<TOL>
0,1,fx0<0,-1, Retrun(x0)); x3 = x1 + (x1-x0) * sign fx0 * fx1 /Sqrt(fx1^2 - fx0*fx2); fx3 = f(x3, data); If( fx1*fx3 < 0, x0 = x1; fx0 = fx1; x2 = x3; fx2 = fx3; , fx3*fx0>0, fx0=fx3; x0=x3, fx3*fx2>0, fx2=fx3; x2=x3, Return(x3) ); ); Return(.); ); /* Solve for x to satisfy f(x, data) = 0 */ f = Function({p, data}, {Default Local}, k = data[1]; n = data[2]; chi = data[3]; p1 = k/n; 2*(k*Log(p1/p)+(n-k)*Log((1-p1)/(1-p))) - chi ); k = 10; n = 500; chi = Chisquare Quantile(1-0.05, 1); data = k||n||chi; /* Upper Confidence Limit */ pl = k/n; pu = 0.99999; ci up rf= Regula Falsi(NameExpr(f), data, pl, pu); ci up = Ridders(NameExpr(f), data, pl, pu); /* Lower Confidence Limit */ pl = 0.00001; pu = k/n; ci low rf= Regula Falsi(NameExpr(f), data, pl, pu); ci low = Ridders(NameExpr(f), data, pl, pu); Show(ci low rf||ci up rf); Show(ci low || ci up);
</TOL>
</TOL>
Ich hoffe, dass diese Informationen auch nur ein bisschen nützlich sind.
