diff options
Diffstat (limited to 'macros/marcumq.sci')
-rw-r--r-- | macros/marcumq.sci | 423 |
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 |