mode_and_variance_to_mean

 This function computes the mean of a distribution given the mode and variance of this distribution.

  INPUTS 
    m                [double]    scalar, mode of the distribution.
    s2               [double]    scalar, variance of the distribution.
    distribution     [integer]   scalar for the distribution shape 
                                    1 gamma
                                    2 inv-gamma-2
                                    3 inv-gamma-1
                                    4 beta    
    lower_bound      [double]    scalar, lower bound of the random variable support (optional).
    upper_bound      [double]    scalar, upper bound of the random variable support (optional).
    
  OUTPUT 
    mu               [double]    scalar, mean of the distribution.
    parameters       [double]    2*1 vector, parameters of the distribution.

0001 function [mu, parameters] = mode_and_variance_to_mean(m,s2,distribution,lower_bound,upper_bound)
0002 % This function computes the mean of a distribution given the mode and variance of this distribution.
0003 %
0004 %  INPUTS
0005 %    m                [double]    scalar, mode of the distribution.
0006 %    s2               [double]    scalar, variance of the distribution.
0007 %    distribution     [integer]   scalar for the distribution shape
0008 %                                    1 gamma
0009 %                                    2 inv-gamma-2
0010 %                                    3 inv-gamma-1
0011 %                                    4 beta
0012 %    lower_bound      [double]    scalar, lower bound of the random variable support (optional).
0013 %    upper_bound      [double]    scalar, upper bound of the random variable support (optional).
0014 %
0015 %  OUTPUT
0016 %    mu               [double]    scalar, mean of the distribution.
0017 %    parameters       [double]    2*1 vector, parameters of the distribution.
0018 %
0019 
0020 % Copyright (C) 2009 Dynare Team
0021 %
0022 % This file is part of Dynare.
0023 %
0024 % Dynare is free software: you can redistribute it and/or modify
0025 % it under the terms of the GNU General Public License as published by
0026 % the Free Software Foundation, either version 3 of the License, or
0027 % (at your option) any later version.
0028 %
0029 % Dynare is distributed in the hope that it will be useful,
0030 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0031 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0032 % GNU General Public License for more details.
0033 %
0034 % You should have received a copy of the GNU General Public License
0035 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0036 
0037 % Check input aruments.
0038 if ~(nargin==3 || nargin==5 || nargin==4 )
0039     error('mode_and_variance_to mean:: 3 or 5 input arguments are needed!')
0040 end
0041 
0042 % Set defaults bounds.
0043 if nargin==3
0044     switch distribution
0045       case 1
0046         lower_bound = 0;
0047         upper_bound = Inf;
0048       case 3
0049         lower_bound = 0;
0050         upper_bound = Inf;
0051       case 2
0052         lower_bound = 0;
0053         upper_bound = Inf;
0054       case 4
0055         lower_bound = 0;
0056         upper_bound = 1;
0057       otherwise
0058         error('Unknown distribution!')
0059     end
0060 end
0061 if nargin==4
0062     switch distribution
0063       case 1
0064         upper_bound = Inf;
0065       case 3
0066         upper_bound = Inf;
0067       case 2
0068         upper_bound = Inf;
0069       case 4
0070         upper_bound = 1;
0071       otherwise
0072         error('Unknown distribution!')
0073     end
0074 end
0075 
0076 
0077 if (distribution==1)% Gamma distribution
0078     if m<lower_bound
0079         error('The mode has to be greater than the lower bound!')
0080     end
0081     if (m-lower_bound)<1e-12
0082         error('The gamma distribution should be specified with the mean and variance.')
0083     end        
0084     m = m - lower_bound ;
0085     beta  = -.5*m*(1-sqrt(1+4*s2/(m*m))) ;
0086     alpha = (m+beta)/beta ;
0087     parameters(1) = alpha;
0088     parameters(2) = beta;
0089     mu = alpha*beta + lower_bound ;
0090     return
0091 end
0092 
0093 if (distribution==2)% Inverse Gamma - 2 distribution
0094     if m<lower_bound+2*eps
0095         error('The mode has to be greater than the lower bound!')
0096     end
0097     m = m - lower_bound ;
0098     if isinf(s2)
0099         nu = 4;
0100         s  = 2*m;
0101     else
0102         delta = 2*(m*m/s2);
0103         poly = [ 1 , -(8+delta) , 20-4*delta , -(16+4*delta) ];
0104         all_roots  = roots(poly);
0105         real_roots = all_roots(find(abs(imag(all_roots))<2*eps));
0106         nu = real_roots(find(real_roots>2));
0107         s  = m*(nu+2);
0108     end
0109     parameters(1) = nu;
0110     parameters(2) = s;
0111     mu = s/(nu-2) + lower_bound;
0112     return
0113 end
0114 
0115 if (distribution==3)% Inverted Gamma 1 distribution
0116     if m<lower_bound+2*eps
0117         error('The mode has to be greater than the lower bound!')
0118     end
0119     m = m - lower_bound ;
0120     if isinf(s2)
0121         nu = 2;
0122         s  = 3*(m*m);
0123     else
0124         [mu, parameters] = mode_and_variance_to_mean(m,s2,2);
0125         nu = sqrt(parameters(1));
0126         nu2 = 2*nu;
0127         nu1 = 2;
0128         err = s2/(m*m) - (nu+1)/(nu-2) + .5*(nu+1)*(gamma((nu-1)/2)/gamma(nu/2))^2;
0129         while abs(nu2-nu1) > 1e-12
0130             if err < 0
0131                 nu1 = nu;
0132                 if nu < nu2
0133                     nu = nu2;
0134                 else
0135                     nu = 2*nu;
0136                     nu2 = nu;
0137                 end
0138             else
0139                 nu2 = nu;
0140             end
0141             nu =  (nu1+nu2)/2;
0142             err = s2/(m*m) - (nu+1)/(nu-2) + .5*(nu+1)*(gamma((nu-1)/2)/gamma(nu/2))^2;
0143         end
0144         s = (nu+1)*(m*m) ;
0145     end
0146     parameters(1) = nu;
0147     parameters(2) = s;
0148     mu = sqrt(.5*s)*gamma(.5*(nu-1))/gamma(.5*nu) + lower_bound ;
0149     return
0150 end
0151 
0152 if (distribution==4)% Beta distribution
0153     if m<lower_bound
0154         error('The mode has to be greater than the lower bound!')
0155     end
0156     if m>upper_bound
0157         error('The mode has to be less than the upper bound!')
0158     end
0159     if (m-lower_bound)<1e-12
0160         error('The beta distribution should be specified with the mean and variance.')
0161     end
0162     if (upper_bound-m)<1e-12
0163         error('The beta distribution should be specified with the mean and variance.')
0164     end
0165     ll = upper_bound-lower_bound;
0166     m  = (m-lower_bound)/ll;
0167     s2 = s2/(ll*ll);
0168     delta = m^2/s2;
0169     poly = NaN(1,4);
0170     poly(1) = 1;
0171     poly(2) = 7*m-(1-m)*delta-3;
0172     poly(3) = 16*m^2-14*m+3-2*m*delta+delta;
0173     poly(4) = 12*m^3-16*m^2+7*m-1;
0174     all_roots = roots(poly);
0175     real_roots = all_roots(find(abs(imag(all_roots))<2*eps));
0176     idx = find(real_roots>1);
0177     if length(idx)>1
0178         error('Multiplicity of solutions for the beta distribution specification.')
0179     elseif isempty(idx)
0180         disp('No solution for the beta distribution specification. You should reduce the variance.')
0181         error();
0182     end
0183     alpha = real_roots(idx);
0184     beta = ((1-m)*alpha+2*m-1)/m;
0185     parameters(1) = alpha;
0186     parameters(2) = beta;
0187     mu = alpha/(alpha+beta)*(upper_bound-lower_bound)+lower_bound;
0188     return
0189 end