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.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 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