#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
NumericVector Cfunc(NumericVector xvec, IntegerVector M, NumericVector beta0, NumericVector alpha) {

    int d = 0, m, i, n = xvec.size();
    NumericVector yvec(n);
    double meanxy = 0.0, meanx = 0.0, meany = 0.0, meanx2 = 0.0, meany2 = 0.0;
    double thresh, num = 0.0, denom = 0.0, tobs, beta1hat, beta0hat, sighat, sighatbeta1hat;
    thresh = R::qt(1.0 - alpha[0] / 2.0, (double)(n - 2), 1, 0);

    for (i = 0; i< n; i++) {
      meanx = meanx + xvec[i];
      meanx2 = meanx2 + R_pow(xvec[i], 2.0);
    }
    meanx = meanx / (double)n;
    meanx2 = meanx2 / (double)n;

    GetRNGstate();
    for (m = 0; m < M[0]; m++) {
      meany = 0;
      meany2 = 0;
      meanxy = 0;
      for (i = 0; i < n; i++) { 
	yvec[i] = beta0[0] + R::runif(0.0, 1.0);
	meany = meany + yvec[i];
	meany2 = meany2 + R_pow(yvec[i], 2.0);
	meanxy = meanxy + xvec[i] * yvec[i]; 
      }
      meany = meany / (double)n;
      meany2 = meany2 / (double)n;
      meanxy = meanxy / (double)n;
      
      num = meanxy - meanx * meany;
      denom = meanx2 - meanx * meanx;
      
      beta1hat = num / denom;
      beta0hat = meany - beta1hat * meanx;
      
      sighat = sqrt((double)n * (meany2 + beta0hat * beta0hat + beta1hat * beta1hat * meanx2 - 2.0 * beta0hat * meany
				 - 2.0 * beta1hat * meanxy + 2.0 * beta0hat * beta1hat * meanx) / (double)(n - 2));
      
      sighatbeta1hat = sighat / sqrt((double)n * denom);
      
      tobs = beta1hat / sighatbeta1hat;
      
      if (fabs(tobs) > thresh) d = d + 1;
    }

    PutRNGstate();
    return Rcpp::wrap((double)d / (double)M[0]);
    
  }	// End Cfunc

