summaryrefslogtreecommitdiff
path: root/macros/marcumq.sci
diff options
context:
space:
mode:
Diffstat (limited to 'macros/marcumq.sci')
-rw-r--r--macros/marcumq.sci423
1 files changed, 387 insertions, 36 deletions
diff --git a/macros/marcumq.sci b/macros/marcumq.sci
index 630ec97..29f4416 100644
--- a/macros/marcumq.sci
+++ b/macros/marcumq.sci
@@ -1,38 +1,389 @@
+// Copyright (C) 2018 - IIT Bombay - FOSSEE
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution. The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
+// Original Source : https://octave.sourceforge.io/
+// Modifieded by: Abinash Singh Under FOSSEE Internship
+// Last Modified on : 3 Feb 2024
+// Organization: FOSSEE, IIT Bombay
+// Email: toolbox@scilab.in
+/*
+Calling Sequence
+ q = marcumq (a, b)
+ q = marcumq (a, b, m)
+ q = marcumq (a, b, m, tol)
+
+Input and Output parameters
+a — Noncentrality parameter --- nonnegative scalar | array of nonnegative numbers
+b — Argument of Marcum Q-function --- nonnegative scalar | array of nonnegative numbers
+m — Order of generalized Marcum Q-function --- positive integer | array of positive integers
+
+Compute the generalized Marcum Q function of order `m` with noncentrality parameter `a` and argument `b`.
+If the order `m` is omitted, it defaults to 1. An optional relative tolerance `tol` may be included,
+and the default value is `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,
+which is 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.
+
+References:
+- Marcum, "Tables of Q Functions", Rand Corporation.
+- R.T. Short, "Computation of Noncentral Chi-squared and Rice Random Variables",
+ www.phaselockedsystems.com/publications
+*/
+
function q = marcumq (a, b, m, tol)
-//This function computes the generalized Marcum Q function of order m with noncentrality parameter a and argument b.
-//Calling Sequence
-//q = marcumq (a, b)
-//q = marcumq (a, b, m)
-//q = marcumq (a, b, m, tol)
-//Parameters
-//a:
-//b:
-//m: default value 1
-//tol: default value eps
-//Description
-//This is an Octave function.
-//This function computes the generalized Marcum Q function of order m with noncentrality parameter a and argument b.
-//The third argument m is the order, which by default is 1.
-//The fourth argument tol is the tolerance, which by default is eps.
-//If input arguments are vectors which correspond in size and degree, the output is a table of values.
-//This function calculates Marcum’s Q function using the infinite Bessel series, which is truncated when the relative error is less than the specified tolerance.
-//Examples
-//marcumq([1,2,3],4)
-//ans =
-// 0.0028895 0.0341348 0.1965122
-
-funcprot(0);
-rhs = argn(2)
-if(rhs<2 | rhs>4)
-error("Wrong number of input arguments.")
-end
- select(rhs)
- case 2 then
- q = callOctave("marcumq",a,b)
- case 3 then
- q = callOctave("marcumq",a,b,m)
- case 4 then
- q = callOctave("marcumq",a,b,m,tol)
- end
+
+ if ((nargin < 2) || (nargin > 5))
+ error(" marcumq : wrong numbers of input arguments ");
+ end
+
+ if nargin < 3 then m = 1 end
+ if nargin < 4 then tol = %eps end
+ if nargin < 5 then max_iter = 100 end
+
+ if (or (a < 0))
+ error ("marcumq: A must be a non-negative value");
+ end
+ if (or (b < 0))
+ error ("marcumq: B must be a non-negative value");
+ end
+ if (or (m < 1) || or (fix (m) ~= m))
+ error ("marcumq: M must be a positive integer");
+ end
+
+ [a, b] = tablify (a, b);
+ q = []
+ for i=1:size(a,1)
+ for j=1:size(a,2)
+ q(i,j) = mq(a(i,j),b(i,j),m,tol)
+ end
+ end
+
endfunction
-
+
+// 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);
+ if (m > 1)
+ for k = 1:m - 1
+ t = (d + 1 / d) * besseli (k, z, 1);
+ S = S + t;
+ d = d * x;
+ end
+ end
+ N = k;
+ k = k + 1
+ else
+ s = -1;
+ c = 1;
+ x = b / a;
+ k = m;
+ d = x^m;
+ S = 0;
+ N = 0;
+ end
+
+ while (abs (t / S) > tol)
+ t = d * besseli (abs (k), z, 1);
+ S = S + t;
+ d = d * x;
+ N = k;
+ k = k + 1 ;
+ end
+ q = c + s * exp (-(a - b)^2 / 2) * S;
+
+endfunction
+
+// Internal helper function to create a table of like dimensions from arguments.
+function [ta , tb] = tablify(a,b)
+ rows = size(a,1)
+ cols = size(b,2)
+ ta=[]
+ for i=1:cols
+ ta = [ta a ]
+ end
+ tb=[]
+ for i=1:rows
+ tb = [ tb ; b]
+ end
+endfunction
+/*
+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_checkalmostequal (q, Q, %eps,1e-4);
+
+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_checkalmostequal(q, Q, %eps,1e-4);
+
+
+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_checkalmostequal (q, Q,%eps,1e-4);
+
+
+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_checkalmostequal(q, Q, %eps,1e-4);
+
+// 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.00];
+ 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_checkalmostequal (q, Q, %eps,1e-4);
+
+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.00];
+ 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_checkalmostequal (q, Q, %eps,1e-4);
+
+test passed
+ 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.00];
+ 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_checkalmostequal (q, Q, %eps,1e-4);
+
+*/ \ No newline at end of file