semip_gcc.c

计量工具箱

		cobs = 0;
		for(i=0;i<m;i++){			
			obsi = mobs[i];
			
			phi = 0.0;
			for(j=cobs;j<cobs+obsi;j++){
				phi = phi + ((z[j] - xb[j])/v[i]);  // form phi
			}		
			awi = 0.0;
			aw2i = 0.0;
			for(j=0;j<m;j++){						// form awi and aw2i
				aw2i = aw2i + W2[i][j]*a[j];		// (W2[i][i] = 0)	
				if(j != i){
					awi = awi + (W[i + j*m] + W[j + i*m])*a[j];
				}
			}
			bi = phi + (rho/sige)*awi - ((rho*rho)/sige)*aw2i;
			di = (1/sige) + ((rho*rho)/sige)*w1[i] + (obsi/v[i]);		
			a[i] = (bi/di) + sqrt(1/di)*snorm();	// Form a[i]           

			cobs = cobs + obsi;
		}
	} // End update of a
	


	// Compute tvec = del*a	
	cobs = 0;
	for(i=0; i<m; i++){
		obsi = mobs[i];
		for(j=cobs; j<cobs + obsi; j++){
			tvec[j] = a[i];
		}
		cobs = cobs + obsi;
	}

//UPDATE: sige

	matvec(Bp,m,m,a,Bpa);			//form Bp*a
	aBBa = 0.0;						//form aBBa = a'*Bp'*Bp*a
	for(i=0; i<m; i++){		
		aBBa = aBBa + Bpa[i]*Bpa[i];
	}
	
	chi = genchi(m + 2.0*nu);

	sige = (aBBa + 2.0*d0)/chi;

	
//UPDATE: v (and vv = del*v)	
	

	for(i=0; i<n; i++){
		e[i] = e0[i] - tvec[i];		//form e  = z - x*bhat - tvec
	}
	cobs = 0;
	for(i=0; i<m; i++){
		obsi = mobs[i];
		chi = genchi(rval + obsi);
		ee = 0.0;
		for(j=cobs; j<cobs + obsi; j++){		// form ee
			ee = ee + e[j]*e[j];
		}
		v[i] = (ee + rval)/chi;					// form v
		for(j=cobs; j<cobs + obsi; j++){
			vv[j] = v[i];						// form vv
		}
		cobs = cobs + obsi;
	}
	


//UPDATE: rval (if necessary)

	if (mm != 0.0){
    rval = gengam(mm,kk);
    }

//UPDATE: rho (using univariate integration)

	      matvec(Wmat,m,m,a,Wa); // form Wa vector

		  	c0 = 0.0;						// form a'*a
		  	c1 = 0.0;                       // form a'*Wa
		  	c2 = 0.0;                       // form (Wa)'*Wa
	        for(i=0; i<m; i++){		
		    c0 = c0 + a[i]*a[i];
		    c1 = c1 + a[i]*Wa[i];
		    c2 = c2 + Wa[i]*Wa[i];
	        }

	        rho = draw_rho(rdet,ldet,c0,c1,c2,sige,ngrid,m,k,rho);



// (4) Update Bp matrix using new rho

	 
	for(i=0; i<m; i++){		                 
        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];
              }
		}
	}

	 
//UPDATE: z 
if (zflag == 1){ // skip this for continuous y-values
	 // (1) Generate vector of means

	 cobs = 0;
	 for(i=0; i<m; i++){
		 obsi = mobs[i];
		 for(j=cobs; j<cobs + obsi; j++){
			 mn[j] = xb[j] + a[i];
		 }			 		 
		 cobs = cobs + obsi;		 
	 }
	 

	 // (2) Sample truncated normal for latent z values

	 normal_truncc(n,mn,vv,limit1,limit2,z);

	 // (3) Compute associated sample y vector: ys

	 for(i=0; i<n; i++){
		 if(z[i]<0){
			 ys[i] = 0.0;
		 }else{
			 ys[i] = 1.0;
		 }
	 }

} // end of if zflag == 1


