/* DIRICHLET.GPR: program to simulate draws from a Dirichlet Distribution
		  John Rust, Yale University, April, 1997 */


   n=1000; /* set number of simulation draws true DGP to construct 
	      posterior distribution in setparm.g */

   nsim=1000; /* set number of simulation draws for computing the
		 probability posterior density is with a ball of
		 radius tol about true theta using simple frequency
		 simulator */

   tol=.01;  /* tolerance for distance between simulated draws
		and true parameter theta */

   {a,cnts,theta}=setparm(n);  /* compute Dirichlet posterior where vector "a"
	         	          defines the hyperparameters of the Dirichlet
				  prior, cnts are the counts of the
				  observed data in each cell, so a+cnts
				  are the parameters of the Dirichlet
				  posterior */

   k=rows(a);              /* dimension of Dirichlet density */

   gam=zeros(k,1);
   dirich=gam;
   distance=zeros(nsim,1);
   distance1=zeros(nsim,1);
   ap=a+cnts;

   j=1; do until j > nsim;   /* simulation loop */

   i=1; do until i > k;

			  /* simulate Gammas as a sum of exponentials */

      gam[i]=sumc(-ln(rndu(ap[i],1)));
      dirich[i]=gam[i];

   i=i+1; endo;
                          /* simulate Dirichlet as normalized Gammas */
   dirich=dirich/sumc(gam);  

                          /* count up whether simulation is within preset
			     tolerance, tol, of theta */
  
   distance[j]=sqrt(sumc((dirich-theta)^2));

   j=j+1; endo;

   "Example of  posterior calculations with Dirichlet conjugate prior";
   "";
   "Prior is Dirichlet with parameters "; a';
   "Posterior is Dirichlet with parameters "; ap';
   "";
   "Simulated posterior probability of being within";
   "a ball of radius "$+ftos(tol,"*.*lf",3,2)$+" about theta* ";;
   sumc((distance .< tol))/nsim;  /* report fraction of simulations within
				     tolerance of theta */
 
   "Normal approximation of posterior probability of being within";
   "a ball of radius "$+ftos(tol,"*.*lf",3,2)$+" about theta* ";;

   /* posterior is approximately normal with mean equal to the MLE of
      theta, cnts/n, with cov matrix equal to inverse of information
      matrix */

   mle=cnts/n;  
   sigma=zeros(3,3);

   sigma[1,1]=mle[1]*(1-mle[1]);
   sigma[2,2]=mle[2]*(1-mle[2]);
   sigma[3,3]=mle[3]*(1-mle[3]);
   sigma[1,2]=-mle[1]*mle[2];
   sigma[2,1]=sigma[1,2];
   sigma[1,3]=-mle[1]*mle[3];
   sigma[3,1]=sigma[1,3];
   sigma[2,3]=-mle[2]*mle[3];
   sigma[3,2]=sigma[2,3];

   sigma=chol(sigma/n);

   j=1; do until j > nsim;

       simthet=mle[1:k-1]+sigma*rndn(k-1,1);
       simthet=simthet|(1-sumc(simthet));
       distance1[j]=sqrt(sumc((simthet-theta)^2));

   j=j+1; endo; 

   sumc((distance1 .< tol))/nsim;  /* report fraction of simulations within
				     tolerance of theta */