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^2) = sqrt(pi/2) L(0.5,0.5*n-1,-0.5*l^2) v = n+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. 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

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^2) 0017 % = sqrt(pi/2) L(0.5,0.5*n-1,-0.5*l^2) 0018 % v = n+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 % 0023 % For n=1, the accuracy is high; for n>1, accuracy improves with increasing n. 0024 % Accuracy is worst when the non-centrality parameter, l, is close to s*sqrt(n). 0025 % Worst case errors as a function of n are: 0026 % n: 1 2 3 5 10 0027 % worst case error in m: 1e-15 0.007 0.004 0.0015 0.0005 0028 0029 % Copyright (C) Mike Brookes 2014 0030 % Version: $Id: v_chimv.m 4969 2014-08-05 18:24:30Z dmb $ 0031 % 0032 % VOICEBOX is a MATLAB toolbox for speech processing. 0033 % Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html 0034 % 0035 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0036 % This program is free software; you can redistribute it and/or modify 0037 % it under the terms of the GNU General Public License as published by 0038 % the Free Software Foundation; either version 2 of the License, or 0039 % (at your option) any later version. 0040 % 0041 % This program is distributed in the hope that it will be useful, 0042 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0043 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0044 % GNU General Public License for more details. 0045 % 0046 % You can obtain a copy of the GNU General Public License from 0047 % http://www.gnu.org/copyleft/gpl.html or by writing to 0048 % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. 0049 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0050 persistent ab pp qq nab 0051 if isempty(ab) 0052 nab=200; % cache a few low values of n 0053 pp=[ 0.595336298258636 -1.213013700592756 -0.018016200037799 1.999986150447582 0]; 0054 qq=[ -0.161514114798972 0.368983655790737 -0.136992134476950 -0.499681107630725 2]; 0055 ni=1./(1:nab); 0056 ab=[polyval(qq,ni);polyval(pp,ni)]; 0057 end 0058 if nargin<3 0059 s=1; 0060 if nargin<2 0061 l=0; 0062 end 0063 end 0064 ls=l/s; 0065 l2=(ls).^2; 0066 s2=s^2; 0067 if n<=nab 0068 if n==1 0069 m=l.*(1-2*normcdf(-ls))+2*s*normpdf(-ls); 0070 else 0071 m=sqrt(l2+n-1+(ab(1,n)+ab(2,n)*l2).^(-1))*s; 0072 end 0073 else 0074 m=sqrt(l2+n-1+(polyval(qq,1/n)+polyval(pp,1/n)*l2).^(-1))*s; 0075 end 0076 v=(n+l2)*s2-m.^2;

Generated on Wed 28-Mar-2018 15:33:24 by