Home > voicebox > v_chimv.m

v_chimv

PURPOSE ^

V_CHIMV approximate mean and variance of non-central chi distribution [m,v]=(n,l,s)

SYNOPSIS ^

function [m,v]=v_chimv(n,l,s)

DESCRIPTION ^

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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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 Tue 19-Sep-2017 12:07:31 by m2html © 2003