octave-signal 1.2.2-1build1 / usr / share / octave / packages / signal-1.2.2 / marcumq.m

## Copyright (C) 2012   Robert T. Short     <rtshort@ieee.org>
##
## This program is free software; you can redistribute it and/or modify it under
## the terms of the GNU General Public License as published by the Free Software
## Foundation; either version 3 of the License, or (at your option) any later
## version.
##
## This program is distributed in the hope that it will be useful, but WITHOUT
## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
## details.
##
## You should have received a copy of the GNU General Public License along with
## this program; if not, see <http://www.gnu.org/licenses/>.

## -*- texinfo -*-
## @deftypefn  {Function File} {@var{Q} = } marcumq (@var{a}, @var{b})
## @deftypefnx {Function File} {@var{Q} = } marcumq (@var{a}, @var{b}, @var{m})
## @deftypefnx {Function File} {@var{Q} = } marcumq (@var{a}, @var{b}, @var{m}, @var{tol})
##
## Compute the generalized Marcum Q function of order @var{M} with
## noncentrality parameter @var{a} and argument @var{b}.  If the order
## @var{M} is omitted it defaults to 1.  An optional relative tolerance
## @var{tol} may be included, the default is @code{eps}.
##
## If the input arguments are commensurate vectors, this function
## will produce a table of values.
##
## This function computes Marcum's Q function using the infinite
## Bessel series, truncated when the relative error is less than
## the specified tolerance.  The accuracy is limited by that of
## the Bessel functions, so reducing the tolerance is probably 
## not useful.
##
## Reference: Marcum, "Tables of Q Functions", Rand Corporation.
##
## Reference: R.T. Short, "Computation of Noncentral Chi-squared
## and Rice Random Variables", www.phaselockedsystems.com/publications
##
## @end deftypefn

function Q = marcumq(a,b,M=1,tol=eps)

  if ( (nargin<2) || (nargin>4) )
    print_usage();
  end
  if ( any(a<0) || any(b<0) )
    error("Parameters to marcumq must be positive");
  end
  if ( any(M<0) || any(floor(M)~=M) )
    error("M must be a positive integer");
  end

  [a,b,M] = tablify(a,b,M);

  Q = arrayfun(@mq, a,b,M,tol);

end

% Subfunction to compute the actual Marcum Q function.
function Q = mq(a,b,M,tol)
  % Special cases.
  if (b==0)
    Q = 1;
    N = 0;
    return;
  end
  if (a==0)
    k = 0:(M-1);
    Q = exp(-b^2/2)*sum(b.^(2*k)./(2.^k .* factorial(k)));
    N = 0;
    return;
  end

  % The basic iteration.  If a<b compute Q_M, otherwise
  % compute 1-Q_M.
  k = M;
  z = a*b;
  t = 1;
  k = 0;
  if (a<b)
    s = +1;
    c =  0;
    x = a/b;
    d = x;
    S = besseli(0,z,1);
    for k=1:M-1
      t = (d+1/d)*besseli(k,z,1);
      S = S + t;
      d = d*x;
    end
    N=k++;
  else
    s = -1;
    c =  1;
    x = b/a;
    k = M;
    d = x^M;
    S = 0;
    N = 0;
  end

  do
    t = d*besseli(abs(k),z,1);
    S = S + t;
    d = d*x;
    N = k++;
  until (abs(t/S)<tol)
  Q = c + s*exp( -(a-b)^2/2 )*S;
end

%  Tests for number and validity of arguments.
%
%!error
%! fail(marcumq(1))
%!error
%! fail(marcumq(-1,1,1,1,1))
%!error
%! fail(marcumq(-1,1))
%!error
%! fail(marcumq(1,-1))
%!error
%! fail(marcumq(1,1,-1))
%!error
%! fail(marcumq(1,1,1.1))

%  Notes on tests and accuracy.
%  -----------------------------------
%  The numbers used as the reference (Q) in the tables below are
%  from J.I. Marcum, "Table of Q Functions", Rand Technical Report
%  RM-339, 1950/1/1.
%
%  There is one discrepency in the tables.  Marcum has
%    Q(14.00,17.10) = 0.001078
%  while we compute 
%    Q(14.00,17.10) = 0.0010785053 = 0.001079
%  This is obviously a non-problem.
%
%  As further tests, I created several different versions of the
%  Q function computation, including a Bessel series expansion and
%  numerical integration.  All of them agree to with 10^(-16).

