function [X, meanZ, sigmaZ, alpha, beta, error, ntry] = afosm2(p1, p2, p3, n_dis)

% AFOSM 拡張１次近似２次モーメント信頼性指標
%       programmed by Kohei TANAKA (tanaka@load.arch.t.u-tokyo.ac.jp)
%
% 書式：
%     [X, meanZ, sigmaZ, alpha, beta, error, ntry] = afosm(funcG, p1, p2, p3, n_dis);
%
% 実行例：
%     [X, meanZ, sigmaZ, alpha, beta, error, ntry] = afosm('X1/X2 - X3',...
%                                                  [40 25 20], [4 5 1], [0 0 0], [1 2 2])
%      X1/X2 - X3   X1:正規分布(40,4)   X2:対数正規分布(25,5)   X3:対数正規分布(20,1)
%
% 入力パラメータ:
%     funcG : 限界状態関数
%             確率変数の数は任意です。変数名は X1, X2, ... Xn です。
%             例） funcG = 'X1 - X2 / X3'
%
%     n_dis,p1,p2,p3 : 各確率変数の分布型とパラメータの与え方
%             1:正規分布  p1:平均値  p2:標準偏差  p3:ダミー変数(0)
%             2:対数正規分布  p1:平均値  p2:標準偏差  p3:ダミー変数(0)　　
%             3:ガンマ分布  p1:平均値  p2:標準偏差  p3:ダミー変数(0)
%             4:グンベル分布  p1:u1(特性最小値)  p2:α1  p3:ダミー変数(0)
%             5:フレッシェ分布  p1:k  p2:sigma  p3:mu(各変数の定義はmatlab help)
%             6:ワイブル分布  p1:k  p2:sigma  p3:mu(各変数の定義はmatlab help)
%
% 出力パラメータ:
%     X : 設計点
%     meanZ : 限界状態関数の平均値
%     sigmaZ : 限界状態関数の標準偏差
%     alfa : 係数α
%     beta : 拡張１次近似２次モーメント信頼性指標
%     error : 収束誤差
%     ntry : 試行回数
%
% この関数では、beta の反復計算を行い、１つ前の beta の値との差が
% 0.01 未満になったところで計算を打ち切ります。また試行回数が100回を超えると自動的に停止します。
%

% 変数の宣言
n = length(p1);                     % 変数の数
X = zeros(1,n);                     % 設計点
alpha = ones(1,n)*sqrt(1/n);         % 自乗和が1になるように初期値を設定
beta = 3.0;                         % 初期値は3.0
ntry = 0;
error = 1;

% 初期値の平均値と標準偏差
% 正規分布、対数正規分布、ガンマ分布についてはそのまま用いる
% 極値分布は変換する必要がある
meanX = zeros(1,n);
sigmaX = zeros(1,n);

for l = 1:n
    if n_dis(l) < 4
        meanX(l) = p1(l);
        sigmaX(l) = p2(l);
    elseif n_dis(l) == 4 % グンベル分布
        [meanX(l) V] = evstat(p1(l),1/p2(l));
        sigmaX(l) = sqrt(V);
    elseif n_dis(l) > 4 % フレッシェ分布、ワイブル分布
        [meanX(l) V] = gevstat(p1(l),p2(l),p3(l));
        sigmaX(l) = sqrt(V);
    end
end
        

