Re: Marcum Q function
- To: mathgroup at smc.vnet.net
- Subject: [mg25567] Re: [mg25526] Marcum Q function
- From: BobHanlon at aol.com
- Date: Sat, 7 Oct 2000 03:36:15 -0400 (EDT)
- Sender: owner-wri-mathgroup at wolfram.com
In a message dated 10/6/2000 12:29:24 AM, a.coccoli at guest.cnuce.cnr.it writes:
>I am trying to evaluate the Marcum Q function with Mathematica, but I am
>finding some problems. I would like to know if you could suggest me the
>right way to do it. (any suggestion to on-line papers about this problem
>are well accepted too).
>
MarcumQ::"usage" =
"MarcumQ[a, b] is the Marcum Q-function. \n\
MarcumQ[m,\[NonBreakingSpace]a,\[NonBreakingSpace]b] is the generalized \
Marcum Q-function. \nMarcumQ[a,\[NonBreakingSpace]b]\[NonBreakingSpace]=\
\[NonBreakingSpace]MarcumQ[1,\[NonBreakingSpace]a,\[NonBreakingSpace]b]. \n\
[Harry L. Van Trees; Detection, Estimation, and Modulation Theory; \nWiley, \
New York; 1968; pp. 394-395, 411]";
Unprotect[MarcumQ];
MarcumQ::invparam = "Invalid parameter. `1` parameter must be `2`. Called
with parameter = `3`";
MarcumQ[m_:1, a_, b_] /; m < 1/2 := Message[MarcumQ::invparam, "m", ">= 1/2",
m];
MarcumQ[m_:1, a_, b_] /; Negative[a] := Message[MarcumQ::invparam, "a",
"nonnegative", a];
MarcumQ[m_:1, a_, b_] /; Negative[b] := Message[MarcumQ::invparam, "b",
"nonnegative", b];
MarcumQ[m_:1, a_, b_ /; b == 0] := 1;
MarcumQ[m_:1, a_, Infinity] := 0;
MarcumQ[m_:1, a_ /; a == 0, b_] := GammaRegularized[m, b^2/2];
MarcumQ[m_ /; m == 1, a_, b_] := MarcumQ[a, b];
MarcumQ[a_ /; NonNegative[a], a_] := 1/2*(1 + Exp[-a^2]*BesselI[0, a^2]);
MarcumQ/: MarcumQ[a_, b_] + MarcumQ[b_, a_] := 1 + Exp[-(1/2)*(a^2 +
b^2)]*BesselI[0, a*b];
MarcumQ[(a_)?NumericQ, (b_)?NumericQ] :=
NIntegrate[x*Exp[-(1/2)*(a^2 + x^2)]*BesselI[0, a*x], {x, b, Infinity}];
MarcumQ[(m_)?NumericQ, (a_)?NumericQ, (b_)?NumericQ] :=
NIntegrate[x*(x/a)^(m - 1)*Exp[-(1/2)*(a^2 + x^2)]*
BesselI[m - 1, a*x], {x, b, Infinity}];
Derivative/: Derivative[0, 0, 1][MarcumQ][m_, a_, b_] := -b^m*a^(1 -
m)*Exp[-(1/2)*(a^2 + b^2)]*
BesselI[m - 1, a*b];
Derivative/: Derivative[0, 1][MarcumQ][a_, b_] := -b*Exp[-(1/2)*(a^2 +
b^2)]*BesselI[0, a*b];
Protect[MarcumQ];
Table[Plot3D[MarcumQ[m, a, b], {a, 0, 4}, {b, 0, 5},
ViewPoint -> {-2.148, -2.156, 1.479}], {m, 1/2, 7/2, 3/2}];
Bob Hanlon