Hola a todos,
¿Cómo calculo las raíces de una función en JSL para intervalos de confianza?
Tengo un caso especial, donde necesito calcular intervalos de confianza que no se implementan en JMP.
Usando el teorema de Wilks, quiero calcular los intervalos de confianza calculando las raíces de la siguiente función:
f (x) = - 2log (Probabilidad (x | Datos)) - ChiSquare Cuantil (1-alfa, 1).
Donde x es un número real y los datos son una muestra univariada
¿Hay alguna construcción en la función JSL que pueda usar? O al menos algo similar?
No pude encontrar nada en la web o en el índice de secuencias de comandos.
Realmente agradecería algunas sugerencias.
Gracias por adelantado
Gracias por la respuesta,
La función Minimizar parece ser un poco inestable en situaciones extremas.
Sin embargo, su consejo me dio un buen punto de partida para encontrar una solución.
Hola, @olithun -san, @rcast15 -san y @DunnTiger764 -san,
Soy Yusuke Ono, tester sénior de JMP Japón. Disculpen de antemano mi bajo nivel de inglés.
Si desea calcular el intervalo de confianza de verosimilitud del perfil en un modelo general estimado por la estimación de máxima verosimilitud, me gustaría recomendarle que utilice la plataforma No Lineal .
Para usar esta función en la plataforma No mentiroso , debe especificar la logaritmo-verosimilitud negativa como función de pérdida. Una vez especificada la función de pérdida en una tabla de datos, la plataforma No lineal calcula automáticamente los límites de confianza de la verosimilitud del perfil con solo hacer clic en el botón [Límites de confianza] .
[Pasos]
(1) Especifique la fórmula de pérdida en una columna de una tabla de datos de JMP. Esta función de pérdida es una log-verosimilitud negativa (= -LogLikelihood(theta|data)), y los parámetros del modelo se definen como Parámetro en la fórmula.
(2) Seleccione [Analizar] > [Modelado especializado] > [No lineal] .
(3) Seleccione la columna que contiene la función de log-verosimilitud negativa en [Pérdida] .
Haga clic en [Aceptar] .
(4) En el informe, haga clic en el botón [Ir] .
(5) Cambie el valor Alfa si es necesario.
Cambie el criterio de convergencia si es necesario.
Click
[Límites de confianza] button.
El siguiente es un ejemplo práctico para calcular el intervalo de confianza de verosimilitud del perfil para el parámetro de probabilidad en una distribución binomial de una muestra. En la fórmula, se omite el término constante (Log(_nC_k)).
Tenga en cuenta que el intervalo de puntuación (método implementado en la plataforma de distribución) debe ser estadísticamente mejor que el intervalo de verosimilitud del perfil en este ejemplo. Este ejemplo solo muestra cómo usar la plataforma no lineal para calcular el intervalo de verosimilitud del perfil.
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
);
El siguiente código verifica que el límite superior calculado por el script anterior sea realmente el que esperamos.
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 plataforma No lineal también puede calcular los intervalos de confianza de probabilidad del perfil cuando su modelo tiene más de un parámetro de modelo.
Tenga en cuenta también que algunas plataformas en JMP y JMP Pro calcularon los intervalos de confianza mediante el método de verosimilitud del perfil.
Si desea calcular la raíz de una ecuación general, creo que la forma de resolverla en JMP depende de su complejidad. JMP no ofrece un solucionador general para ecuaciones generales. El Profiler puede ser útil en algunas situaciones. Puede obtener una solución razonable especificando (f(x))^2 como función objetivo y utilizando la función Minimize() de JSL o la plataforma no lineal. Además, creo que no es muy difícil, aunque no muy fácil, programar un código específico del método Regula-Falsi, el método de Ridders o el método de Newton utilizando el lenguaje de scripting JMP (JSL). Si su ecuación general es muy compleja, puede ser mejor usar el código Python de JMP. En JMP 18, puede usar Python desde JMP sin problemas.
El siguiente es un ejemplo aproximado para programar el método Regula-Falsi y el método de Ridder. También puede usar su función JSL como argumento para otra función JSL. El siguiente código calcula el intervalo de confianza de la verosimilitud del perfil para el parámetro de probabilidad en la distribución binomial con k = 10 y n = 500.
Este script es solo un ejemplo preliminar para mostrar la dificultad de programar algoritmos con JSL. No he revisado los detalles y el siguiente código no incluye gestión de errores.
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>
Me gustaría esperar que esta información sea útil aunque sea un poquito.
