![[Makino's Page]](/~t/img/BannerLeft24.png)
![[Between Man and the Moon]](/~t/img/BannerRight.png)
Approximating sum of log-normals
It seems not easy to calculate the sum of (correlated) lognormal random variables. We can find various approximation methods, but they all have their pros and cons; and many of them are very complex. I tried to implement some of them in MATLAB. Be warned that I'm not confident that the implementation is correct. If you find anything wrong, please don't hesitate to tell me.
- FWSumLogNormal.m: Fenton-Wilkinson Approximation
- Simple, Quick. Good approximation for tail region, but usually bad for body region. Based on Pirinen's paper
- SYSumLogNormal.m: Ho's Approximation for Schwartz-Yeh Approximation
- Complex. Good approximation for body region, but bad for tail region. Based on Ho's paper
- MehtaSumLogNormal.m: Mehta et al.'s Approximation
- The focus of approximation can be chosen by parameter s. Currently the calculation complexity is exponential to the number of RVs. Uses fsolve in optimization toolbox. Based on Mehta et al.'s paper
FWSumLogNormal.m: Fenton-Wilkinson Approximation
function [ out ] = FWSumLogNormal( mu, sigma )
% mu: N-by-1 mean vector of the lognormal RVs
% sigma: N-by-N covariance matrix of the lognormal RVs
% ret: 2-by-1 vector, [mean; variance], of the merged lognormal
% Based on Fenton-Wilkinson Approximation
% Fenton, L.F. (1960). The sum of log-normal probability distibutions in
% scattered transmission systems.
% IEEE Trans. Commun. Systems 8: 57-67.
% Correlated version from
% Pekka Pirinen. STATISTICAL POWER SUM ANALYSIS FOR NONIDENTICALLY
% DISTRIBUTED CORRELATED LOGNORMAL SIGNALS
% Copyright (c) 2011 Takaki Makino.
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License http://www.gnu.org/licenses/gpl.html
% for more details.
sz = size(mu, 1);
musigmad = mu + diag(sigma) * 0.5 ;
%firstmoment = exp( musigmad );
%sumfm = sum(firstmoment);
% logsumfm = log(sumfm);
logsumfm = log(sum(exp( musigmad )));
%muxsig = repmat( musigmad, 1, size(sigma, 2) );
%secondmoment = exp( muxsig + muxsig' + 0.5*(sigma + sigma') );
%sumsm = sum( sum( secondmoment ) );
% logsumsm = log(sumsm);
logsumsm = log( sum( sum( exp( \
repmat( musigmad, 1, sz ) + repmat( musigmad', sz, 1 ) + 0.5*(sigma + sigma') \
) ) ) );
omu = 2*logsumfm - 0.5*logsumsm;
osigma2 = logsumsm - 2*logsumfm;
out = [omu; osigma2];
end
SYSumLogNormal.m: Ho's Approximation for Schwartz-Yeh Approximation
function [ ret ] = SYSumLogNormal( mu, sigma )
% mu: N-by-1 mean vector of the lognormal RVs
% sigma: N-by-N covariance matrix of the lognormal RVs
% ret: 2-by-1 vector, [mean; variance], of the merged lognormal
% Based on Ho's Approximation for Schwartz-Yeh Approximation
% IEEE TRANSACTIONS ON VEHICULAR TECHNOLOGY, VOL. 44. NO. 4, NOVEMBER 1995
% Calculating the Mean and Variance of Power Sums with Two Log-Normal Component
% Chia-Lu Ho
% Copyright (c) 2011 Takaki Makino.
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License http://www.gnu.org/licenses/gpl.html
% for more details.
% prepare constants used in trape_integral
delta = 0.3;
xval = -3.6:delta:3.6;
uval = exp( -2 * sinh(xval) );
vval = 1 ./ (1+uval);
zval = (2 * delta) * vval .* uval .* cosh(xval);
% sort inputs
[sortmu, idx] = sort(mu, 1, 'descend');
sortsigma = sigma(idx, idx);
% merge two-by-two
while size(sortmu,1) >= 2
mw = sortmu(2) - sortmu(1);
sigy12 = sortsigma(1,1) ;
sigy1 = sqrt( sortsigma(1,1) );
sigy22 = sortsigma(2,2) ;
sigy2 = sqrt( sortsigma(2,2) );
rhoy1y2sigy1sigy2 = sortsigma(1,2) ;
sigw2 = sigy12 + sigy22 - 2*rhoy1y2sigy1sigy2;
sigw = sqrt( sigw2 );
% rhoy1y2 = rhoy1y2sigy1sigy2 / (sigy1 * sigy2)
eta = - mw / sigw;
fw = @(w)( exp(-(w - mw).^2/(2*sigw2) ) / sqrt(2*pi*sigw2) );
h0 = @(v)( exp( -0.5*(-log(v)+eta).^2 ) * sqrt(0.5/pi) );
h1 = @(v)( (fw(log(v))+fw(-log(v))).* log(1+v) );
h2 = @(v)( (fw(log(v))-fw(-log(v)))./ (1+1./v) );
h3 = @(v)( (fw(log(v))+fw(-log(v))).* log(1+v).^2 );
h4 = @(v)( -fw(-log(v)) .* log(v) .* log(1+v) );
h5 = @(v)( fw(-log(v)) ./ (1+1./v) );
h6 = @(v)( fw(-log(v)) .* log(1+v) );
h1234_old = @(v)( [h1(v); h2(v); h3(v); h4(v)] );
%hall = @(v)( [h1(v); h2(v); h3(v); h5(v); h6(v) ] );
I0 = erfc( sqrt(eta) );
I = h1234(vval) * zval';
% I = trape_integral( @h1234 );
% I1 = I(1);
% I2 = I(2);
% I3 = I(3);
% I4 = I(4);
% I5 = Iall(4);
% I6 = Iall(5);
% I4 = sigw2 * (fw(0) * log(2) - I5) + mw * I6;
A0 = sigw / sqrt(2 * pi) * exp( - eta^2 / 2 ) + mw * I0;
G1 = A0 + I(1);
G2 = I(3) + 2*I(4) + sigw2 * I0 + mw * A0;
G3 = sigw2 * (I(2) + I0);
mz = sortmu(1) + G1;
sigz2 = sigy12 - G1^2 + G2 + 2*(rhoy1y2sigy1sigy2 - sigy12)*(G3 / sigw2);
sigz = sqrt( sigz2 );
diag3 = diag(sortsigma); diag3 = sqrt(diag3(3:size(diag3,1)));
rhos = sortsigma(3:size(sortsigma,1),1) ./ diag3;
rhosy = sortsigma(3:size(sortsigma,1),2) ./ diag3;
newrhos = rhos + (rhosy - rhos) * (G3 / sigw2);
newsig = newrhos .* diag3;
newmu = vertcat( mz, sortmu(3:size(sortmu)));
newsig0 = sortsigma( 3:size(sortsigma,1), 3:size(sortsigma,2));
newsig1 = horzcat( sigz2, newsig');
newsig2 = horzcat( newsig, newsig0);
newsigma = vertcat( newsig1, newsig2 );
sortmu = newmu;
sortsigma = newsigma;
end
ret = [ sortmu ; sortsigma ];
%function [out] = trape_integral( func )
% hval = func( vval );
% out = hval * zval';
%end
function [out] = h1234(v)
logv = log(v);
fwlogv = fw(logv);
fwmlogv = fw(-logv);
log1v = log(1+v);
h1_ = (fwlogv+fwmlogv).* log1v;
out = [ h1_ ;
(fwlogv-fwmlogv)./ (1+1./v) ;
h1_ .* log1v ;
-fwmlogv .* log1v .* logv ];
end
end
MehtaSumLogNormal.m: Mehta et al.'s approximation
function [ xopt ] = MehtaSumLogNormal( mu, sigma, s )
% mu: N-by-1 mean vector of the lognormal RVs
% sigma: N-by-N covariance matrix of the lognormal RVs
% s: 2-by-1 vector that specifies the points to be matched (default [0.2;1.0])
% ret: 2-by-1 vector, [mean; variance], of the merged lognormal
% Based on Mehta et al.'s approximation
% "Approximating a Sum of Random Variables with a Lognormal"
% Neelesh Mehta, Jingxian Wu, Andreas Molisch, Jin Zhang
% IEEE TRANSACTIONS ON WIRELESS COMMUNICATIONS, VOL. 6, NO. 7, pp.2690-2699, JULY 2007
% doi:10.1109/TWC.2007.051000.
% Copyright (c) 2011 Takaki Makino.
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License http://www.gnu.org/licenses/gpl.html
% for more details.
% prepare values used in HermiteIntegration
% Values taken from Table 25.10 of
% M. Abramowitz and I. Stegun, Handbook of mathematical functions with formulas,
% graphs, and mathematical tables. Dover, 9 ed., 1972.
% 12-point approximation
awarray = [
0.314240376254259, 5.701352352525e-1 ;
-0.314240376254259, 5.701352352525e-1 ;
0.947788391240164, 2.604923102642e-1 ;
-0.947788391240164, 2.604923102642e-1 ;
1.597682635152605, 5.160798561588e-2 ;
-1.597682635152605, 5.160798561588e-2 ;
2.279507080501060, 3.905390584629e-3 ;
-2.279507080501060, 3.905390584629e-3 ;
3.020637025120890, 8.573687043588e-5 ;
-3.020637025120890, 8.573687043588e-5 ;
3.889724897869782, 2.658551684356e-7 ;
-3.889724897869782, 2.658551684356e-7 ];
% May use cheaper approximation; replace awarray
% 5-point approximation
% awarray = [
% 0, 9.453087204829e-1 ;
% 0.958572464613819, 3.936193231522e-1 ;
% -0.958572464613819, 3.936193231522e-1 ;
% -2.020182870456086, 1.995324205905e-2 ;
% 2.020182870456086, 1.995324205905e-2 ];
% 8-point approximation
% awarray = [
% 0.381186990207322, 6.611470125582e-1 ;
% -0.381186990207322, 6.611470125582e-1 ;
% 1.157153712446780, 2.078023258149e-1 ;
% -1.157153712446780, 2.078023258149e-1 ;
% 1.981656756695843, 1.707798300741e-2 ;
% -1.981656756695843, 1.707798300741e-2 ;
% 2.930637420257244, 1.996040722114e-4 ;
% -2.930637420257244, 1.996040722114e-4 ];
a = awarray(:,1)' ;
w = awarray(:,2)' ;
N = size(awarray, 1);
if nargin < 3, s = [0.2; 1.0]; end
% Perform Hermite Integration to find MGFs
val = MultiHermiteIntegration( s, mu, sigma );
% Determine initial point with a quick approximation
x0 = FWSumLogNormal( mu, sigma );
f = @( x ) ( SingleHermiteIntegration(s, x(1), x(2)) - val );
option = optimset('Display', 'off');
ret = fsolve( f, x0, option );
function [ ret ] = SingleHermiteIntegration( s, mu, sigma )
%s1 = repmat(-s, [N,1])
%s2 = repmat( exp( (sqrt(2* sigma) * a + mu ) ), size(s) )
%s3 = exp(s1 .* s2)
%s4 = repmat(w, size(s))
%s5 = s4 .* s3
% ret = sum(s5)
ret = sum( repmat(w, size(s)) .* \
exp( repmat(-s, [1, N]) .* \
repmat( exp( (sqrt(2*sigma) * a + mu ) ), size(s) ) ), 2 ) \
/ sqrt(pi);
end
function [ ret ] = MultiHermiteIntegration( s, mu, sigma )
K = prod(size(mu));
sums = reshape(mu,[K,1]);
dim1n = [N];
dimnn = [1];
dimkn = [K];
[U, Lambda] = eig(sigma);
C = U * sqrt(Lambda);
for j=1:K
dim1n = horzcat( 1, dim1n );
dimnn = horzcat( dimnn, N );
anj = repmat( reshape(a, dim1n), [dimkn, 1] );
ckj = repmat( C(:,j), dimnn );
sums = repmat( sums, dim1n ) + sqrt(2) * ckj .* anj;
dimkn = horzcat( dimkn, N );
dimsn = horzcat( dimkn, N );
end
% sums is KxNxNx...N matrix
sums = exp( repmat(-s, dimnn) .* repmat( sum( exp( sums ), 1 ), size(s) ) );
while( size(sums, 2) > 1 )
sizesums = size(sums);
dimsn = [ sizesums( 1:ndims(sums)-1 ), 1];
dim1n = [ ones( 1, ndims(sums)-1), N ];
weight = repmat( reshape(w/sqrt(pi), dim1n), dimsn );
sums = sum( sums .* weight, ndims(sums) );
end
ret = sums;
end
end