皆さんこんにちは、
信頼区間のJSLで関数の根を計算するにはどうすればよいですか?
JMPに実装されていない信頼区間を計算する必要がある特別なケースがあります。
Wilks の定理を使用して、次の関数の根を計算することで信頼区間を計算します。
f(x)=-2log(尤度(x|データ))-ChiSquare Quantile(1-α, 1).
x は実数で、データは単変量サンプルです
使用できるビルドインJSL関数はありますか? または、少なくとも似たようなものですか？
Web やスクリプト インデックスで何も見つかりませんでした。
いくつかの提案をいただければ幸いです。
前もって感謝します
この提案をウィッシュ リストのディスカッションに追加してみませんか? 最初に検索して、他の誰かが既に提案しているかどうかを確認し、投票してください (kudo)。 誰も持っていない場合は、検討のために追加してください。 ありがとう！
Minimize(Abs(function)) は、場合によっては不適切な動作をすることがあります。JMP は、R の uniroot() の行に沿って、ブラケット ルート ソルバーを実際に追加する必要があります。それまではリダースの方法(リダース法 - ウィキペディア) はjslでプログラムできます。ブレントの方法とほぼ同じくらい効率的です (ブレント法 - ウィキペディア）、しかしプログラムはより簡単です。Numerical Recipes には実装が含まれています。
答えてくれてありがとう、
最小化機能は極端な状況では少し不安定なようです。
ただし、あなたのヒントは、解決策を考え出すための良い出発点を与えてくれました。
見るヘルプ関して最小化()関数。 あなたの式は Abs( -2log(Likelihood(x|Data))-ChiSquare Quantile(1-アルファ, 1) ）。 この数値最適化は、根を解くことができるはずです。 最も簡単な方法は、スクリプト エディタを使用することです。 選択するファイル > 新規 > スクリプト. 右クリックして選択埋め込みログを表示. このようにして、関数呼び出しの結果が同じウィンドウに表示されます。
こんにちは、 @オリサンさん、 @rcast15さん、 @ダンタイガー764さん、
JMPジャパンのシニアテスター、小野裕介です。英語が下手なことをあらかじめお詫び申し上げます。
最大尤度推定によって推定された一般的なモデルにおけるプロファイル尤度信頼区間を計算する場合は、非線形プラットフォームを使用することをお勧めします。
非線形プラットフォームでこの機能を使用するには、損失関数として負の対数尤度を自分で指定する必要があります。データテーブルで損失関数を指定すると、 [信頼限界]ボタンをクリックするだけで、非線形プラットフォームはプロファイル尤度の信頼限界を自動的に計算します。
[手順]
(1) JMPデータテーブルの列に損失関数を指定します。この損失関数は負の対数尤度（= -LogLikelihood(theta|data)）であり、モデルパラメータは式中のパラメータとして定義されます。
（２） [解析] > [特殊なモデリング] > [非線形]を選択します。
（３） [損失]に対する負の対数尤度関数を含む列を選択します。
[わかりました]をクリックします。
（4）レポートで[行く]ボタンをクリックします。
（５）必要に応じてアルファ値を変更します。
必要に応じて収束基準を変更します。
Click
[信頼限界] button.
以下は、1標本二項分布における確率パラメータのプロファイル尤度信頼区間を計算する簡単な例です。式では定数項(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には一般的な方程式を解くためのソルバーは用意されていません。状況によっては、プロファイラーが役立つかもしれません。目的関数に(f(x))^2を指定してMinimize() JSL関数や非線形プラットフォームを使うことで、妥当な解が得られるかもしれません。また、JMPスクリプト言語(JSL)を使って、Regula-Falsi法、Ridders法、Newton法などの具体的なコードをプログラムするのは、簡単ではないものの、それほど難しくはないと思います。もし、一般的な方程式が非常に複雑な場合は、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>
この情報が少しでも役に立つことを期待したいと思います。
