octave_packages/signal-1.1.3/marcumq.m

   1 ## Copyright (C) 2012   Robert T. Short     <rtshort@ieee.org>
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify it under
   4 ## the terms of the GNU General Public License as published by the Free Software
   5 ## Foundation; either version 3 of the License, or (at your option) any later
   6 ## version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful, but WITHOUT
   9 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  10 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  11 ## details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License along with
  14 ## this program; if not, see <http://www.gnu.org/licenses/>.
  15
  16 ## -*- texinfo -*-
  17 ## @deftypefn  {Function File} {@var{Q} = } marcumq (@var{a}, @var{b})
  18 ## @deftypefnx {Function File} {@var{Q} = } marcumq (@var{a}, @var{b}, @var{m})
  19 ## @deftypefnx {Function File} {@var{Q} = } marcumq (@var{a}, @var{b}, @var{m}, @var{tol})
  20 ##
  21 ## Compute the generalized Marcum Q function of order @var{M} with
  22 ## noncentrality parameter @var{a} and argument @var{b}.  If the order
  23 ## @var{M} is omitted it defaults to 1.  An optional relative tolerance
  24 ## @var{tol} may be included, the default is @code{eps}.
  25 ##
  26 ## If the input arguments are commensurate vectors, this function
  27 ## will produce a table of values.
  28 ##
  29 ## This function computes Marcum's Q function using the infinite
  30 ## Bessel series, truncated when the relative error is less than
  31 ## the specified tolerance.  The accuracy is limited by that of
  32 ## the Bessel functions, so reducing the tolerance is probably
  33 ## not useful.
  34 ##
  35 ## Reference: Marcum, "Tables of Q Functions", Rand Corporation.
  36 ##
  37 ## Reference: R.T. Short, "Computation of Noncentral Chi-squared
  38 ## and Rice Random Variables", www.phaselockedsystems.com/publications
  39 ##
  40 ## @end deftypefn
  41
  42 function Q = marcumq(a,b,M=1,tol=eps)
  43
  44   if ( (nargin<2) || (nargin>4) )
  45     print_usage();
  46   end
  47   if ( any(a<0) || any(b<0) )
  48     error("Parameters to marcumq must be positive");
  49   end
  50   if ( any(M<0) || any(floor(M)~=M) )
  51     error("M must be a positive integer");
  52   end
  53
  54   nr = max([size(a,1) size(b,1) size(M,1)]);
  55   nc = max([size(a,2) size(b,2) size(M,2)]);
  56   a = padarray(a, [nr - size(a,1) nc - size(a,2)], "replicate", "post");
  57   b = padarray(b, [nr - size(b,1) nc - size(b,2)], "replicate", "post");
  58   M = padarray(M, [nr - size(M,1) nc - size(M,2)], "replicate", "post");
  59
  60   Q = arrayfun(@mq, a,b,M,tol);
  61
  62 end
  63
  64 % Subfunction to compute the actual Marcum Q function.
  65 function Q = mq(a,b,M,tol)
  66   % Special cases.
  67   if (b==0)
  68     Q = 1;
  69     N = 0;
  70     return;
  71   end
  72   if (a==0)
  73     k = 0:(M-1);
  74     Q = exp(-b^2/2)*sum(b.^(2*k)./(2.^k .* factorial(k)));
  75     N = 0;
  76     return;
  77   end
  78
  79   % The basic iteration.  If a<b compute Q_M, otherwise
  80   % compute 1-Q_M.
  81   k = M;
  82   z = a*b;
  83   t = 1;
  84   k = 0;
  85   if (a<b)
  86     s = +1;
  87     c =  0;
  88     x = a/b;
  89     d = x;
  90     S = besseli(0,z,1);
  91     for k=1:M-1
  92       t = (d+1/d)*besseli(k,z,1);
  93       S = S + t;
  94       d = d*x;
  95     end
  96     N=k++;
  97   else
  98     s = -1;
  99     c =  1;
 100     x = b/a;
 101     k = M;
 102     d = x^M;
 103     S = 0;
 104     N = 0;
 105   end
 106
 107   do
 108     t = d*besseli(abs(k),z,1);
 109     S = S + t;
 110     d = d*x;
 111     N = k++;
 112   until (abs(t/S)<tol)
 113   Q = c + s*exp( -(a-b)^2/2 )*S;
 114 end
 115
 116 %  Tests for number and validity of arguments.
 117 %
 118 %!error
 119 %! fail(marcumq(1))
 120 %!error
 121 %! fail(marcumq(-1,1,1,1,1))
 122 %!error
 123 %! fail(marcumq(-1,1))
 124 %!error
 125 %! fail(marcumq(1,-1))
 126 %!error
 127 %! fail(marcumq(1,1,-1))
 128 %!error
 129 %! fail(marcumq(1,1,1.1))
 130
 131 %  Notes on tests and accuracy.
 132 %  -----------------------------------
 133 %  The numbers used as the reference (Q) in the tables below are
 134 %  from J.I. Marcum, "Table of Q Functions", Rand Technical Report
 135 %  RM-339, 1950/1/1.
 136 %
 137 %  There is one discrepency in the tables.  Marcum has
 138 %    Q(14.00,17.10) = 0.001078
 139 %  while we compute
 140 %    Q(14.00,17.10) = 0.0010785053 = 0.001079
 141 %  This is obviously a non-problem.
 142 %
 143 %  As further tests, I created several different versions of the
 144 %  Q function computation, including a Bessel series expansion and
 145 %  numerical integration.  All of them agree to with 10^(-16).
 146
 147 %!test
 148 %! 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;
 149 %!      11.00;12.00;13.00;14.00;15.00;16.00;17.00;18.00;19.00;20.00;
 150 %!      21.00;22.00;23.00;24.00];
 151 %! b = [0.000000, 0.100000, 1.100000, 2.100000, 3.100000, 4.100000];
 152 %! Q = [1.000000, 0.995012, 0.546074, 0.110251, 0.008189, 0.000224;
 153 %!      1.000000, 0.995019, 0.546487, 0.110554, 0.008238, 0.000226;
 154 %!      1.000000, 0.996971, 0.685377, 0.233113, 0.034727, 0.002092;
 155 %!      1.000000, 0.999322, 0.898073, 0.561704, 0.185328, 0.027068;
 156 %!      1.000000, 0.999944, 0.985457, 0.865241, 0.526735, 0.169515;
 157 %!      1.000000, 0.999998, 0.999136, 0.980933, 0.851679, 0.509876;
 158 %!      1.000000, 1.000000, 0.999979, 0.998864, 0.978683, 0.844038;
 159 %!      1.000000, 1.000000, 1.000000, 0.999973, 0.998715, 0.977300;
 160 %!      1.000000, 1.000000, 1.000000, 1.000000, 0.999969, 0.998618;
 161 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.999966;
 162 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
 163 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
 164 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
 165 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
 166 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
 167 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
 168 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
 169 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
 170 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
 171 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
 172 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
 173 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
 174 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
 175 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
 176 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
 177 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000];
 178 %! q = marcumq(a,b);
 179 %! assert(q,Q,1e-6);
 180
 181 %!test
 182 %! 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;
 183 %!      11.00;12.00;13.00;14.00;15.00;16.00;17.00;18.00;19.00;20.00;
 184 %!      21.00;22.00;23.00;24.00];
 185 %! b = [5.100000, 6.100000, 7.100000, 8.100000, 9.100000, 10.10000];
 186 %! Q = [0.000002, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000;
 187 %!      0.000002, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000;
 188 %!      0.000049, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000;
 189 %!      0.001606, 0.000037, 0.000000, 0.000000, 0.000000, 0.000000;
 190 %!      0.024285, 0.001420, 0.000032, 0.000000, 0.000000, 0.000000;
 191 %!      0.161412, 0.022812, 0.001319, 0.000030, 0.000000, 0.000000;
 192 %!      0.499869, 0.156458, 0.021893, 0.001256, 0.000028, 0.000000;
 193 %!      0.839108, 0.493229, 0.153110, 0.021264, 0.001212, 0.000027;
 194 %!      0.976358, 0.835657, 0.488497, 0.150693, 0.020806, 0.001180;
 195 %!      0.998549, 0.975673, 0.833104, 0.484953, 0.148867, 0.020458;
 196 %!      0.999965, 0.998498, 0.975152, 0.831138, 0.482198, 0.147437;
 197 %!      1.000000, 0.999963, 0.998458, 0.974742, 0.829576, 0.479995;
 198 %!      1.000000, 1.000000, 0.999962, 0.998426, 0.974411, 0.828307;
 199 %!      1.000000, 1.000000, 1.000000, 0.999961, 0.998400, 0.974138;
 200 %!      1.000000, 1.000000, 1.000000, 1.000000, 0.999960, 0.998378;
 201 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.999960;
 202 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
 203 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
 204 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
 205 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
 206 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
 207 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
 208 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
 209 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
 210 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
 211 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000];
 212 %! q = marcumq(a,b);
 213 %! assert(q,Q,1e-6);
 214
 215
 216 %!test
 217 %! 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;
 218 %!      11.00;12.00;13.00;14.00;15.00;16.00;17.00;18.00;19.00;20.00;
 219 %!      21.00;22.00;23.00;24.00];
 220 %! b = [11.10000, 12.10000, 13.10000, 14.10000, 15.10000, 16.10000];
 221 %! Q = [0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000;
 222 %!      0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000;
 223 %!      0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000;
 224 %!      0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000;
 225 %!      0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000;
 226 %!      0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000;
 227 %!      0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000;
 228 %!      0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000;
 229 %!      0.000026, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000;
 230 %!      0.001155, 0.000026, 0.000000, 0.000000, 0.000000, 0.000000;
 231 %!      0.020183, 0.001136, 0.000025, 0.000000, 0.000000, 0.000000;
 232 %!      0.146287, 0.019961, 0.001120, 0.000025, 0.000000, 0.000000;
 233 %!      0.478193, 0.145342, 0.019778, 0.001107, 0.000024, 0.000000;
 234 %!      0.827253, 0.476692, 0.144551, 0.019625, 0.001096, 0.000024;
 235 %!      0.973909, 0.826366, 0.475422, 0.143881, 0.019494, 0.001087;
 236 %!      0.998359, 0.973714, 0.825607, 0.474333, 0.143304, 0.019381;
 237 %!      0.999959, 0.998343, 0.973546, 0.824952, 0.473389, 0.142803;
 238 %!      1.000000, 0.999959, 0.998330, 0.973400, 0.824380, 0.472564;
 239 %!      1.000000, 1.000000, 0.999958, 0.998318, 0.973271, 0.823876;
 240 %!      1.000000, 1.000000, 1.000000, 0.999958, 0.998307, 0.973158;
 241 %!      1.000000, 1.000000, 1.000000, 1.000000, 0.999957, 0.998297;
 242 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.999957;
 243 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
 244 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
 245 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
 246 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000];
 247 %! q = marcumq(a,b);
 248 %! assert(q,Q,1e-6);
 249
 250
 251 %!test
 252 %! 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;
 253 %!      11.00;12.00;13.00;14.00;15.00;16.00;17.00;18.00;19.00;20.00;
 254 %!      21.00;22.00;23.00;24.00];
 255 %! b = [17.10000, 18.10000, 19.10000];
 256 %! Q = [0.000000, 0.000000, 0.000000;
 257 %!      0.000000, 0.000000, 0.000000;
 258 %!      0.000000, 0.000000, 0.000000;
 259 %!      0.000000, 0.000000, 0.000000;
 260 %!      0.000000, 0.000000, 0.000000;
 261 %!      0.000000, 0.000000, 0.000000;
 262 %!      0.000000, 0.000000, 0.000000;
 263 %!      0.000000, 0.000000, 0.000000;
 264 %!      0.000000, 0.000000, 0.000000;
 265 %!      0.000000, 0.000000, 0.000000;
 266 %!      0.000000, 0.000000, 0.000000;
 267 %!      0.000000, 0.000000, 0.000000;
 268 %!      0.000000, 0.000000, 0.000000;
 269 %!      0.000000, 0.000000, 0.000000;
 270 %!      0.000024, 0.000000, 0.000000;
 271 %!      0.001078, 0.000024, 0.000000;
 272 %!      0.019283, 0.001071, 0.000023;
 273 %!      0.142364, 0.019197, 0.001065;
 274 %!      0.471835, 0.141976, 0.019121;
 275 %!      0.823429, 0.471188, 0.141630;
 276 %!      0.973056, 0.823030, 0.470608;
 277 %!      0.998289, 0.972965, 0.822671;
 278 %!      0.999957, 0.998281, 0.972883;
 279 %!      1.000000, 0.999957, 0.998274;
 280 %!      1.000000, 1.000000, 0.999956;
 281 %!      1.000000, 1.000000, 1.000000];
 282 %! q = marcumq(a,b);
 283 %! assert(q,Q,1e-6);
 284
 285 %  The tests for M>1 were generating from Marcum's tables by
 286 %  using the formula
 287 %    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)
 288 %
 289 %!test
 290 %! M = 2;
 291 %! 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];
 292 %!
 293 %! b = [    0.00,     0.10,     2.10,     7.10,    12.10,    17.10];
 294 %! Q = [1.000000, 0.999987, 0.353353, 0.000000, 0.000000, 0.000000;
 295 %!      1.000000, 0.999988, 0.353687, 0.000000, 0.000000, 0.000000;
 296 %!      1.000000, 0.999992, 0.478229, 0.000000, 0.000000, 0.000000;
 297 %!      1.000000, 0.999999, 0.745094, 0.000001, 0.000000, 0.000000;
 298 %!      1.000000, 1.000000, 0.934771, 0.000077, 0.000000, 0.000000;
 299 %!      1.000000, 1.000000, 0.992266, 0.002393, 0.000000, 0.000000;
 300 %!      1.000000, 1.000000, 0.999607, 0.032264, 0.000000, 0.000000;
 301 %!      1.000000, 1.000000, 0.999992, 0.192257, 0.000000, 0.000000;
 302 %!      1.000000, 1.000000, 1.000000, 0.545174, 0.000000, 0.000000;
 303 %!      1.000000, 1.000000, 1.000000, 0.864230, 0.000040, 0.000000;
 304 %!      1.000000, 1.000000, 1.000000, 0.981589, 0.001555, 0.000000;
 305 %!      1.000000, 1.000000, 1.000000, 0.998957, 0.024784, 0.000000;
 306 %!      1.000000, 1.000000, 1.000000, 0.999976, 0.166055, 0.000000;
 307 %!      1.000000, 1.000000, 1.000000, 1.000000, 0.509823, 0.000000;
 308 %!      1.000000, 1.000000, 1.000000, 1.000000, 0.846066, 0.000032;
 309 %!      1.000000, 1.000000, 1.000000, 1.000000, 0.978062, 0.001335;
 310 %!      1.000000, 1.000000, 1.000000, 1.000000, 0.998699, 0.022409;
 311 %!      1.000000, 1.000000, 1.000000, 1.000000, 0.999970, 0.156421;
 312 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.495223;
 313 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.837820;
 314 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.976328;
 315 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.998564;
 316 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.999966;
 317 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
 318 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
 319 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000];
 320 %! q = marcumq(a,b,M);
 321 %! assert(q,Q,1e-6);
 322
 323 %!test
 324 %! M = 5;
 325 %! 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];
 326 %!
 327 %! b = [    0.00,     0.10,     2.10,     7.10,    12.10,    17.10];
 328 %! Q = [1.000000, 1.000000, 0.926962, 0.000000, 0.000000, 0.000000;
 329 %!      1.000000, 1.000000, 0.927021, 0.000000, 0.000000, 0.000000;
 330 %!      1.000000, 1.000000, 0.947475, 0.000001, 0.000000, 0.000000;
 331 %!      1.000000, 1.000000, 0.980857, 0.000033, 0.000000, 0.000000;
 332 %!      1.000000, 1.000000, 0.996633, 0.000800, 0.000000, 0.000000;
 333 %!      1.000000, 1.000000, 0.999729, 0.011720, 0.000000, 0.000000;
 334 %!      1.000000, 1.000000, 0.999990, 0.088999, 0.000000, 0.000000;
 335 %!      1.000000, 1.000000, 1.000000, 0.341096, 0.000000, 0.000000;
 336 %!      1.000000, 1.000000, 1.000000, 0.705475, 0.000002, 0.000000;
 337 %!      1.000000, 1.000000, 1.000000, 0.933009, 0.000134, 0.000000;
 338 %!      1.000000, 1.000000, 1.000000, 0.993118, 0.003793, 0.000000;
 339 %!      1.000000, 1.000000, 1.000000, 0.999702, 0.045408, 0.000000;
 340 %!      1.000000, 1.000000, 1.000000, 0.999995, 0.238953, 0.000000;
 341 %!      1.000000, 1.000000, 1.000000, 1.000000, 0.607903, 0.000001;
 342 %!      1.000000, 1.000000, 1.000000, 1.000000, 0.896007, 0.000073;
 343 %!      1.000000, 1.000000, 1.000000, 1.000000, 0.987642, 0.002480;
 344 %!      1.000000, 1.000000, 1.000000, 1.000000, 0.999389, 0.034450;
 345 %!      1.000000, 1.000000, 1.000000, 1.000000, 0.999988, 0.203879;
 346 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.565165;
 347 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.876284;
 348 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.984209;
 349 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.999165;
 350 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.999983;
 351 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
 352 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
 353 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000];
 354 %! q = marcumq(a,b,M);
 355 %! assert(q,Q,1e-6);
 356
 357 %!test
 358 %! M = 10;
 359 %! 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];
 360 %!
 361 %! b = [    0.00,     0.10,     2.10,     7.10,    12.10,    17.10];
 362 %! Q = [1.000000, 1.000000, 0.999898, 0.000193, 0.000000, 0.000000;
 363 %!      1.000000, 1.000000, 0.999897, 0.000194, 0.000000, 0.000000;
 364 %!      1.000000, 1.000000, 0.999931, 0.000416, 0.000000, 0.000000;
 365 %!      1.000000, 1.000000, 0.999980, 0.002377, 0.000000, 0.000000;
 366 %!      1.000000, 1.000000, 0.999997, 0.016409, 0.000000, 0.000000;
 367 %!      1.000000, 1.000000, 0.999999, 0.088005, 0.000000, 0.000000;
 368 %!      1.000000, 1.000000, 1.000000, 0.302521, 0.000000, 0.000000;
 369 %!      1.000000, 1.000000, 1.000000, 0.638401, 0.000000, 0.000000;
 370 %!      1.000000, 1.000000, 1.000000, 0.894322, 0.000022, 0.000000;
 371 %!      1.000000, 1.000000, 1.000000, 0.984732, 0.000840, 0.000000;
 372 %!      1.000000, 1.000000, 1.000000, 0.998997, 0.014160, 0.000000;
 373 %!      1.000000, 1.000000, 1.000000, 0.999972, 0.107999, 0.000000;
 374 %!      1.000000, 1.000000, 1.000000, 1.000000, 0.391181, 0.000000;
 375 %!      1.000000, 1.000000, 1.000000, 1.000000, 0.754631, 0.000004;
 376 %!      1.000000, 1.000000, 1.000000, 1.000000, 0.951354, 0.000266;
 377 %!      1.000000, 1.000000, 1.000000, 1.000000, 0.995732, 0.006444;
 378 %!      1.000000, 1.000000, 1.000000, 1.000000, 0.999843, 0.065902;
 379 %!      1.000000, 1.000000, 1.000000, 1.000000, 0.999998, 0.299616;
 380 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.676336;
 381 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.925312;
 382 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.992390;
 383 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.999679;
 384 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.999995;
 385 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
 386 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000;
 387 %!      1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000];
 388 %! q = marcumq(a,b,M);
 389 %! assert(q,Q,1e-6);