Home > voicebox > pdfmoments.m

pdfmoments

PURPOSE ^

PDFMOMENTS convert between central moments, raw moments and cumulants [C,R,K]=(T,M,B,A)

SYNOPSIS ^

function [c,r,k]=pdfmoments(t,m,b,a)

DESCRIPTION ^

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) pdfmoments('k',[m,s^2,0,0,0])=pdfmoments('k',[0,1,0,0,0],m,s) gives a Gaussian with mean m and std dev s

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [c,r,k]=pdfmoments(t,m,b,a)
0002 %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) pdfmoments('k',[m,s^2,0,0,0])=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: pdfmoments.m 9678 2017-04-04 14:28:37Z 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

Generated on Tue 10-Oct-2017 08:30:10 by m2html © 2003