Home > matlab > draw_prior_density.m

draw_prior_density

PURPOSE ^

Computes values of prior densities at many points (before plotting)

SYNOPSIS ^

function [x,f,abscissa,dens,binf,bsup] = draw_prior_density(indx,bayestopt_);

DESCRIPTION ^

 Computes values of prior densities at many points (before plotting)

 INPUTS
    indx          [integer]    Parameter number.
    bayestopt_    [structure]  Describes the prior beliefs.
    
 OUTPUTS
    x             [double]     Row vector, subset of 'abscissa' such as the density is less than 10
    f             [double]     Row vector, subset of 'dens' such as the density is less than 10
    abscissa      [double]     Row vector, abscissa 
    dens          [double]     Row vector, density
    binf:         [double]     Scalar, first element of x
    bsup:         [double]     Scalar, last element of x

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [x,f,abscissa,dens,binf,bsup] = draw_prior_density(indx,bayestopt_);
0002 % Computes values of prior densities at many points (before plotting)
0003 %
0004 % INPUTS
0005 %    indx          [integer]    Parameter number.
0006 %    bayestopt_    [structure]  Describes the prior beliefs.
0007 %
0008 % OUTPUTS
0009 %    x             [double]     Row vector, subset of 'abscissa' such as the density is less than 10
0010 %    f             [double]     Row vector, subset of 'dens' such as the density is less than 10
0011 %    abscissa      [double]     Row vector, abscissa
0012 %    dens          [double]     Row vector, density
0013 %    binf:         [double]     Scalar, first element of x
0014 %    bsup:         [double]     Scalar, last element of x
0015 
0016 
0017 % Copyright (C) 2004-2011 Dynare Team
0018 %
0019 % This file is part of Dynare.
0020 %
0021 % Dynare is free software: you can redistribute it and/or modify
0022 % it under the terms of the GNU General Public License as published by
0023 % the Free Software Foundation, either version 3 of the License, or
0024 % (at your option) any later version.
0025 %
0026 % Dynare is distributed in the hope that it will be useful,
0027 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0028 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0029 % GNU General Public License for more details.
0030 %
0031 % You should have received a copy of the GNU General Public License
0032 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0033 
0034 pshape  = bayestopt_.pshape;
0035 p3      = bayestopt_.p3;
0036 p4      = bayestopt_.p4;
0037 p6      = bayestopt_.p6;
0038 p7      = bayestopt_.p7;
0039 
0040 truncprior = 1e-3;
0041 steps = 200;
0042 
0043 switch pshape(indx)
0044   case 1% Beta prior
0045     density = @(x,a,b,aa,bb) betapdf((x-aa)/(bb-aa), a, b)/(bb-aa);
0046     infbound = betainv(truncprior,p6(indx),p7(indx))*(p4(indx)-p3(indx))+p3(indx);
0047     supbound = betainv(1-truncprior,p6(indx),p7(indx))*(p4(indx)-p3(indx))+p3(indx);
0048     abscissa = linspace(infbound,supbound,steps);
0049     dens = density(abscissa,p6(indx),p7(indx),p3(indx),p4(indx));
0050   case 2% Generalized Gamma prior
0051     density = @(x,a,b,c) gampdf(x-c,a,b);
0052     try
0053         infbound = gaminv(truncprior,p6(indx),p7(indx))+p3(indx);
0054         supbound = gaminv(1-truncprior,p6(indx),p7(indx))+p3(indx);
0055     catch
0056         % Workaround for ticket #161
0057         if exist('OCTAVE_VERSION')
0058             error(['Due to a bug in Octave, you must choose other values for mean and/or variance of your prior on ' bayestopt_.name{indx} ', or use another shape'])
0059         else
0060             rethrow(lasterror)
0061         end
0062     end
0063     abscissa = linspace(infbound,supbound,steps);
0064     dens = density(abscissa,p6(indx),p7(indx),p3(indx));
0065   case 3% Gaussian prior
0066     infbound = norminv(truncprior,p6(indx),p7(indx)); 
0067     supbound = norminv(1-truncprior,p6(indx),p7(indx));
0068     abscissa = linspace(infbound,supbound,steps);
0069     dens = normpdf(abscissa,p6(indx),p7(indx));  
0070   case 4% Inverse-gamma of type 1 prior
0071     try
0072         infbound = 1/sqrt(gaminv(1-10*truncprior, p7(indx)/2, 2/p6(indx)))+p3(indx);
0073         supbound = 1/sqrt(gaminv(10*truncprior, p7(indx)/2, 2/p6(indx)))+p3(indx);
0074     catch
0075         % Workaround for ticket #161
0076         if exist('OCTAVE_VERSION')
0077             error(['Due to a bug in Octave, you must choose other values for mean and/or variance of your prior on ' bayestopt_.name{indx} ', or use another shape'])
0078         else
0079             rethrow(lasterror)
0080         end
0081     end
0082     abscissa = linspace(infbound,supbound,steps);
0083     dens = exp(lpdfig1(abscissa-p3(indx),p6(indx),p7(indx)));  
0084   case 5% Uniform prior
0085     infbound = p6(indx);
0086     supbound = p7(indx);
0087     abscissa = linspace(infbound,supbound,steps);
0088     dens = ones(1, steps) / (supbound-infbound);
0089   case 6% Inverse-gamma of type 2 prior
0090     try
0091         infbound = 1/(gaminv(1-10*truncprior, p7(indx)/2, 2/p6(indx)))+p3(indx);
0092         supbound = 1/(gaminv(10*truncprior, p7(indx)/2, 2/p6(indx)))+p3(indx);
0093     catch
0094         % Workaround for ticket #161
0095         if exist('OCTAVE_VERSION')
0096             error(['Due to a bug in Octave, you must choose other values for mean and/or variance of your prior on ' bayestopt_.name{indx} ', or use another shape'])
0097         else
0098             rethrow(lasterror)
0099         end
0100     end
0101     abscissa = linspace(infbound,supbound,steps);
0102     dens = exp(lpdfig2(abscissa-p3(indx),p6(indx),p7(indx)));
0103   otherwise
0104     error(sprintf('draw_prior_density: unknown distribution shape (index %d, type %d)', indx, pshape(indx)));
0105 end 
0106 
0107 if pshape(indx) ~= 5 
0108     [junk,k1] = max(dens);
0109     if k1 == 1 || k1 == length(dens)
0110         k = find(dens > 10);
0111         dens(k) = NaN;
0112     end
0113 end
0114 binf = abscissa(1);
0115 bsup = abscissa(end);
0116 x = abscissa;
0117 f = dens;
0118 f(find(x<bayestopt_.lb(indx)))=0;
0119 f(find(x>bayestopt_.ub(indx)))=0;

Generated on Mon 21-May-2012 02:42:43 by m2html © 2005