対数正規分布の和の分布が必要なのだが、どうも簡単ではないらしい。 世の中にはいろいろな近似法が提案されているが、どれも一長一短らしい。 しかし、近似法は論文には載っているが、なにやらやたら複雑。 というわけで、いくつかの近似法を MATLAB で書いてみました。 ほんとに正しく実装できているか、とっても謎なので、間違いを見つけたら教えてください。
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
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
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