Salut tout le monde,
Comment calculer les racines d'une fonction en JSL pour les intervalles de confiance ?
J'ai un cas particulier, où j'ai besoin de calculer des intervalles de confiance qui ne sont pas implémentés dans JMP.
En utilisant le théorème de Wilks, je souhaite calculer des intervalles de confiance en calculant les racines de la fonction suivante :
f(x)=-2log(Probabilité(x|Données))-Quantile Chicarré( 1-alpha, 1).
Où x est un nombre réel et les données sont un échantillon univarié
Existe-t-il une fonction JSL intégrée que je peux utiliser ? Ou au moins quelque chose de similaire?
Je n'ai rien trouvé sur le Web ou dans l'index de script.
J'apprécierais vraiment quelques suggestions.
Merci d'avance
Bonjour, @olithun -san, @rcast15 -san et @DunnTiger764 -san,
Je suis Yusuke Ono, testeur senior chez JMP Japon. Je tiens à vous informer par avance de mon faible niveau d'anglais.
Si vous souhaitez calculer l'intervalle de confiance de vraisemblance du profil dans un modèle général estimé par l'estimation du maximum de vraisemblance, je vous recommande d'utiliser la plateforme Non linéaire .
Pour utiliser cette fonctionnalité dans la plateforme Non-linéaire , vous devez spécifier vous-même la log-vraisemblance négative comme fonction de perte. Une fois la fonction de perte spécifiée dans une table de données, la plateforme Non-linéaire calcule automatiquement les limites de confiance de vraisemblance du profil en cliquant simplement sur le bouton [Limites de confiance] .
[Mesures]
(1) Spécifiez la formule de perte dans une colonne d'une table de données JMP. Cette fonction de perte est une log-vraisemblance négative (= -LogLikelihood(theta|data)), et les paramètres du modèle sont définis comme paramètre dans la formule.
(2) Sélectionnez [Analyser] > [Modélisation spécialisée] > [Non linéaire] .
(3) Sélectionnez la colonne qui contient la fonction de vraisemblance logarithmique négative pour [Perte] .
Cliquez sur [D'ACCORD] .
(4) Dans le rapport, cliquez sur le bouton [Aller] .
(5) Modifiez la valeur Alpha si nécessaire.
Modifiez le critère de convergence si nécessaire.
Click
[Limites de confiance] button.
Voici un exemple pratique de calcul de l'intervalle de confiance de vraisemblance du profil pour le paramètre de probabilité dans une distribution binomiale à un échantillon. Dans la formule, je supprime le terme constant (Log(_nC_k)).
Notez que l'intervalle de score (méthode implémentée dans la plateforme de distribution) doit être statistiquement supérieur à l'intervalle de vraisemblance du profil dans cet exemple. Cet exemple montre uniquement comment utiliser la plateforme non linéaire pour calculer l'intervalle de vraisemblance du profil.
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
);
Le code suivant vérifie que la limite supérieure calculée par le script ci-dessus est vraiment ce que nous attendons.
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 plateforme Non linéaire peut également calculer les intervalles de confiance de vraisemblance du profil lorsque votre modèle comporte plusieurs paramètres de modèle.
Notez également que certaines plateformes dans JMP et JMP Pro calculent les intervalles de confiance par la méthode de vraisemblance du profil.
Si vous souhaitez calculer la racine d'une équation générale, la résolution de l'équation sur JMP dépend de sa complexité. JMP ne propose pas de solveur général pour les équations générales. Le profileur peut être utile dans certaines situations. Vous pouvez obtenir une solution raisonnable en spécifiant (f(x))^2 comme fonction objectif et en utilisant la fonction JSL Minimize() ou la plateforme non linéaire. De plus, je pense qu'il n'est pas très difficile, bien que difficile, de programmer un code spécifique pour la méthode de Regula-Falsi, la méthode de Ridders ou la méthode de Newton en utilisant le langage de script JMP (JSL). Si votre équation générale est très complexe, il peut être préférable d'utiliser le code Python de JMP. Dans JMP 18, vous pouvez utiliser Python de JMP sans problème.
Voici un exemple très approximatif de programmation des méthodes Regula-Falsi et Ridder. Vous pouvez également utiliser votre fonction JSL comme argument pour une autre fonction JSL. Le code suivant calcule l'intervalle de confiance de vraisemblance du profil pour le paramètre de probabilité en distribution binomiale avec k = 10 et n = 500.
Ce script n'est qu'un exemple très approximatif pour illustrer la difficulté de programmer les algorithmes avec JSL. Je n'ai pas vérifié les détails et le code suivant ne contient aucune gestion des erreurs.
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>
J’aimerais que ces informations soient utiles, même un tout petit peu.
