X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?p=CreaPhase.git;a=blobdiff_plain;f=octave_packages%2Fsignal-1.1.3%2Fmarcumq.m;fp=octave_packages%2Fsignal-1.1.3%2Fmarcumq.m;h=d9dac9e823b2f988868386cca1d9216fef29b3bb;hp=0000000000000000000000000000000000000000;hb=c880e8788dfc484bf23ce13fa2787f2c6bca4863;hpb=1705066eceaaea976f010f669ce8e972f3734b05

diff --git a/octave_packages/signal-1.1.3/marcumq.m b/octave_packages/signal-1.1.3/marcumq.m
new file mode 100644
index 0000000..d9dac9e
--- /dev/null
+++ b/octave_packages/signal-1.1.3/marcumq.m
@@ -0,0 +1,389 @@
+## 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
+
+  nr = max([size(a,1) size(b,1) size(M,1)]);
+  nc = max([size(a,2) size(b,2) size(M,2)]);
+  a = padarray(a, [nr - size(a,1) nc - size(a,2)], "replicate", "post");
+  b = padarray(b, [nr - size(b,1) nc - size(b,2)], "replicate", "post");
+  M = padarray(M, [nr - size(M,1) nc - size(M,2)], "replicate", "post");
+
+  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);