大家好,
如何在 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>
我希望這些資訊有一點用處。