%!test
%! a = [0.00;0.05;1.00;2.00;3.00;4.00;5.00;6.00;7.00;8.00;9.00;10.00;
%!      11.00;12.00;13.00;14.00;15.00;16.00;17.00;18.00;19.00;20.00;
%!      21.00;22.00;23.00;24.00];
%! b = [0.000000, 0.100000, 1.100000, 2.100000, 3.100000, 4.100000]; 
%! Q = [1.000000, 0.995012, 0.546074, 0.110251, 0.008189, 0.000224;
%!      1.000000, 0.995019, 0.546487, 0.110554, 0.008238, 0.000226;
%!      1.000000, 0.996971, 0.685377, 0.233113, 0.034727, 0.002092;
%!      1.000000, 0.999322, 0.898073, 0.561704, 0.185328, 0.027068;
%!      1.000000, 0.999944, 0.985457, 0.865241, 0.526735, 0.169515;
%!      1.000000, 0.999998, 0.999136, 0.980933, 0.851679, 0.509876;
%!      1.000000, 1.000000, 0.999979, 0.998864, 0.978683, 0.844038;
%!      1.000000, 1.000000, 1.000000, 0.999973, 0.998715, 0.977300;
%!      1.000000, 1.000000, 1.000000, 1.000000, 0.999969, 0.998618;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.999966;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000]; 
%! q = marcumq(a,b);
%! assert(q,Q,1e-6);

%!test
%! a = [0.00;0.05;1.00;2.00;3.00;4.00;5.00;6.00;7.00;8.00;9.00;10.00;
%!      11.00;12.00;13.00;14.00;15.00;16.00;17.00;18.00;19.00;20.00;
%!      21.00;22.00;23.00;24.00];
%! b = [5.100000, 6.100000, 7.100000, 8.100000, 9.100000, 10.10000]; 
%! Q = [0.000002, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000;
%!      0.000002, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000;
%!      0.000049, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000;
%!      0.001606, 0.000037, 0.000000, 0.000000, 0.000000, 0.000000;
%!      0.024285, 0.001420, 0.000032, 0.000000, 0.000000, 0.000000;
%!      0.161412, 0.022812, 0.001319, 0.000030, 0.000000, 0.000000;
%!      0.499869, 0.156458, 0.021893, 0.001256, 0.000028, 0.000000;
%!      0.839108, 0.493229, 0.153110, 0.021264, 0.001212, 0.000027;
%!      0.976358, 0.835657, 0.488497, 0.150693, 0.020806, 0.001180;
%!      0.998549, 0.975673, 0.833104, 0.484953, 0.148867, 0.020458;
%!      0.999965, 0.998498, 0.975152, 0.831138, 0.482198, 0.147437;
%!      1.000000, 0.999963, 0.998458, 0.974742, 0.829576, 0.479995;
%!      1.000000, 1.000000, 0.999962, 0.998426, 0.974411, 0.828307;
%!      1.000000, 1.000000, 1.000000, 0.999961, 0.998400, 0.974138;
%!      1.000000, 1.000000, 1.000000, 1.000000, 0.999960, 0.998378;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.999960;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000]; 
%! q = marcumq(a,b);
%! assert(q,Q,1e-6);


%!test
%! a = [0.00;0.05;1.00;2.00;3.00;4.00;5.00;6.00;7.00;8.00;9.00;10.00;
%!      11.00;12.00;13.00;14.00;15.00;16.00;17.00;18.00;19.00;20.00;
%!      21.00;22.00;23.00;24.00];
%! b = [11.10000, 12.10000, 13.10000, 14.10000, 15.10000, 16.10000]; 
%! Q = [0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000;
%!      0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000;
%!      0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000;
%!      0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000;
%!      0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000;
%!      0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000;
%!      0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000;
%!      0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000;
%!      0.000026, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000;
%!      0.001155, 0.000026, 0.000000, 0.000000, 0.000000, 0.000000;
%!      0.020183, 0.001136, 0.000025, 0.000000, 0.000000, 0.000000;
%!      0.146287, 0.019961, 0.001120, 0.000025, 0.000000, 0.000000;
%!      0.478193, 0.145342, 0.019778, 0.001107, 0.000024, 0.000000;
%!      0.827253, 0.476692, 0.144551, 0.019625, 0.001096, 0.000024;
%!      0.973909, 0.826366, 0.475422, 0.143881, 0.019494, 0.001087;
%!      0.998359, 0.973714, 0.825607, 0.474333, 0.143304, 0.019381;
%!      0.999959, 0.998343, 0.973546, 0.824952, 0.473389, 0.142803;
%!      1.000000, 0.999959, 0.998330, 0.973400, 0.824380, 0.472564;
%!      1.000000, 1.000000, 0.999958, 0.998318, 0.973271, 0.823876;
%!      1.000000, 1.000000, 1.000000, 0.999958, 0.998307, 0.973158;
%!      1.000000, 1.000000, 1.000000, 1.000000, 0.999957, 0.998297;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.999957;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000]; 
%! q = marcumq(a,b);
%! assert(q,Q,1e-6);


