V_PDFMOMENTS convert between central moments, raw moments and cumulants [C,R,K]=(T,M,B,A)
Inputs: t text string containing:
'm','r','k' Imput is central moments, raw moments, cumulants [default 'm']
'M','R','K' Ouptut c is central moments, raw moments, cumulants [default 'M']
m vector of input moments; m(r) is moment r, m(1) is always the mean
a,b If input moments are for x, output moments are for a*x+b [defaulta=1, b=0]
Outputs: c central moments (or as determined by 'R' or 'K' options)
r raw moments
k cumulants
(a) For all formats, the first element is the mean (i.e. not the first central moment or cumulant which equal zero).
(b) v_pdfmoments('k',[m,s^2,0,0,0])=v_pdfmoments('k',[0,1,0,0,0],m,s) gives a Gaussian with mean m and std dev s0001 function [c,r,k]=v_pdfmoments(t,m,b,a) 0002 %V_PDFMOMENTS convert between central moments, raw moments and cumulants [C,R,K]=(T,M,B,A) 0003 % Inputs: t text string containing: 0004 % 'm','r','k' Imput is central moments, raw moments, cumulants [default 'm'] 0005 % 'M','R','K' Ouptut c is central moments, raw moments, cumulants [default 'M'] 0006 % m vector of input moments; m(r) is moment r, m(1) is always the mean 0007 % a,b If input moments are for x, output moments are for a*x+b [defaulta=1, b=0] 0008 % 0009 % Outputs: c central moments (or as determined by 'R' or 'K' options) 0010 % r raw moments 0011 % k cumulants 0012 % 0013 % (a) For all formats, the first element is the mean (i.e. not the first central moment or cumulant which equal zero). 0014 % (b) v_pdfmoments('k',[m,s^2,0,0,0])=v_pdfmoments('k',[0,1,0,0,0],m,s) gives a Gaussian with mean m and std dev s 0015 0016 % Copyright (c) 1998 Mike Brookes, mike.brookes@ic.ac.uk 0017 % Version: $Id: v_pdfmoments.m 10865 2018-09-21 17:22:45Z dmb $ 0018 % 0019 % VOICEBOX is a MATLAB toolbox for speech processing. 0020 % Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html 0021 % 0022 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0023 % This program is free software; you can redistribute it and/or modify 0024 % it under the terms of the GNU General Public License as published by 0025 % the Free Software Foundation; either version 2 of the License, or 0026 % (at your option) any later version. 0027 % 0028 % This program is distributed in the hope that it will be useful, 0029 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0030 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0031 % GNU General Public License for more details. 0032 % 0033 % You can obtain a copy of the GNU General Public License from 0034 % http://www.gnu.org/copyleft/gpl.html or by writing to 0035 % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. 0036 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0037 persistent n0 bc mk fa 0038 if isempty(n0) 0039 n0=3; % order 3 0040 bc=[1 1 0 0; 1 2 1 0; 1 3 3 1]; % binomial coefficients 0041 mk={[1 1 1];[0 1 1 1]}; % cumulant coefficients: one row per term, powers of k(2:i+1), coef k->m, coef m->k 0042 mn=[1 0; 1 1]; % [i,j] = number of terms in m(i+1) whose lowest moment is >= j+1 0043 fa=1; % factorial list 0044 end 0045 % check arguments 0046 if nargin<4 || isempty(a) 0047 a=1; 0048 end 0049 if nargin<3 || isempty(b) 0050 b=0; 0051 end 0052 if isempty(t) 0053 t=''; 0054 end 0055 n=length(m); % number of moments required 0056 if n>n0 % check if need to update coefficient arrays 0057 if fix(n/2-1)>length(fa) % we need factorials up to fix(n/2-1) 0058 fal=length(fa); 0059 fa(fix(n/2-1),1)=0; % enlarge factorial vector 0060 for i=fal+1:fix(n/2-1) 0061 fa(i)=i*fa(i-1); % create new factorials 0062 end 0063 end 0064 bc(n,n+1)=0; % enlarge binomial coefficient array 0065 mk{n-1,1}=[]; % enlarge cumulant coefficient array 0066 mn(n-1,n-1)=0; % enlarge cumulant coefficient counts 0067 for i=n0+1:n 0068 bc(i,1:i+1)=[1 bc(i-1,1:i)+bc(i-1,2:i+1)]; % update binomial coefficients 0069 j=fix((i+1)/2); % first coefficient row to sum 0070 nr=1+sum(mn(((j-1:i-3)+(n-1)*(i-j-2:-1:0)))); % number of terms 0071 mki=zeros(nr,i+1); % coefficient matrix 0072 ix=1; 0073 mki(1,i-1:i+1)=1; % first term always has a coefficient of 1 0074 for r=j:i-2 % previous coefficients to use 0075 nk=mn(r-1+(n-1)*(i-r-2)); % number of new coefficients for this value of r 0076 mkk=mk{r-1}; % old coefficients for this value of r 0077 mkik=mkk(1:nk,1:r-1); % extract just the list of powers for each term 0078 mkik(:,i-r-1)=mkik(:,i-r-1)+1; % increment the power of moment i-r 0079 mki(ix+1:ix+nk,1:r-1)=mkik; % and save as new terms 0080 mki(ix+1:ix+nk,i)=mkk(1:nk,r)*bc(i,i-r+1)./mkik(:,i-r-1); % calculate coefficient for r->m 0081 rho=sum(mkik,2)-1; % rho is one less than the sum of the moment powers 0082 mki(ix+1:ix+nk,i+1)= mki(ix+1:ix+nk,i).*fa(rho).*(-1).^rho; % calculate coefficient for m->r 0083 ix=ix+nk; % update the number of terms so far 0084 end 0085 mki=sortrows(mki); % sort according to the lowest moment that is used 0086 mn(i-1,1:i-1)=[nr sum(cumprod(mki(:,1:i-2)==0,2),1)]; % update count of terms with lowest moment >= j+1 0087 mk{i-1}=mki; % save in persistent cell array 0088 end 0089 n0=n; % coefficients are now calculated up to order n 0090 end 0091 % apply scaling if input type is 'c' or 'k' 0092 mu=a*m(1)+b; % calculate new mean 0093 c=m; % initialize output shapes 0094 r=m; 0095 k=m; 0096 m=m(:)'; % now force the input to be a row vector 0097 if any(t=='k') 0098 tin=3; % set input type 0099 k(:)=k(:)'.*a.^(1:n); 0100 k(1)=0; % first cumulant is actually zero 0101 elseif any(t=='r') 0102 tin=2; 0103 else 0104 tin=1; 0105 c(:)=c(:)'.*a.^(1:n); 0106 c(1)=0; % first cenral moment is actually zero 0107 end 0108 tout=[(~any(t=='K') && ~any(t=='R')) (nargout>=2 || any(t=='R')) (nargout>=3 || any(t=='K'))]; % outputs required 0109 for il=1:2 % loop through conversion routines twice 0110 % first convert between moments 0111 if il==1 % convert unscaled R -> C 0112 v=[1 m.*a.^(1:n)]; 0113 bb=b-mu; 0114 doit=tin==2 && (tout(1) || tout(3)); 0115 else % convert C -> R or unscaled R -> R 0116 if tin==2 % input type was 'r' (v is OK from previous iteration) 0117 bb=b; 0118 else % input type was 'c' or 'k' 0119 v=[1 c(:)']; 0120 bb=mu; 0121 end 0122 doit=tout(2); % convert if 'R' output required 0123 end 0124 if doit 0125 y=v(2:end); 0126 if bb~=0 % don't bother if the constant term is zero 0127 for i=1:n 0128 y(i)=polyval(bc(i,1:i+1).*v(1:i+1),bb); 0129 end 0130 end 0131 if il==1 % convert unscaled R -> C 0132 c(:)=y; 0133 else % convert C -> R or unscaled R -> R 0134 r(:)=y; 0135 end 0136 end 0137 % now convert cumulants to/from moments 0138 if il==1 % convert K -> C 0139 x=k(:)'; 0140 doit=tin==3 && (tout(1) || tout(2)); 0141 else % convert C -> K 0142 x=c(:)'; 0143 doit=(tin<3) && tout(3); 0144 end 0145 if doit 0146 y=x; 0147 for i=4:n 0148 mki=mk{i-1}; % get coefficient matrix 0149 y(i)=mki(:,i-1+il)'*prod(repmat(x(2:i),size(mki,1),1).^mki(:,1:i-1),2); % calculate moment/cumulant (neat but not efficient) 0150 end 0151 if il==1 % converted K -> C 0152 c(:)=y; 0153 else % converted C -> K 0154 k(:)=y; 0155 end 0156 end 0157 end 0158 c(1)=mu; % restore the means 0159 k(1)=mu; 0160 if any(t=='R') 0161 c=r; 0162 elseif any(t=='K') 0163 c=k; 0164 end