# v_maxgauss

## PURPOSE

V_MAXGAUSS determine gaussian approximation to max of a gaussian vector [p,u,v,r]=(m,c,d)

## SYNOPSIS

function [u,v,p,r] = v_maxgauss(m,c,d)

## DESCRIPTION

    V_MAXGAUSS determine gaussian approximation to max of a gaussian vector [p,u,v,r]=(m,c,d)

     Inputs:
           m(N,1) is the mean vector of length N
           c(N,N) is Covariance matrix (or c(N,1) is vector of covariances)
           d(N,K) is Covariance w.r.t some other variables

     Outputs:
           u is mean(max(x)) where x is the random variable of length N
           v is var(max(x))
           p(N,1) is prob of each element being the max
           r(1,K) is covariance between max(x) and the variables corresponding to the columns of d

     To find the min instead of max, just negate m and u.

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

## SOURCE CODE

    0001 function [u,v,p,r] = v_maxgauss(m,c,d)
    0002 %V_MAXGAUSS determine gaussian approximation to max of a gaussian vector [p,u,v,r]=(m,c,d)
    0003 %
    0004 % Inputs:
    0005 %       m(N,1) is the mean vector of length N
    0006 %       c(N,N) is Covariance matrix (or c(N,1) is vector of covariances)
    0007 %       d(N,K) is Covariance w.r.t some other variables
    0008 %
    0009 % Outputs:
    0010 %       u is mean(max(x)) where x is the random variable of length N
    0011 %       v is var(max(x))
    0012 %       p(N,1) is prob of each element being the max
    0013 %       r(1,K) is covariance between max(x) and the variables corresponding to the columns of d
    0014 %
    0015 % To find the min instead of max, just negate m and u.
    0016
    0017 % The algorithm combines the elements of m in pairs in sequence as suggested
    0018 % in [1]. Errors are introduced because max(x(i),x(j)) is wrongly assumed to be
    0019 % gaussian. To minimize the errors, we use a greedy algorithm (approximately as
    0020 % in [2]) in which the chosen pair has the greatest imbalance in selection probability.
    0021 %
    0022 %
    0023 % Refs: [1] C. E. Clark. The greatest of a finite set of random variables.
    0024 %           Operations Research, 9(2):145�162, March 1961.
    0025 %       [2] D. Sinha, H. Zhou, and N. V. Shenoy.
    0026 %           Advances in Computation of the Maximum of a Set of Gaussian Random Variables.
    0027 %           IEEE Trans Computer-Aided Design of Integrated Circuits and Systems,
    0028 %           26(8):1522�1533, Aug. 2007. doi: 10.1109/TCAD.2007.893544.
    0029
    0030 %      Copyright (C) Mike Brookes 1997
    0031 %      Version: $Id: v_maxgauss.m 10865 2018-09-21 17:22:45Z dmb $
    0032 %
    0033 %   VOICEBOX is a MATLAB toolbox for speech processing.
    0034 %   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
    0035 %
    0036 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0037 %   This program is free software; you can redistribute it and/or modify
    0038 %   it under the terms of the GNU Lesser General Public License as published by
    0039 %   the Free Software Foundation; either version 3 of the License, or
    0040 %   (at your option) any later version.
    0041 %
    0042 %   This program is distributed in the hope that it will be useful,
    0043 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
    0044 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    0045 %   GNU Lesser General Public License for more details.
    0046 %
    0047 %   You can obtain a copy of the GNU Lesser General Public License from
    0048 %   https://www.gnu.org/licenses/ .
    0049 %
    0050 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0051
    0052 nrh=-sqrt(0.5);
    0053 kpd=sqrt(0.5/pi);
    0054 m=m(:);
    0055 nm=length(m);
    0056 if nargin<2
    0057     c=eye(nm);
    0058 elseif numel(c)==nm
    0059     c=diag(c);
    0060 end
    0061 p=eye(nm);
    0062 ix=find(m==-Inf);       % remove negative infinities immediately
    0063 if ~isempty(ix)
    0064     m(ix)=[];
    0065     c(ix,:)=[];
    0066     c(:,ix)=[];
    0067     p(:,ix)=[];
    0068     nm=length(m);
    0069 end
    0070 ix=0:nm^2-1;                        % index for nm x nm matrix elements (from 0)
    0071 ix = ix(ix<(nm+1)*floor(ix/nm));    % only keep the strict upper triangle
    0072 mij=zeros(2,1);                     % store for means
    0073 while nm>1
    0074
    0075     % calculate scaled differences in means
    0076
    0077     gm=m(:,ones(1,nm));
    0078     gm = gm-gm.';               % find the difference between all pairs of means
    0079     gm=gm(ix+1);                % only keep the strict upper triangular elements
    0080     cd=diag(c);
    0081     gv=cd(:,ones(1,nm))-c;
    0082     gv=gv+gv.';
    0083     gv=gv(ix+1);                % These are the corresponding variance sums
    0084     jx=find(gv<=0);
    0085     if ~isempty(jx)             % special case: two variables differ by a constant
    0086         jx=jx(1);               % take first pair for which this is true
    0087         j=floor(ix(jx)/nm);
    0088         i=ix(jx)-nm*j+1;
    0089         j=j+1;
    0090         dm=gm(jx);
    0091         if dm>0                 % if x(i)>x(j) then abolish j
    0092             m(j)=[];
    0093             c(j,:)=[];
    0094             c(:,j)=[];
    0095             p(:,j)=[];
    0096         else                    % if x(i)<=x(j) then abolish i
    0097             m(i)=[];
    0098             c(i,:)=[];
    0099             c(:,i)=[];
    0100             p(:,i)=[];
    0101         end
    0102     else
    0103         % select the pair of variables with the the highest ratio of
    0104         % squared difference in means to sum of variances and combine into a single variable
    0105         [gg,jx]=max(gm.^2./gv);         % jx indicates which pair to combine
    0106         j=floor(ix(jx)/nm);             % convert jx into (i,j) pair
    0107         i=ix(jx)-nm*j+1;
    0108         j=j+1;                          % combine variables i and j
    0109         dm=gm(jx);                      % mean of x(i)-x(j)
    0110         ds=sqrt(gv(jx));                % std dev of x(i)-x(j)
    0111         dms=dm/ds;
    0112         q = 0.5 * erfc(nrh*dms);        % q =normcdf(dm,0,ds) = Phi(dms) = prob x(i) > x(j)
    0113         f=kpd*exp(-0.5*dms^2);          % f=phi(dms)=ds*normpdf(dm,0,ds)
    0114         mij(1)=m(i);
    0115         mij(2)=m(j);
    0116         u=dm*q+mij(2)+ds*f;                                         % mean of max{x(i),x(j)}
    0117         v=(mij(1)+mij(2)-u)*u +cd(i)*q+cd(j)*(1-q)-mij(1)*mij(2);   % variance of max{x(i),x(j)}
    0118
    0119         % replace x(i) with max{x(i),x(j)}
    0120
    0121         m(i)=u;
    0122         c(i,:)=q*c(i,:)+(1-q)*c(j,:);
    0123         c(:,i)=c(i,:).';
    0124         c(i,i)=v;
    0125         p(:,i)=q*p(:,i)+(1-q)*p(:,j);
    0126
    0127         % now abolish x(j)
    0128
    0129         m(j)=[];
    0130         c(j,:)=[];
    0131         c(:,j)=[];
    0132         p(:,j)=[];
    0133     end
    0134
    0135     ix=ix(1:(nm-1)*(nm-2)/2);
    0136     ix=ix-floor(ix/nm);
    0137     nm=nm-1;
    0138 end         % main while loop
    0139 u=m(1);
    0140 v=c(1);
    0141 p=p/sum(p);     % force sum=1
    0142 if nargin>2
    0143     r=p.'*d;
    0144 end