// ===================
// Save sample draws 	
// ===================


	 
	 if (iter > nomit-1){
		 pdraw[iter - nomit] = rho;							//save rho-draws
		 sdraw[iter - nomit] = sige;					    //save sige-draws	 
		 if(mm != 0.0)
			rdraw[iter - nomit] = rval;						//save r-draws (if necessary)
	     for(j=0; j<k; j++)
			bdraw[(iter - nomit) + j*nsave] = bhat[j];		//save beta-draws
		 for(i=0; i<m; i++){
			vmean[i] = vmean[i] + v[i]/((double) (nsave));	//save mean v-draws
			adraw[(iter - nomit) + i*nsave] = a[i];			//save a-draws
			amean[i] = amean[i] + a[i]/((double) (nsave));	//save mean a-draws
		 }
	     for(i=0; i<n; i++){
			zmean[i] = zmean[i] + z[i]/((double) (nsave));	//save mean z-draws
			yhat[i] = yhat[i] + mn[i]/((double) (nsave));   //save mean y-values
		 }
	 }

	 } // End iteration loop
	 

// ===============================
// END SAMPLER
// ===============================


// Free up allocated vectors

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


	free_dvector(z,0);
	free_dvector(tvec,0);   
	free_dvector(e,0);   
	free_dvector(e0,0);
	free_dvector(xb,0);
	free_dvector(mn,0);
	free_dvector(ys,0);      
    free_dvector(limit1,0);			
	free_dvector(limit2,0);
	free_dvector(zmt,0);
	free_dvector(vv,0);
	free_dvector(bhat,0);
	free_dvector(b0,0);
	free_dvector(b0tmp,0);
	free_dvector(v,0);
	free_dvector(b1,0);
	free_dvector(b1tmp,0);
	free_dvector(Bpa,0);
	free_dvector(rdet,0);
	free_dvector(ldet,0);
	free_dvector(w1,0);
	free_dvector(Wa,0);


} // End of semip_gc



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

// START GATEWAY FUNCTION

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

