-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathmarcumq.sci
147 lines (134 loc) · 5.1 KB
/
marcumq.sci
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
// Copyright (C) Scilab Enterprises - 20xx-2012 - Kartik PATEL <[email protected]>
//
// 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
function q_vec = marcumq(a_vec, b_vec, varargin)
// q = marcumq(a, b, m)
// This function calculates Marcum's Q function.
// Q_M(a,b) = 1/a^(m-1) \int_{b}^{\inf} x^m exp(-(x^2 + a^2)/2) I_{m-1} (ax) dx
// Where I_{m-1} (x) is modified bessel function.
//
// The calculation has been done using the infinite Bessel series. The infinte series has
// truncated when the relative error is been less than specified tolerance.
//
// Input Arguments:
// a_vec : Scalar or Vector or Matrix of values of a
// b_vec : Scalar or Vector or Matrix of values of b
// m : Scalar or Vector or Matrix of values of m
// Conditions on input:
// If size(a_vec) == size(b_vec) == size(m_vec) : The returned vector will be of size(m_vec)
// where q_vec(i) = marcumq(a_vec(i), b_vec(i), m_vec(i))
// If size of inputs a_vec, m_vec or b_vec is equal to [1 1] then
// it will be converted to size of other inputs that is not [1 1].
// If more than one input has size not equal to [1 1] and size of vectors are not same then
// error will be thrown.
// Output:
// q : Vector of size of inputs with values calculated by taking each individual elements
// from the inputs. i.e. q(i) = marcumq(a(i), b(i), m(i), tol)
//
// Reference:
// [1] R. T. Short, "Computation of Noncentral Chi-squared and Rice Random Variables",
// http://www.phaselockedsystems.com/NoncentralChiSquared.pdf
if argn(2)>3 | argn(2)<2 then
error(msprintf(gettext("Wrong number of Input argument\n")));
end
if or(imag(a_vec)~=0) then
error(msprintf(gettext("a must be non-negative real number\n")));
elseif or(a_vec < 0) then
error(msprintf(gettext("a must be non-negative real number\n")));
end
if or(imag(b_vec)~=0) then
error(msprintf(gettext("b must be non-negative real number\n")));
elseif or(b_vec < 0) then
error(msprintf(gettext("b must be non-negative real number\n")));
end
if size(a_vec) == [1 1] then
a_vec = repmat(a_vec, size(b_vec));
elseif size(b_vec) == [1 1] then
b_vec = repmat(b_vec, size(a_vec));
elseif ~isequal(size(a_vec), size(b_vec)) then
error(msprintf(gettext("Dimensions of input vectors must be same\n")));
end
if argn(2) == 3 then
m_vec = varargin(1);
if or(imag(m_vec)~=0) then
error(msprintf(gettext("m_vec must be positive integer.\n")));
elseif or(m_vec <= 0) then
error(msprintf(gettext("m_vec must be positive integer.\n")));
elseif or(floor(m_vec)~=m_vec) then
error(msprintf(gettext("m_vec must be positive integer.\n")));
end
if size(m_vec) == [1 1] then
m_vec = repmat(m_vec, size(a_vec));
else
if size(a_vec) == [1 1] then
a_vec = repmat(a_vec, size(m_vec));
b_vec = repmat(b_vec ,size(m_vec));
else
error(msprintf(gettext("Dimensions of input vectors must be same\n")));
end
end
else
m_vec = ones(size(a_vec));
end
tol = 1e-9;
//Variable check completed
q_vec = zeros(a_vec);
for i = 1:prod(size(a_vec))
b = b_vec(i);
a = a_vec(i);
m = m_vec(i);
if b==0 then
q = 1;
return;
end
if a==0 then
k = 0:m-1;
q = sum((b^(2.*k)./(2^(k).*factorial(k))))
q = exp(-(b^2)/2)*q;
return;
end
z = a * b;
if a < b then
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
k = k + 1;
while(1==1)
t = d*besseli(k,z,1);
S = S + t;
d = d*x;
k = k + 1;
if (abs(t/S) < tol) then
break;
end
end
q = S*exp(-((a-b)^2)/2);
else
x = b/a;
k = m;
d = x^m;
S = 0;
while(1==1)
t = d*besseli(k,z,1);
S = S + t;
d = d*x;
k = k + 1;
if (abs(t/S) < tol) then
break;
end
end
q = 1 - S*exp(-((a-b)^2)/2);
end
q_vec(i) = q;
end
endfunction