%!test
%! a = [0.00;0.05;1.00;2.00;3.00;4.00;5.00;6.00;7.00;8.00;9.00;10.00;
%!      11.00;12.00;13.00;14.00;15.00;16.00;17.00;18.00;19.00;20.00;
%!      21.00;22.00;23.00;24.00];
%! b = [17.10000, 18.10000, 19.10000]; 
%! Q = [0.000000, 0.000000, 0.000000;
%!      0.000000, 0.000000, 0.000000;
%!      0.000000, 0.000000, 0.000000;
%!      0.000000, 0.000000, 0.000000;
%!      0.000000, 0.000000, 0.000000;
%!      0.000000, 0.000000, 0.000000;
%!      0.000000, 0.000000, 0.000000;
%!      0.000000, 0.000000, 0.000000;
%!      0.000000, 0.000000, 0.000000;
%!      0.000000, 0.000000, 0.000000;
%!      0.000000, 0.000000, 0.000000;
%!      0.000000, 0.000000, 0.000000;
%!      0.000000, 0.000000, 0.000000;
%!      0.000000, 0.000000, 0.000000;
%!      0.000024, 0.000000, 0.000000;
%!      0.001078, 0.000024, 0.000000;
%!      0.019283, 0.001071, 0.000023;
%!      0.142364, 0.019197, 0.001065;
%!      0.471835, 0.141976, 0.019121;
%!      0.823429, 0.471188, 0.141630;
%!      0.973056, 0.823030, 0.470608;
%!      0.998289, 0.972965, 0.822671;
%!      0.999957, 0.998281, 0.972883;
%!      1.000000, 0.999957, 0.998274;
%!      1.000000, 1.000000, 0.999956;
%!      1.000000, 1.000000, 1.000000]; 
%! q = marcumq(a,b);
%! assert(q,Q,1e-6);

%  The tests for M>1 were generating from Marcum's tables by
%  using the formula
%    Q_M(a,b) = Q(a,b) + exp(-(a-b)^2/2)*sum_{k=1}^{M-1}(b/a)^k*exp(-ab)*I_k(ab)
%
%!test
%! M = 2;
%! a = [0.00;0.05;1.00;2.00;3.00;4.00;5.00;6.00;7.00;8.00;9.00;10.00;11.00;12.00;13.00;14.00;15.00;16.00;17.00;18.00;19.00;20.00;21.00;22.00;23.00;24.000000];
%! 
%! b = [    0.00,     0.10,     2.10,     7.10,    12.10,    17.10]; 
%! Q = [1.000000, 0.999987, 0.353353, 0.000000, 0.000000, 0.000000;
%!      1.000000, 0.999988, 0.353687, 0.000000, 0.000000, 0.000000;
%!      1.000000, 0.999992, 0.478229, 0.000000, 0.000000, 0.000000;
%!      1.000000, 0.999999, 0.745094, 0.000001, 0.000000, 0.000000;
%!      1.000000, 1.000000, 0.934771, 0.000077, 0.000000, 0.000000;
%!      1.000000, 1.000000, 0.992266, 0.002393, 0.000000, 0.000000;
%!      1.000000, 1.000000, 0.999607, 0.032264, 0.000000, 0.000000;
%!      1.000000, 1.000000, 0.999992, 0.192257, 0.000000, 0.000000;
%!      1.000000, 1.000000, 1.000000, 0.545174, 0.000000, 0.000000;
%!      1.000000, 1.000000, 1.000000, 0.864230, 0.000040, 0.000000;
%!      1.000000, 1.000000, 1.000000, 0.981589, 0.001555, 0.000000;
%!      1.000000, 1.000000, 1.000000, 0.998957, 0.024784, 0.000000;
%!      1.000000, 1.000000, 1.000000, 0.999976, 0.166055, 0.000000;
%!      1.000000, 1.000000, 1.000000, 1.000000, 0.509823, 0.000000;
%!      1.000000, 1.000000, 1.000000, 1.000000, 0.846066, 0.000032;
%!      1.000000, 1.000000, 1.000000, 1.000000, 0.978062, 0.001335;
%!      1.000000, 1.000000, 1.000000, 1.000000, 0.998699, 0.022409;
%!      1.000000, 1.000000, 1.000000, 1.000000, 0.999970, 0.156421;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.495223;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.837820;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.976328;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.998564;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.999966;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000]; 
%! q = marcumq(a,b,M);
%! assert(q,Q,1e-6);

