#include 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