# v_pdfmoments

## PURPOSE

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

## SYNOPSIS

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

## DESCRIPTION

    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 s

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

## SOURCE CODE

    0001 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 Lesser General Public License as published by
    0025 %   the Free Software Foundation; either version 3 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 Lesser General Public License for more details.
    0032 %
    0033 %   You can obtain a copy of the GNU Lesser General Public License from
    0034 %   https://www.gnu.org/licenses/ .
    0035 %
    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