対数正規分布の和の近似

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