From: <mtm...@us...> - 2012-05-04 23:46:59
|
Revision: 10360 http://octave.svn.sourceforge.net/octave/?rev=10360&view=rev Author: mtmiller Date: 2012-05-04 23:46:50 +0000 (Fri, 04 May 2012) Log Message: ----------- marcumq: improved rewrite added to signal package, submitted by Robert T. Short <rt...@ie...>. Closes artifact #3522119 Modified Paths: -------------- trunk/octave-forge/main/signal/NEWS Added Paths: ----------- trunk/octave-forge/main/signal/inst/marcumq.m Modified: trunk/octave-forge/main/signal/NEWS =================================================================== --- trunk/octave-forge/main/signal/NEWS 2012-05-04 09:09:56 UTC (rev 10359) +++ trunk/octave-forge/main/signal/NEWS 2012-05-04 23:46:50 UTC (rev 10360) @@ -4,6 +4,10 @@ signal-x.y.z Release Date: yyyy-mm-dd Release Manager: =============================================================================== + ** The function `marcumq' was imported from the communications package and has + been completely rewritten to improve performance and fix computational + errors. + ** Package is no longer automatically loaded ** The functions `__ellip_ws' and `__ellip_ws_min' have been removed (they Added: trunk/octave-forge/main/signal/inst/marcumq.m =================================================================== --- trunk/octave-forge/main/signal/inst/marcumq.m (rev 0) +++ trunk/octave-forge/main/signal/inst/marcumq.m 2012-05-04 23:46:50 UTC (rev 10360) @@ -0,0 +1,391 @@ +## Copyright (C) 2012 Robert T. Short <rt...@ie...> +## +## This 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 2 of the License, or +## (at your option) any later version. +## +## This software 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}, @var{m}, @var{tol}=eps) +## +## @noindent +## Compute the generalized Marcum Q function of order @var{M} with +## noncentrality parameter @var{a} and argument @var{b}. An optional +## relative tolerance @var{tol} may be included. +## +## @noindent +## If the input arguments are commensurate vectors, this function +## will produce a table of values (see tablify). +## +## @noindent +## 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. +## +## @noindent +## Reference: Marcum, "Tables of Q Functions", Rand Corporation. +## +## @noindent +## Reference: R.T. Short, "Computation of Noncentral Chi-squared +## and Rice Random Variables", www.phaselockedsystems.com/publications +## +## @seealso{imarcumq,tablify} +## +## @end deftypefn + +## Author: Robert T. Short (rt...@ie...) +## Description: Generalized Marcum Q function. +function [Q] = marcumq(a,b,M=1,tol=eps) + + if ( (nargin<2) || (nargin>4) ) + usage('[Q] = marcumq(a,b,M=1,tol=eps)'); + 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); This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |