semip_gcc.c

计量工具箱

#include <stdio.h>
#include <math.h>
#include <string.h>
#include <malloc.h>
#include <mex.h>
#include "randomlib.h"
#include "matrixjpl.h"

// *****************************************************
// mex semip_gcc file (semip_gcc.c)

// *****************************************************

// routine to draw rho from the conditional posterior distribution
double draw_rho(double *rvec,
				double *ldet,
				double epe0,
				double eped,
				double epe0d,
				double sige,
				int nrho,
				int n,
				int k,
     			double rho)
{

	int i;
	double adj, nmk, nmk2, rnd;
	double *s, *den, *num;
	double dsum, rsum, z;
	

// compute denominator (integrating constant)
	      s    = dvector(0,nrho-1);
	      den  = dvector(0,nrho-1);
	      num  = dvector(0,nrho-1);

		  nmk = (double)(n-k);
		  nmk2 = nmk/2;	  

	  for(i=0; i<nrho; i++){
		      z = epe0 - 2.0*rvec[i]*eped + rvec[i]*rvec[i]*epe0d;
              s[i] = ldet[i] - (z/2*sige);
	  }

// adjustment for scaling that subtracts the maximum value
	adj = s[0];
	for(i=1; i<nrho; i++){
		if (s[i] > adj)
			adj = s[i];
	}
	
	for(i=0; i<nrho; i++){
		den[i] = (s[i] - adj);
		den[i] = exp(den[i]);
		}

	dsum = 0.0; // trapezoid rule
    for(i=0; i<nrho-1; i++){
		dsum = dsum + (rvec[i+1]+rvec[i])*(den[i+1]-den[i])/2;
	}

	    rsum = 0.0;
		for(i=0; i<nrho; i++){
		s[i] = dabs(den[i]/dsum); // normalize rho post
        rsum = rsum + s[i];
		if (i == 0){
		den[i] = s[i];
		} else {
		den[i] = den[i-1] + s[i]; // cumulative sum
		}
	}
	
    rnd = rsum*ranf();

	// create rho draw via inversion
	for(i=0; i<nrho; i++){
	  if (rnd <= den[i]){
	  rho = rvec[i];
	  break;
	  }
	  }
	  	  
	free_dvector(s,0);
	free_dvector(den,0);
	free_dvector(num,0);

	return (rho);

}




// *****************************************************
// MAIN SAMPLER
// *****************************************************

// Estimates robust spatial probit model using MCMC

// semip_gc(bdraw,adraw,pdraw,sdraw,rdraw,vmean,amean,zmean,yhat,cntr,
// y,x,W,ndraw,nomit,nsave,n,k,m,mobs,a,nu,d0,rval,mm,kk,detval,ngrid,TI,TIc)
     

void semip_gc(

		//Output Arguments

        double *bdraw, // draws for beta (ndraw - nomit,k) matrix
		double *adraw, // draws for regional effects, a, (m,1) vector
        double *pdraw, // draws for rho (ndraw - nomit,1) vector
		double *sdraw, // draws for sige (ndraw - nomit,1) vector
		double *rdraw, // draws for rval (ndraw - nomit,1) vector (if mm != 0)
		double *vmean, // mean of vi draws (n,1) vector			
        double *amean, // mean of a-draws (m,1) vector
        double *zmean, // mean of latent z-draws (n,1) vector
        double *yhat,  // mean of posterior predicted y (n,1) vector
            
		//Input Arguments

		double *y,		// (n,1) lhs vector with 0,1 values
 		double *x,		// (n,k) matrix of explanatory variables   
 		double *W,		// (m,m) weight matrix
        int ndraw,		// # of draws
        int nomit,		// # of burn-in draws to omit
		int nsave,		// # of draws saved (= ndraw - nomit)
        int n,			// # of observations
		int k,			// # of explanatory variables
        int m,			// # of regions
        int *mobs,		// (m,1) vector of obs numbers in each region
		double *a,		// (m,1) vector of regional effects
		double nu,		// prior parameter for sige	 
        double d0,		// prior parameter for sige
        double rval,	// hyperparameter r            
        double mm,		// prior parameter for rval
        double kk,		// prior parameter for rval
        double *detval, // (ngrid,2) matrix with [rho , log det values]            
        int ngrid,		// # of values in detval (rows)
        double *TI,		// prior var-cov for beta (inverted in matlab)
		double *TIc)	// prior var-cov * prior mean	

