Home > voicebox > maxgauss.m

maxgauss

PURPOSE ^

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

SYNOPSIS ^

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

DESCRIPTION ^

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] = maxgauss(m,c,d)
0002 %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):145162, 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):15221533, Aug. 2007. doi: 10.1109/TCAD.2007.893544.
0029 
0030 %      Copyright (C) Mike Brookes 1997
0031 %      Version: $Id: maxgauss.m 713 2011-10-16 14:45:43Z 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 General Public License as published by
0039 %   the Free Software Foundation; either version 2 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 General Public License for more details.
0046 %
0047 %   You can obtain a copy of the GNU General Public License from
0048 %   http://www.gnu.org/copyleft/gpl.html or by writing to
0049 %   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
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

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