大家好,
如何在 JSL 中计算置信区间的函数根?
我有一个特殊情况,我需要计算 JMP 中未实现的置信区间。
使用威尔克斯定理,我想通过计算以下函数的根来计算置信区间:
f(x)=-2log(Likelihood(x|Data))-ChiSquare Quantile(1-alpha, 1)。
其中 x 是实数,数据是单变量样本
我可以使用任何内置的 JSL 函数吗? 或者至少是类似的东西?
我在网上或脚本索引中找不到任何东西。
我真的很感激一些建议。
提前致谢
为什么不将此建议添加到愿望清单讨论中? 首先搜索以查看是否其他人已经建议并投票支持它(kudo)。 如果没有其他人有,那么添加它以供考虑。 谢谢!
Minimize(Abs(function)) 在某些情况下可能表现不佳。JMP 应该真正添加一个带括号的根求解器,就像 R 中的 uniroot() 一样。在此之前,Riders 的方法 ( Ridders 方法 - 维基百科) 可以在jsl中编程。它几乎和布伦特的方法一样有效(布伦特的方法 - 维基百科),但更容易编程。数字食谱包括一个实现。
感谢你的回答,
Minimize 函数在极端情况下似乎有点不稳定。
但是,您的提示为我提供了一个很好的起点来提出解决方案。
看帮助有关最小化()功能。 你的表情是 Abs( -2log(Likelihood(x|Data))-卡方分位数(1-alpha, 1) )。 这种数值优化应该能够求解根。 最简单的方法是使用脚本编辑器。 选择文件 > 新建 > 脚本. 右键单击并选择显示嵌入式日志. 这样,函数调用的结果将显示在同一窗口中。
你好, @olithun -san、 @rcast15 -san 和@DunnTiger764 -san,
我是小野裕介(Yusuke Ono),JMP日本的高级测试员。我先承认我的英语不太好。
如果您想计算由最大似然估计估计的一般模型中的轮廓似然置信区间,我建议您使用非线性平台。
要在非线性平台中使用此功能,您需要自行指定负对数似然作为损失函数。在数据表中指定损失函数后,非线性平台只需点击[置信限度]按钮即可自动计算轮廓似然置信限度。
[步骤]
(1)在 JMP 数据表的某一列中指定损失公式。此损失函数为负对数似然函数 (= -LogLikelihood(theta|data)),模型参数定义为公式中的参数。
(2)选择[分析] > [专业建模] > [非线性] 。
(3)选择包含负对数似然函数的列作为[损失] 。
点击[确定] 。
(4)在报告中,点击[前往]按钮。
(5)如果需要,更改 Alpha 值。
如果需要,更改收敛标准。
Click
[置信限度] button.
以下是一个计算单样本二项分布中概率参数的轮廓似然置信区间的小例子。公式中我去掉了常数项 (Log(_nC_k))。
请注意,本例中的得分区间(即“分布”平台中实现的方法)在统计上必须优于轮廓似然区间。本例仅演示如何使用“非线性”平台计算轮廓似然区间。
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
);
下面的代码检查上面的脚本计算出的上限是否确实是我们所期望的。
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;
*/
当您的模型具有多个模型参数时,非线性优化平台还可以计算概况似然置信区间。
还要注意,JMP 和 JMP Pro 中的一些平台通过概况似然法计算置信区间。
如果您想计算一个一般方程的根,我认为如何在 JMP 上求解该方程取决于该方程的复杂程度。JMP 不提供一般方程的通用求解器。Profiler 在某些情况下可能会有用。您可以通过指定 (f(x))^2 作为目标函数并使用 Minimize() JSL 函数或非线性平台来获得合理的解。而且,我认为使用 JMP 脚本语言 (JSL) 编写 Regula-Falsi 方法、Ridders 方法或牛顿法的具体代码虽然不太容易,但并不难。如果您的一般方程非常复杂,最好使用 JMP 中的 Python 代码。在 JMP 18 中,您可以非常流畅地从 JMP 使用 Python。
以下是一个非常粗略的 Regula-Falsi 方法和 Ridder 方法的编程示例。您也可以将 JSL 函数用作其他 JSL 函数的参数。以下代码计算二项分布中 k = 10 和 n = 500 的概率参数的剖面似然置信区间。
这个脚本只是一个粗略的示例,用来展示用 JSL 编写算法有多难。我没有检查细节,而且下面的代码中也没有错误处理。
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>
我希望这些信息有一点点用处。