%!test
%! M = 5;
%! a = [0.00;0.05;1.00;2.00;3.00;4.00;5.00;6.00;7.00;8.00;9.00;10.00;11.00;12.00;13.00;14.00;15.00;16.00;17.00;18.00;19.00;20.00;21.00;22.00;23.00;24.000000];
%! 
%! b = [    0.00,     0.10,     2.10,     7.10,    12.10,    17.10]; 
%! Q = [1.000000, 1.000000, 0.926962, 0.000000, 0.000000, 0.000000;
%!      1.000000, 1.000000, 0.927021, 0.000000, 0.000000, 0.000000;
%!      1.000000, 1.000000, 0.947475, 0.000001, 0.000000, 0.000000;
%!      1.000000, 1.000000, 0.980857, 0.000033, 0.000000, 0.000000;
%!      1.000000, 1.000000, 0.996633, 0.000800, 0.000000, 0.000000;
%!      1.000000, 1.000000, 0.999729, 0.011720, 0.000000, 0.000000;
%!      1.000000, 1.000000, 0.999990, 0.088999, 0.000000, 0.000000;
%!      1.000000, 1.000000, 1.000000, 0.341096, 0.000000, 0.000000;
%!      1.000000, 1.000000, 1.000000, 0.705475, 0.000002, 0.000000;
%!      1.000000, 1.000000, 1.000000, 0.933009, 0.000134, 0.000000;
%!      1.000000, 1.000000, 1.000000, 0.993118, 0.003793, 0.000000;
%!      1.000000, 1.000000, 1.000000, 0.999702, 0.045408, 0.000000;
%!      1.000000, 1.000000, 1.000000, 0.999995, 0.238953, 0.000000;
%!      1.000000, 1.000000, 1.000000, 1.000000, 0.607903, 0.000001;
%!      1.000000, 1.000000, 1.000000, 1.000000, 0.896007, 0.000073;
%!      1.000000, 1.000000, 1.000000, 1.000000, 0.987642, 0.002480;
%!      1.000000, 1.000000, 1.000000, 1.000000, 0.999389, 0.034450;
%!      1.000000, 1.000000, 1.000000, 1.000000, 0.999988, 0.203879;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.565165;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.876284;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.984209;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.999165;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.999983;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000]; 
%! q = marcumq(a,b,M);
%! assert(q,Q,1e-6);

%!test
%! M = 10;
%! a = [0.00;0.05;1.00;2.00;3.00;4.00;5.00;6.00;7.00;8.00;9.00;10.00;11.00;12.00;13.00;14.00;15.00;16.00;17.00;18.00;19.00;20.00;21.00;22.00;23.00;24.000000];
%! 
%! b = [    0.00,     0.10,     2.10,     7.10,    12.10,    17.10]; 
%! Q = [1.000000, 1.000000, 0.999898, 0.000193, 0.000000, 0.000000;
%!      1.000000, 1.000000, 0.999897, 0.000194, 0.000000, 0.000000;
%!      1.000000, 1.000000, 0.999931, 0.000416, 0.000000, 0.000000;
%!      1.000000, 1.000000, 0.999980, 0.002377, 0.000000, 0.000000;
%!      1.000000, 1.000000, 0.999997, 0.016409, 0.000000, 0.000000;
%!      1.000000, 1.000000, 0.999999, 0.088005, 0.000000, 0.000000;
%!      1.000000, 1.000000, 1.000000, 0.302521, 0.000000, 0.000000;
%!      1.000000, 1.000000, 1.000000, 0.638401, 0.000000, 0.000000;
%!      1.000000, 1.000000, 1.000000, 0.894322, 0.000022, 0.000000;
%!      1.000000, 1.000000, 1.000000, 0.984732, 0.000840, 0.000000;
%!      1.000000, 1.000000, 1.000000, 0.998997, 0.014160, 0.000000;
%!      1.000000, 1.000000, 1.000000, 0.999972, 0.107999, 0.000000;
%!      1.000000, 1.000000, 1.000000, 1.000000, 0.391181, 0.000000;
%!      1.000000, 1.000000, 1.000000, 1.000000, 0.754631, 0.000004;
%!      1.000000, 1.000000, 1.000000, 1.000000, 0.951354, 0.000266;
%!      1.000000, 1.000000, 1.000000, 1.000000, 0.995732, 0.006444;
%!      1.000000, 1.000000, 1.000000, 1.000000, 0.999843, 0.065902;
%!      1.000000, 1.000000, 1.000000, 1.000000, 0.999998, 0.299616;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.676336;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.925312;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.992390;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.999679;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.999995;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
%!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000]; 
%! q = marcumq(a,b,M);
%! assert(q,Q,1e-6);