void mexFunction( int nlhs, mxArray *plhs[],
                  int nrhs, const mxArray *prhs[] )
{
  
// Local variables for the gateway function

// semip_gc(bdraw,adraw,pdraw,sdraw,rdraw,vmean,amean,zmean,yhat,
// y,x,W,ndraw,nomit,nsave,n,k,m,mobs,a,nu,d0,rval,mm,kk,detval,ngrid,TI,TIc)
     
	// (The inputs are passed to the subroutine and the outputs
	//  are extracted from the subroutine and passed back to Matlab.)

  double *bdraw,*adraw,*pdraw,*sdraw,*rdraw,*vmean,*amean,*zmean,*yhat;
  double *y,*x,*W,*a,nu,d0,rval,mm,kk,*detval,ngrid,*TI,*TIc,*m_obs;
  int i,ndraw,nomit,nsave,n,k,m,*mobs;
  
  long *Seed1, *Seed2;
  int buflen, status, flag;
  char *phrase;
  static char phrase2[8];
  
    phrase2[0] = 'h';
    phrase2[1] = 'h';
    phrase2[2] = 'i';
    phrase2[3] = 't';
    phrase2[4] = 'h';
    phrase2[5] = 'e';
    phrase2[6] = 'r';
    phrase2[7] = 'e';



  /* Check for proper number of arguments. */

  if(nrhs == 20) {  
     flag = 0;
    phrase = phrase2;
  } else if (nrhs == 21){
    flag = 1;
  } else {
    mexErrMsgTxt("semip_gc: 20 or 21 inputs required.");
  }

 
  if(nlhs != 9) {
    mexErrMsgTxt("semip_gcc: 9 output arguments needed");
  }


    if (flag == 1) {

    // input must be a string
    if ( mxIsChar(prhs[20]) != 1)
      mexErrMsgTxt("semip_gc: seed must be a string.");
    // input must be a row vector
    if (mxGetM(prhs[20])!=1)
      mexErrMsgTxt("semip_gc: seed input must be a row vector.");

    // get the length of the input string
    buflen = (mxGetM(prhs[20]) * mxGetN(prhs[20])) + 1;

    // allocate memory for input string
    phrase = mxCalloc(buflen, sizeof(char));

    // copy the string data from prhs[0] into a C string input_ buf.
    // If the string array contains several rows, they are copied,
    // one column at a time, into one long string array.
    //
    status = mxGetString(prhs[20], phrase, buflen);
    if(status != 0)
      mexWarnMsgTxt("semip_gc: Not enough space. seed string truncated.");
    }
    
    // allow the user to set a seed or rely on clock-based seed
   if (flag == 0) {
    setSeedTimeCore(phrase);
   } else {
	phrtsd(phrase,&Seed1,&Seed2);
    setall(Seed1,Seed2);
   }


//  Parse input arguments

// semip_gcc(y,x,W,ndraw,nomit,nsave,n,k,m,mobs,a,nu,d0,rval,mm,kk,detval,ngrid,TI,TIc);

     y = mxGetPr(prhs[0]);
     x = mxGetPr(prhs[1]);
	 W = mxGetPr(prhs[2]);
	 ndraw = (int) mxGetScalar(prhs[3]);
	 nomit = (int) mxGetScalar(prhs[4]);
	 nsave = (int) mxGetScalar(prhs[5]);
	 n = (int) mxGetScalar(prhs[6]);
	 k = (int) mxGetScalar(prhs[7]);
	 m = (int) mxGetScalar(prhs[8]);
	 m_obs = mxGetPr(prhs[9]);
     a = mxGetPr(prhs[10]);
	 nu = mxGetScalar(prhs[11]);
	 d0 = mxGetScalar(prhs[12]);
	 rval = mxGetScalar(prhs[13]);
	 mm = mxGetScalar(prhs[14]);
	 kk = mxGetScalar(prhs[15]);
	 detval = mxGetPr(prhs[16]);
	 ngrid = (int) mxGetScalar(prhs[17]);	 
	 TI = mxGetPr(prhs[18]);
	 TIc = mxGetPr(prhs[19]);
	 

	// no need for error checking on inputs
	// since this was done in the matlab function

	  // Transfer m_obs into an integer vector

	 //(1) Make integer vector storage
	 
	 mobs = ivector(0,m-1);

	 //(2) Put m_obs into mobs

	for(i=0; i<m; i++){
		mobs[i] = (int) m_obs[i];
	}
			

    /* Create matrix mxArrays for the return arguments */
	 //[bdraw,adraw,pdraw,sdraw,rdraw,vmean,amean,zmean,yhat,cntr]

    plhs[0] = mxCreateDoubleMatrix(nsave,k, mxREAL);	// beta draws
	plhs[1] = mxCreateDoubleMatrix(nsave,m, mxREAL);	// a-draws
	plhs[2] = mxCreateDoubleMatrix(nsave,1, mxREAL);	// rho draws
	plhs[3] = mxCreateDoubleMatrix(nsave,1, mxREAL);	// sige draws
	plhs[4] = mxCreateDoubleMatrix(nsave,1, mxREAL);	// r-draws
	plhs[5] = mxCreateDoubleMatrix(m,1, mxREAL);		// v-means
	plhs[6] = mxCreateDoubleMatrix(m,1, mxREAL);		// a-means
	plhs[7] = mxCreateDoubleMatrix(n,1, mxREAL);		// z-means
	plhs[8] = mxCreateDoubleMatrix(n,1, mxREAL);		// yhat-means


	// Get pointers to the data in the mxArrays for the C-subroutine


    bdraw = mxGetPr(plhs[0]);
	adraw = mxGetPr(plhs[1]);
	pdraw = mxGetPr(plhs[2]);	
	sdraw = mxGetPr(plhs[3]);
	rdraw = mxGetPr(plhs[4]);
	vmean = mxGetPr(plhs[5]);
	amean = mxGetPr(plhs[6]);
	zmean = mxGetPr(plhs[7]);
	yhat  = mxGetPr(plhs[8]);

// Call the Subroutine

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

// free up mobs
	
	free_ivector(mobs,0);


}