{
    
// Local Variables

	int i,j,l,iter,invt,accept,rcount,cobs,obsi,zflag;
    double rmin,rmax,sige,chi,ee,ratio,p,rho,rho_new;
	double junk,aBBa,vsqrt,ru,rhox,rhoy,phi,awi,aw2i,bi,di;
    double *z,*tvec,*e,*e0,*xb,*mn,*ys,*limit1,*limit2,*zmt;
	double *bhat,*b0,*b0tmp,*v,*vv,*b1,*b1tmp,*Bpa,*rdet,*ldet,*w1;
	double **xs,**xst,**xsxs,**A0,**A1,**Bpt,**Bp,**xmat,**W2; 	
	double *Wa, **Wmat;
	double c0, c1,c2;

		   
// Allocate Matrices and Vectors


	  xs     = dmatrix(0,n-1,0,k-1);	  
	  xst    = dmatrix(0,k-1,0,n-1);      
	  xsxs   = dmatrix(0,k-1,0,k-1);
	  A0     = dmatrix(0,k-1,0,k-1);	  
	  A1     = dmatrix(0,m-1,0,m-1);
	  Bp     = dmatrix(0,m-1,0,m-1);
	  Bpt    = dmatrix(0,m-1,0,m-1);
	  xmat	 = dmatrix(0,n-1,0,k-1);
	  W2     = dmatrix(0,m-1,0,m-1);
	  Wmat   = dmatrix(0,m-1,0,m-1);

	   z     = dvector(0,n-1);
	   tvec  = dvector(0,n-1);
	   e     = dvector(0,n-1);
	   e0    = dvector(0,n-1);
	   xb    = dvector(0,n-1);
	   mn    = dvector(0,n-1);
	   ys    = dvector(0,n-1);      
      limit1 = dvector(0,n-1);			
	  limit2 = dvector(0,n-1);
	  zmt	 = dvector(0,n-1);
	  vv	 = dvector(0,n-1);
	  bhat   = dvector(0,k-1);
	  b0     = dvector(0,k-1);
	  b0tmp  = dvector(0,k-1);
	  v      = dvector(0,m-1);
	  b1     = dvector(0,m-1);
	  b1tmp  = dvector(0,m-1);
	  Bpa    = dvector(0,m-1);
	  w1     = dvector(0,m-1);
	  rdet   = dvector(0,ngrid-1);
	  ldet   = dvector(0,ngrid-1);	 
	  Wa     = dvector(0,m-1);
   
// Initializations


    junk = 0.0;  // Placeholder for mean in meanvar()
  	rho = 0.5;  
	sige = 1.0;
	zflag = 0;  // a flag for 0,1 y-values

// Initialize (z,vv,limits,tvec) 

	for(i=0; i<n; i++){
        z[i] = y[i];
		vv[i] = 1.0;
		tvec[i] = 1.0;
	    if (y[i] == 0.0){
			limit1[i] = -10000;
			limit2[i] = 0.0;
			zflag = 1; // we have 0,1 y-values so sample z
	    }else{
			limit1[i] = 0.0;
			limit2[i] = 10000;
		}
	}
        
	
	// Initialize v, b1tmp, Bp


	for(i=0; i<m; i++){
		v[i] = 1.0;
	    b1tmp[i] = 1.0;		                 
        for(j=0; j<m; j++){
              if (j == i){
                 Bp[i][j] = 1 - rho*W[i + j*m];
			  }else{
                 Bp[i][j] = - rho*W[i + j*m];
              }
              Wmat[i][j] = W[i + j*m];
		}
	}


	
	// Parse detval into rdet and ldet vectors 
	//       and define (rmin,rmax)
	
	for(i=0; i<ngrid; i++){
	    j=0;
		rdet[i] = detval[i + j*ngrid];
	    j=1;
		ldet[i] = detval[i + j*ngrid];
	}
	
	rmin = rdet[0];
    rmax = rdet[ngrid-1];  

	
	// Put x into xmat

	for(i=0; i<n; i++){
		for(j=0; j<k; j++)
            xmat[i][j] = x[i + j*n];
	}


	// Compute matrices to be used for updating a-vector

	if(m > 100){  // Used only for large m

		for(i=0; i<m; i++){

			w1[i] = 0.0;							// Form w1[i]
			for(j=0;j<m;j++){
				w1[i] = w1[i] + (W[j + i*m]*W[j + i*m]);
			}

			for(j=0;j<m;j++){						// Form W2[i][*]
				W2[i][j] = 0.0;
				if(j != i){					
					for(l=0;l<m;l++){
						W2[i][j] = W2[i][j] + W[l + i*m]*W[l + j*m];
					}
				}		// Note that W2[i][i] = 0.0 by construction
			}			
		}
	} // End if(m > 10)



// ======================================
// Start the Sampler
// ======================================

	for(iter=0; iter<ndraw; iter++){


// UPDATE: beta

// (1) Apply stnd devs (vsqrt) to x, z, z - tvec
	
	for(i=0; i<n; i++){
		vsqrt = sqrt(vv[i]);
	    zmt[i] = (z[i] - tvec[i])/vsqrt;
		for(j=0; j<k; j++)
			xs[i][j] = x[i + j*n]/vsqrt;
	}
	      

//  (2) Construct A0 matrix


	transpose(xs,n,k,xst);				// form xs'
	  matmat(xst,k,n,xs,k,xsxs);		// form xs'*xs
      for(i=0; i<k; i++){				// form xs'*xs + TI
            for(j=0; j<k; j++)
                 A0[i][j] = xsxs[i][j] + TI[i + j*k];
          }      

      invt = inverse(A0, k);			// replace A0 = inv(A0)
	  if (invt != 1)
		 mexPrintf("semip_gc: Inversion error in beta conditional \n");

// (3) Construct b0 vector

	matvec(xst,k,n,zmt,b0tmp);			//form xs'*zmt
      for(i=0; i<k; i++){
		b0tmp[i] = b0tmp[i] + TIc[i];		//form b0tmp = xs'*zmt + TIc
      }
	matvec(A0,k,k,b0tmp,b0);				//form b0 = A0*b0tmp

// (4) Do multivariate normal draw from N(b0,A0)

	normal_rndc(A0, k, bhat);			// generate N(0,A0) vector
	for(i=0; i<k; i++){
		bhat[i] = bhat[i] + b0[i];		// add mean vector, b0
	}

// (5) Now update related values:


	matvec(xmat,n,k,bhat,xb);			// form xb = xmat*bhat


	for(i=0; i<n; i++){					// form e0 = z - x*bhat
		e0[i] = z[i] - xb[i];	
	}

	cobs = 0;
	for(i=0; i<m; i++){					// form b1tmp = e0i/v[i]
		obsi = mobs[i];
		b1tmp[i] = 0.0;
		for(j=cobs; j<cobs + obsi; j++){
			b1tmp[i] = b1tmp[i] + (e0[j]/v[i]);	
		}
		cobs = cobs + obsi;
	}
	
        

// UPDATE: a 
	
	if(m <= 100){  // For small m use straight inverse

		// (1) Define A1 and b1

		transpose(Bp,m,m,Bpt);			    // form Bp'
		matmat(Bpt,m,m,Bp,m,A1);			// form Bp'*Bp
	
		for(i=0; i<m; i++){					// form A1 = (1/sige)*Bp'Bp + diag(mobs/v)
			for(j=0; j<m; j++){
				if (j == i){
					A1[i][j] = (1/sige)*A1[i][j] + ((double)mobs[i])/v[i];
				}else{
					A1[i][j] = (1/sige)*A1[i][j];
				}
			}
		}
	

		inverse(A1,m);						// set A1 = inv(A1)
	
		matvec(A1,m,m,b1tmp,b1);			// form b1

	
		// (2) Do multivariate normal draw from N(b1,A1)

		normal_rndc(A1,m,a);				// generate N(0,A1) vector
		for(i=0; i<m; i++){
			a[i] = a[i] + b1[i];			// add mean vector, b1
		}

	}else{   // For large m use marginal distributions