[x,y]=cmsbounds(m)

function [x,y]=cmsbounds(m)
% cmsbounds: calculate bounds on a distribution given a sequence of moments
% usage: [x,y]=cmsbounds(m)
%
% Given a sequence M of the first 2*n raw moments of a distribution,
% CMSBOUNDS calculates arrays X and Y of length n and n+1 respectively
% such that Y(i) <= F(X(i)) <= Y(i+1) for any distribution F with the given
% moments M (Chebyshev-Markov-Stieltjes inequalities).
%
% The raw moments are M(k) = E[X^k] = Integral of x^k dF(x).
%
% This version returns a warning if the moment sequence is inconsistent
% (e.g. m(2) < m(1)^2) or degenerate (e.g. 2n moments of a discrete
% distribution on less than 2n points).
%
% Example: Normal distribution
%  mu=0;
%  sig=1;
%  M(1)=mu;
%  M(2)=sig^2+mu^2;
%  for k=3:40
%    M(k)=mu*M(k-1)+(k-1)*sig^2*M(k-2);
%  end;
%  [x,y] = cmsbounds(M)  % returns 1x20 and 1x21 arrays
%  all(normcdf(x,mu,sig)<=y(2:end)) % true
%  all(normcdf(x,mu,sig)>=y(1:end-1)) % true
%
%
% Theory:
%
% CMSBOUNDS uses the Chebyshev-Markov-Stieltjes inequalities to calculate
% these bounds.  According to the theory, let F(x) be any probability
% measure (cumulative distribution function) satisfying E(x^k)=M(k) for
% 1<=k<=2*n and let Pk(x) be the unique set of orthogonal polynomials for
% k<=n such that E[Pj(x)*Pk(x)] = 1 if j=k and 0 otherwise.
%
% Let x1,...,xn be the zeros of Pn(x), and let
%   qn(x) = 1 / (P1(x)^2+...+P(n-1)(x)^2)
%
% The Chebyshev-Markov-Stieltjes theorm states that
%   qn(x1)+...+qn(x(k-1)) <= F(xk) <= qn(x1)+...+qn(xk)
%
% The difficult part is calculating the orthogonal polynomials, however we
% don't need to do it directly.  Use the Hankel matrix:
%   Hn = [1,M1,M2,...,Mn; M1,...,M(n+1); ... Mn,...,M(2n)];
% (for positive measures other than probability measures, replace the 1
% with F(infinity))
%
% Let Hn(x) be the matrix obtained by replacing the last row with
%   [1,x,...,x^n]
% It turns out that
%   Pn(x) = det(Hn(x))/sqrt(det(Hn)*det(H(n-1)).
%
% But it gets better: suppose Un is Hn reduced to upper triangular form.
% The set of coefficients are valid if no diagonal element of U are
% negative.
%
% We know that all systems of orthogonal polynomials satisfy
%   x*Pn(x) = bn*P(n+1)(x) + an*Pn(x) + b(n-1)*P(n-1)(x)
% and by partially expanding the determinants we find, after some algebra:
%   an = (U(n,n+1)/U(n,n) - U(n-1,n)/U(n-1,n-1))
%   bn = sqrt(U(n+1,n+1)/U(n,n))
%
% By writing the set of recursion relations as a linear system, we know
% that all roots of Pn(x) are eigenvalues of the tridiagonal matrix
% with diagonal a0,..., a(n-1) and offdiagonal both b0,...,b(n-2).
%
% The eigenvectors corresponding to xj are proportional to
% [P0(xj)=1, P1(xj), ..., P(n-1)(xj)].
%
%
%
% Implementation notes:
%
% CMSBOUNDS can theoretically accept as inputs any array objects that are
% compatible with doubles, e.g. fractions, vpi-based fractions,
% multiple-precision numbers etc.
%
% This can be useful if double precision is insufficient to determine the
% upper-triangular factor of Hn, although this operation does appear to be
% numerically stable.
%
% After the upper-triangular factor has been determined, the values are
% converted to compute the eigenvalues/vectors of the tridiagonal matrix
% (this part is numerically very stable).
%
% References:
%   http://en.wikipedia.org/wiki/Chebyshev-Markov-Stieltjes_inequalities
%   http://www.williams.edu/go/math/sjmiller/public_html/book/papers/jcmp.pdf
%   http://www.mathworks.com/matlabcentral/fileexchange/24878 (fractions toolbox)

% Author: Ben Petschel 21/8/2009
% Version history:
%   21/8/2009 - first release


% construct Hankel matrix
m=m(:);
n=floor(numel(m)/2); % H is n+1 x n+1

H = repmat(m(1,1),n+1,n+1); % ensures H is same object type as m
H(1,1)=1;
H(2:end,1)=m(1:n);
for k=2:n+1
  H(:,k)=m(k-1:k+n-1); % H(n+1,n+1)=m(n)
end;

% do row reduction on H
for k=1:n
  if H(k,k)<=0
    warning('cmsbounds:nonpositive','moment sequence is degenerate or inconsistent');
    x=[];
    y=[];
    return
  end;
  for j=k+1:n+1
    H(j,:)=H(j,:)-(H(j,k)/H(k,k))*H(k,:);
  end;
end;


% get the coefficients of the recursion for P0,P1,...,P(n-1)
% symmetrix tridiagonal matrix with diag elements ak, off-diag bk
A=zeros(n,n); % nxn matrix
A(1,1)=double(H(1,2)/H(1,1));
for k=2:n
  A(k,k)=double(H(k,k+1)/H(k,k)-H(k-1,k)/H(k-1,k-1));
  A(k,k-1)=sqrt(double(H(k,k)/H(k-1,k-1)));
  A(k-1,k)=A(k,k-1);
end;

[V,D]=eig(A); % eigenvalues are x-values; eigenvectors are P evaluated at x times constant
x=diag(D)'; % row vector
Pxsq=sum(V.^2)./(V(1,:).^2); % Px^2
y=[0,cumsum(1./Pxsq)];