Needs persistent p1 and p2;

global estim_params abd bayestopt_.
Debugged generalised beta prior.
index (i) missing somewhere


git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@1104 ac1d8469-bf42-47a9-8791-bf33cf982152

matlab/prior_draw.m

@@ -19,8 +19,8 @@ function pdraw = prior_draw(init,cc)
 %  
 % part of DYNARE, copyright S. Adjemian, M. Juillard (2006)
 % Gnu Public License.
-global M_ options_ estim_params_ bayestopt_ 
-persistent fname npar bounds pshape pmean pstd a b p3 p4 condition
+global estim_params_  bayestopt_
+persistent fname npar bounds pshape pmean pstd a b p1 p2 p3 p4 condition
   
 if init
     nvx  = estim_params_.nvx;
@@ -29,8 +29,6 @@ if init
     ncn  = estim_params_.ncn;
     np   = estim_params_.np ;
     npar = nvx+nvn+ncx+ncn+np;
-%    MhDirectoryName = CheckPath('metropolis');
-%    fname = [ MhDirectoryName '/' M_.fname];
     pshape = bayestopt_.pshape;
     pmean = bayestopt_.pmean;
     pstd  = bayestopt_.pstdev;
@@ -43,22 +41,22 @@ if init
     if nargin == 2
         bounds = cc;
     else
-        bounds = [-Inf Inf];
+        bounds = kron(ones(npar,1),[-Inf Inf]);
     end
     for i = 1:npar
         switch pshape(i)
           case 3% Gaussian prior
-%            b(i) = pstd(i)^2/(pmean(i)-p3(i));
-%            a(i) = (pmean(i)-p3(i))/b(i);
+            b(i) = pstd(i)^2/(pmean(i)-p3(i));
+            a(i) = (pmean(i)-p3(i))/b(i);
           case 1% Beta prior
             mu = (p1(i)-p3(i))/(p4(i)-p3(i));
             stdd = p2(i)/(p4(i)-p3(i));
             a(i) = (1-mu)*mu^2/stdd^2 - mu;
-	  b(i) = a*(1/mu - 1);        
+            b(i) = a(i)*(1/mu - 1);        
           case 2;%Gamma prior
             mu = p1(i)-p3(i);
             b(i) = p2(i)^2/mu;
-	  a(i) = mu/b;
+            a(i) = mu/b(i);
           case {5,4,6}
             % Nothing to do here
             %
@@ -109,7 +107,8 @@ for i = 1:npar
             y2 = gamma_draw(b(i),1,0);
             tmp = y1/(y1+y2);
             if tmp >= bounds(i,1) && tmp <= bounds(i,2)
-                pdraw(i) = pmean(i)+tmp*pstd(i);
+                %pdraw(i) = pmean(i)+tmp*pstd(i);
+                pdraw(i) = p3(i)+tmp*(p4(i)-p3(i));
                 break
             end
         end
@@ -138,13 +137,13 @@ end
 
   
   
-function gamma_draw(a,b,c)
+function g = gamma_draw(a,b,c)
 % Bauwens, Lubrano & Richard (page 316)
 if a >30
     z = randn;
     g = b*(z+sqrt(4*a-1))^2/4 + c; 
 else
-    condi = 1
+    condi = 1;
     while condi
         x = -1;
         while x<0