![[Makino's Page]](/~t/img/xBannerLeft24.png.pagespeed.ic.yx6lsWjDp1.png)
![[Between Man and the Moon]](/~t/img/xBannerRight.png.pagespeed.ic.DOFE24z6u6.png)
対数正規分布の和の近似
対数正規分布の和の分布が必要なのだが、どうも簡単ではないらしい。 世の中にはいろいろな近似法が提案されているが、どれも一長一短らしい。 しかし、近似法は論文には載っているが、なにやらやたら複雑。 というわけで、いくつかの近似法を MATLAB で書いてみました。 ほんとに正しく実装できているか、とっても謎なので、間違いを見つけたら教えてください。
- FWSumLogNormal.m: Fenton-Wilkinson Approximation
- 1次、2次モーメントをあわせるだけ。単純で高速、tail 領域 (大きい変数値に対する確率分布) ではよい近似となることが知られていますが、 メインの部分ではだいぶずれます。 Based on Pirinen's paper
- SYSumLogNormal.m: Ho's Approximation for Schwartz-Yeh Approximation
- log 領域での1次、2次モーメントを合わせます。複雑。つーかなにやってるのか理解してません。一般的に FW 法よりよいといわれていますが、tail 領域では確率を過小評価する傾向があるようです。 Based on Ho's paper
- MehtaSumLogNormal.m: Mehta et al.'s Approximation
- MGF を2つ計算して合わせます。パラメータ s によって、近似をよくする場所を選べます。確率変数の数の指数オーダーの計算時間。Optimization toolbox の fsolve を使用。 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