V_CHIMV approximate mean and variance of non-central chi distribution [m,v]=(n,l,s)
Inputs: n = degrees of freedom
l = non-centrality parameter = sqrt(sum(mean^2)) [default 0]
(can be a vector or matrix to calculate many different values at once)
s = standard deviation of Gaussian [default 1]
Outputs: m = mean of chi distribution
v = variance of chi distribution
If x=c+randn(n,1) is a column vector of Gaussian random numbers with mean vector c, then
z=sqrt(x'*x) has a chi distributon with n degrees of freedom and non-centrality parameter
l=sqrt(c'*c). The mean and variance of a chi distribution are given precisely by
m = sqrt(2)*exp(gammaln(0.5*n+0.5)-gammaln(0.5*n))*hypergeom(-0.5,0.5*n,-0.5*(l/s)^2)*s
= sqrt(pi/2) L(0.5,0.5*n-1,-0.5*(l/s)^2)*s
v = n*s^2+l^2-m^2
where L(n,a,x) is the generalized Laguerre polynomial L_n^{(a)}(x) but this is very slow
to calculate so this routine approximates these expressions.
The expressions are from [1] together with the following identities taken
from [2] or [3]: 13.2.39, 13.6.19, 5.2.5, 5.4.1, 5.5.1, 5.4.6.
For n=1, the accuracy is high; for n>1, accuracy improves with increasing n.
Accuracy is worst when the non-centrality parameter, l, is close to s*sqrt(n).
Worst case errors as a function of n are:
n: 1 2 3 5 10
worst case error in m: 1e-15 0.007 0.004 0.0015 0.0005
References:
[1] J. H. Park. Moments of the generalized Rayleigh distribution. Quarterly of Applied Mathematics, 19: 45–49, 1961.
[2] F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark, editors. NIST Handbook of Mathematical Functions. CUP, 2010. ISBN 978-0-521-14063-8.
[3] F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark, editors. NIST Digital Library of Mathematical Functions. 2010-2014. URL http://dlmf.nist.gov/.0001 function [m,v]=v_chimv(n,l,s) 0002 %V_CHIMV approximate mean and variance of non-central chi distribution [m,v]=(n,l,s) 0003 % 0004 % Inputs: n = degrees of freedom 0005 % l = non-centrality parameter = sqrt(sum(mean^2)) [default 0] 0006 % (can be a vector or matrix to calculate many different values at once) 0007 % s = standard deviation of Gaussian [default 1] 0008 % 0009 % Outputs: m = mean of chi distribution 0010 % v = variance of chi distribution 0011 % 0012 % If x=c+randn(n,1) is a column vector of Gaussian random numbers with mean vector c, then 0013 % z=sqrt(x'*x) has a chi distributon with n degrees of freedom and non-centrality parameter 0014 % l=sqrt(c'*c). The mean and variance of a chi distribution are given precisely by 0015 % 0016 % m = sqrt(2)*exp(gammaln(0.5*n+0.5)-gammaln(0.5*n))*hypergeom(-0.5,0.5*n,-0.5*(l/s)^2)*s 0017 % = sqrt(pi/2) L(0.5,0.5*n-1,-0.5*(l/s)^2)*s 0018 % v = n*s^2+l^2-m^2 0019 % 0020 % where L(n,a,x) is the generalized Laguerre polynomial L_n^{(a)}(x) but this is very slow 0021 % to calculate so this routine approximates these expressions. 0022 % The expressions are from [1] together with the following identities taken 0023 % from [2] or [3]: 13.2.39, 13.6.19, 5.2.5, 5.4.1, 5.5.1, 5.4.6. 0024 % 0025 % For n=1, the accuracy is high; for n>1, accuracy improves with increasing n. 0026 % Accuracy is worst when the non-centrality parameter, l, is close to s*sqrt(n). 0027 % Worst case errors as a function of n are: 0028 % n: 1 2 3 5 10 0029 % worst case error in m: 1e-15 0.007 0.004 0.0015 0.0005 0030 % 0031 % References: 0032 % [1] J. H. Park. Moments of the generalized Rayleigh distribution. Quarterly of Applied Mathematics, 19: 45–49, 1961. 0033 % [2] F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark, editors. NIST Handbook of Mathematical Functions. CUP, 2010. ISBN 978-0-521-14063-8. 0034 % [3] F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark, editors. NIST Digital Library of Mathematical Functions. 2010-2014. URL http://dlmf.nist.gov/. 0035 0036 % Copyright (C) Mike Brookes 2014-2020 0037 % Version: $Id: v_chimv.m 11260 2020-07-18 20:07:58Z dmb $ 0038 % 0039 % VOICEBOX is a MATLAB toolbox for speech processing. 0040 % Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html 0041 % 0042 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0043 % This program is free software; you can redistribute it and/or modify 0044 % it under the terms of the GNU General Public License as published by 0045 % the Free Software Foundation; either version 2 of the License, or 0046 % (at your option) any later version. 0047 % 0048 % This program is distributed in the hope that it will be useful, 0049 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0050 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0051 % GNU General Public License for more details. 0052 % 0053 % You can obtain a copy of the GNU General Public License from 0054 % http://www.gnu.org/copyleft/gpl.html or by writing to 0055 % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. 0056 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0057 persistent ab pp qq nab 0058 if isempty(ab) 0059 nab=200; % cache a few low values of n 0060 pp=[ 0.595336298258636 -1.213013700592756 -0.018016200037799 1.999986150447582 0]; 0061 qq=[ -0.161514114798972 0.368983655790737 -0.136992134476950 -0.499681107630725 2]; 0062 ni=1./(1:nab); 0063 ab=[polyval(qq,ni);polyval(pp,ni)]; 0064 end 0065 if nargin<3 0066 s=1; 0067 if nargin<2 0068 l=0; 0069 end 0070 end 0071 ls=l/s; 0072 l2=(ls).^2; 0073 s2=s^2; 0074 if n<=nab 0075 if n==1 0076 m=l.*(1-2*normcdf(-ls))+2*s*normpdf(-ls); 0077 else 0078 m=sqrt(l2+n-1+(ab(1,n)+ab(2,n)*l2).^(-1))*s; 0079 end 0080 else 0081 m=sqrt(l2+n-1+(polyval(qq,1/n)+polyval(pp,1/n)*l2).^(-1))*s; 0082 end 0083 v=(n+l2)*s2-m.^2;