Salve a tutti,
Come faccio a calcolare le radici di una funzione in JSL per gli intervalli di confidenza?
Ho un caso speciale, in cui devo calcolare gli intervalli di confidenza che non sono implementati in JMP.
Usando il teorema di Wilks, voglio calcolare gli intervalli di confidenza calcolando le radici della seguente funzione:
f(x)=-2log(Likelihood(x|Data))-ChiSquare Quantile( 1-alpha, 1).
Dove x è un numero reale e Data è un campione univariato
C'è qualche build nella funzione JSL che posso usare? O almeno qualcosa di simile?
Non sono riuscito a trovare nulla sul Web o nell'indice di scripting.
Apprezzerei davvero alcuni suggerimenti.
Grazie in anticipo
Ciao, @olithun -san, @rcast15 -san e @DunnTiger764 -san,
Sono Yusuke Ono, Senior Tester presso JMP Giappone. Mi scuso in anticipo per il mio scarso inglese.
Se si desidera calcolare l'intervallo di confidenza della verosimiglianza del profilo in un modello generale stimato dalla stima della massima verosimiglianza, si consiglia di utilizzare la piattaforma Non lineare .
Per utilizzare questa funzionalità nella piattaforma Non bugiardo , è necessario specificare autonomamente la logaritmo-verosimiglianza negativa come funzione di perdita. Una volta specificata la funzione di perdita in una tabella dati, la piattaforma Non lineare calcola automaticamente i limiti di confidenza della verosimiglianza del profilo solo cliccando sul pulsante [Limiti di confidenza] .
[Passaggi]
(1) Specificare la formula di perdita in una colonna di una tabella dati JMP. Questa funzione di perdita è una logaritmo-verosimiglianza negativa (= -LogLikelihood(theta|data)) e i parametri del modello sono definiti come parametro nella formula.
(2) Selezionare [Analizza] > [Modellazione specializzata] > [Non lineare] .
(3) Selezionare la colonna che contiene la funzione di verosimiglianza logaritmica negativa per [Perdita] .
Fare clic su [OK] .
(4) Nel report, fare clic sul pulsante [Vai] .
(5) Se necessario, modificare il valore Alpha.
Se necessario, modificare il criterio di convergenza.
Click
[Limiti di confidenza] button.
Di seguito è riportato un esempio pratico per calcolare l'intervallo di confidenza della verosimiglianza del profilo per il parametro di probabilità nella distribuzione binomiale a un campione. Nella formula, elimino il termine costante (Log(_nC_k)).
Si noti che l'intervallo di punteggio (che è il metodo implementato nella piattaforma di distribuzione) deve essere statisticamente migliore dell'intervallo di verosimiglianza del profilo in questo esempio. Questo esempio mostra solo come utilizzare la piattaforma non lineare per calcolare l'intervallo di verosimiglianza del profilo.
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
);
Il codice seguente verifica che il limite superiore calcolato dallo script precedente sia effettivamente quello che ci aspettiamo.
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;
*/
La piattaforma non lineare può anche calcolare gli intervalli di confidenza della verosimiglianza del profilo quando il modello ha più di un parametro.
Si noti inoltre che alcune piattaforme in JMP e JMP Pro calcolano gli intervalli di confidenza con il metodo della verosimiglianza del profilo.
Se si desidera calcolare la radice di un'equazione generale, credo che il modo in cui risolverla su JMP dipenda dalla complessità dell'equazione stessa. JMP non offre un risolutore generale per equazioni generali. Il Profiler potrebbe essere utile in alcune situazioni. Si potrebbe ottenere una soluzione ragionevole specificando (f(x))^2 come funzione obiettivo e utilizzando la funzione JSL Minimize() o la piattaforma Nonlinear. E penso che non sia molto difficile, sebbene non facilissimo, programmare un codice specifico per il metodo Regula-Falsi, il metodo di Ridders o il metodo di Newton utilizzando il linguaggio di scripting JMP (JSL). Se l'equazione generale è molto complessa, potrebbe essere meglio utilizzare il codice Python di JMP. In JMP 18, è possibile utilizzare Python di JMP in modo molto fluido.
Quello che segue è un esempio molto approssimativo di programmazione del metodo Regula-Falsi e del metodo di Ridder. È anche possibile utilizzare la funzione JSL come argomento per l'altra funzione JSL. Il codice seguente calcola l'intervallo di confidenza della verosimiglianza del profilo per il parametro di probabilità nella distribuzione binomiale con k = 10 e n = 500.
Questo script è solo un esempio molto approssimativo per mostrare quanto sia difficile programmare gli algoritmi con JSL. Non ho verificato i dettagli e non c'è alcuna gestione degli errori nel codice seguente.
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>
Mi auguro che queste informazioni siano utili anche solo un po'.