% ループ
while error > 0.01 && ntry < 100
    
    % 設計点の計算
    X = meanX - alpha.* beta .* sigmaX;
    
    % 限界状態関数の計算
    [output_G grad_G] = cal_Func_G(X);
    
    % 任意の確率分布に対する正規化近似（平均値、標準偏差の変換）
    for m = 1:n
        
        switch n_dis(m)
            case 1 % 確率分布が正規分布
                meanX(m) = p1(m);  sigmaX(m) = p2(m);
            case 2 % 確率分布が対数正規分布
                [meanX(m) sigmaX(m)] = convert_logn(p1(m), p2(m), X(m));                
            case 3 % ガンマ分布
                [meanX(m) sigmaX(m)] = convert_gamma(p1(m), p2(m), X(m));    
            case 4 % 第1種極値分布（グンベル分布）
                [meanX(m) sigmaX(m)] = convert_gumbel(p1(m), p2(m), X(m));    
            case 5 % 第2種極値分布（フレッシェ分布）
                [meanX(m) sigmaX(m)] = convert_frecht(p1(m), p2(m), p3(m), X(m));    
            case 6 % 第3種極値分布（ワイブル分布）
                [meanX(m) sigmaX(m)] = convert_weibull(p1(m), p2(m), p3(m), X(m));    
        end
        
    end
                
    % 係数αの計算
    alpha = grad_G.*sigmaX./sqrt(sum(grad_G.^2.*sigmaX.^2));

    % 限界状態関数Zの平均値、標準偏差の計算
    meanZ = output_G + sum(grad_G.*(meanX - X));
    sigmaZ = sum(alpha .* grad_G .* sigmaX);

    % 信頼性指標βの算出
    new_beta = meanZ/sigmaZ;
    
    % エラーの計算
    error = abs((beta - new_beta)/beta);
    
    % βの更新
    beta = new_beta;
    
    ntry = ntry + 1;
    
end
    
return

% 限界状態関数の設定
function [output_G grad_G] = cal_Func_G(X)
    
    % 限界状態関数を与える 
    output_G = X(1)/X(2) - X(3);
    
    % 限界状態関数の偏微分を計算する
    % 前からX(1)の偏微分、X(2)の偏微分、X(3)の偏微分
    grad_G = [1/X(2), -X(1)/X(2)^2, -1];

return

% 正規化近似　対数正規分布
% 対数正規分布のパラメータはλとζ
% ζ=√log(1+(σ^2/μ^2))
% λ=logμ-0.5*ζ^2
function [meanX sigmaX] = convert_logn(mean, sigma, X)

% 変数変換
zeta = sqrt( log(1 + sigma.^2 / mean.^2) );
lambda = log(mean) - 0.5 * zeta^2;

sigmaX = normpdf(norminv(logncdf(X, lambda, zeta),0,1),0,1)/lognpdf(X, lambda, zeta);
meanX = X - norminv(logncdf(X, lambda, zeta),0,1)*sigmaX;

return

% 正規化近似　ガンマ分布
% ガンマ分布のパラメータはνとk
% μ=k/ν、σ^2=k/ν^2
% ν=σ^2/μ、k=μ^2/σ^2
function [meanX sigmaX] = convert_gamma(mean, sigma, X)

% 変数変換
k = mean^2/sigma^2;
nyu = sigma^2/mean;

sigmaX = normpdf(norminv(gamcdf(X, k, nyu)))/gampdf(X, k, nyu);
meanX = X - norminv(gamcdf(X, k, nyu))*sigmaX;

return

% 正規化近似　グンベル分布
% u1=p1 alpha1 = p2
% matlab内での変数の与え方
% MU=u1, sigma=1/p2
function [meanX sigmaX] = convert_gumbel(u1, alpha1, X)

sigmaX = normpdf(norminv(evcdf(X, u1, 1/alpha1)))/evpdf(X, u1, 1/alpha1);
meanX = X - norminv(evcdf(X, u1, 1/alpha1))*sigmaX;

return

% 正規化近似　フレッシェ分布
function [meanX sigmaX] = convert_frecht(k, sigma, mu, X)

sigmaX = normpdf(norminv(gevcdf(X, k, sigma, mu)))/gevpdf(X, k, sigma, mu);
meanX = X - norminv(gevcdf(X, k, sigma, mu))*sigmaX;

return

% 正規化近似　ワイブル分布
function [meanX sigmaX] = convert_weibull(k, sigma, mu, X)

sigmaX = normpdf(norminv(gevcdf(X, k, sigma, mu)))/gevpdf(X, k, sigma, mu);
meanX = X - norminv(gevcdf(X, k, sigma, mu))*sigmaX;

return