Exercise 38 solutions

% Question 1
%
x = [3 5 4 3 6 7 -1 2 -3 -4 2 -5 3 2 -1];
y = [2 7 4 7 6 9 -1 3 -2 -1 3 -1 2 4 2];
nboot = 1e6;
d_boot = ones(1,nboot) * NaN;
nx = length(x);
ny = length(y);
d_obs = mean(y)-mean(x);
xy = [x,y];
nxy = length(xy);

tic
for i=1:nboot
        % generate new x ix and new y iy
        ixy = randi(nxy,1,nxy);
        ix = xy(ixy(1:nx));
        iy = xy(ixy(nx+1:end));
        % recompute d statistic
        d_boot(i) = mean(iy)-mean(ix);
end
toc
p = length(find(d_boot>=d_obs)) / nboot;
disp(['p=',num2str(round(p*1000)/1000)])

% Question 2
%
matlabpool open
x = [3 5 4 3 6 7 -1 2 -3 -4 2 -5 3 2 -1];
y = [2 7 4 7 6 9 -1 3 -2 -1 3 -1 2 4 2];
nboot = 1e6;
d_boot = ones(1,nboot) * NaN;
nx = length(x);
ny = length(y);
d_obs = mean(y)-mean(x);
xy = [x,y];
nxy = length(xy);

tic
parfor i=1:nboot
        % generate new x ix and new y iy
        ixy = randi(nxy,1,nxy);
        ix = xy(ixy(1:nx));
        iy = xy(ixy(nx+1:end));
        % recompute d statistic
        d_boot(i) = mean(iy)-mean(ix);
end
toc
p = length(find(d_boot>=d_obs)) / nboot;
disp(['p=',num2str(round(p*1000)/1000)])
matlabpool close

// compile with: gcc -o e38_q1 e38_q1.c

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

// helper functions (see below)
//
double mean(double x[], int n);
void randperm(int perm[], int n);

int main(int argc, char *argv[]) {

  // initialize random seed based on current time
  srand((unsigned) time(NULL));

  // Question 1

  double x[] = {3,5,4,3,6,7,-1,2,-3,-4,2,-5,3,2,-1};
  double y[] = {2,7,4,7,6,9,-1,3,-2,-1,3,-1,2,4,2};
  double d = mean(y,15) - mean(x,15);
  int bootcount = 0;
  int i, j;
  int nboot = 1e6;
  double xy[30];
  double di;
  for (i=0; i<30; i++) {
    if (i<15) xy[i] = x[i];
    else xy[i] = y[i-15];
  }
  int seq[30];
  double xi[15], yi[15];

  for (i=0; i<nboot; i++) {
    randperm(seq, 30);
    for (j=0; j<30; j++) {
      if (j<15) xi[j] = xy[seq[j]];
      else yi[j-15] = xy[seq[j]];
    }
    di = mean(yi,15) - mean(xi,15);
    if (di >= d) bootcount = bootcount + 1;
  }
  double pboot = (double)bootcount / (double)nboot;
  printf("pboot = %7.5f\n", pboot);

  return 0;
}

/* ---------- HELPER FUNCTIONS ---------- */

double mean(double x[], int n) {
  int i;
  double m = 0.0;
  for (i=0; i<n; i++) {
    m = m + x[i];
  }
  m = m /n;
  return m;
}

void randperm(int perm[], int n) {
  int i, j, t;
  for(i=0; i<n; i++)
    perm[i] = i;
  for(i=0; i<n; i++) {
    j = rand()%(n-i)+i;
    t = perm[j];
    perm[j] = perm[i];
    perm[i] = t;
  }
}

// compile with: gcc -o e38_q2 e38_q2.c -fopenmp

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <omp.h>

// helper functions (see below)
//
double mean(double x[], int n);
void randperm(int perm[], int n);

int main(int argc, char *argv[]) {

  // initialize random seed based on current time
  srand((unsigned) time(NULL));

  // Question 2

  double x[] = {3,5,4,3,6,7,-1,2,-3,-4,2,-5,3,2,-1};
  double y[] = {2,7,4,7,6,9,-1,3,-2,-1,3,-1,2,4,2};
  double d = mean(y,15) - mean(x,15);
  int bootcount = 0;
  int i, j;
  int nboot = 1e6;
  double xy[30];
  double di;
  for (i=0; i<30; i++) {
    if (i<15) xy[i] = x[i];
    else xy[i] = y[i-15];
  }
  int seq[30];
  double xi[15], yi[15];

  #pragma omp parallel for private(i,j,xi,yi,seq)
  for (i=0; i<nboot; i++) {
    randperm(seq, 30);
    for (j=0; j<30; j++) {
      if (j<15) xi[j] = xy[seq[j]];
      else yi[j-15] = xy[seq[j]];
    }
    di = mean(yi,15) - mean(xi,15);
    if (di >= d) bootcount = bootcount + 1;
  }
  double pboot = (double)bootcount / (double)nboot;
  printf("pboot = %7.5f\n", pboot);

  return 0;
}

/* ---------- HELPER FUNCTIONS ---------- */

double mean(double x[], int n) {
  int i;
  double m = 0.0;
  for (i=0; i<n; i++) {
    m = m + x[i];
  }
  m = m /n;
  return m;
}

void randperm(int perm[], int n) {
  int i, j, t;
  for(i=0; i<n; i++)
    perm[i] = i;
  for(i=0; i<n; i++) {
    j = rand()%(n-i)+i;
    t = perm[j];
    perm[j] = perm[i];
    perm[i] = t;
  }
}