set.seed(12345); #set the random number generator
x = rnorm(100);
y = rnorm(100,mean=0.5)

# t-test, although they are clearly not normal
t.test( x, y, alternative="two.sided"); 
t.test( x, y, alternative="less"); #one-sided test

tt = t.test( x, y, alternative="less");

################
#alternatively, let's try simulation 

t.stat = function( x1, x2 ) {
  ( mean(x1) - mean(x2)) / sqrt( var(x1)/length(x1) + var(x2)/length(x2));
}

t.obs = t.stat(x,y);

n=1000 #1K permutations
permu = data.frame( matrix(NA, nrow=n, ncol=3) );
names(permu) = c("xbar","ybar","t");
head(permu);

## Reset:
#RNGkind(ok[1])

i=1;
while( i <= n ) {
  new  = sample( c(x,y), replace = TRUE); 
  newX = new[1:100];
  newY = new[101:200];
  permu$t[i] = t.stat( newX, newY);
  permu$xbar = mean(newX);
  permu$ybar = mean(newY);
  i = i+1; 
}

sub = permu[permu$t< t.obs, ]
my.p = length(sub$t) / length(permu$t);
my.p

results = rbind(my.p, tt$p.value);
rownames(results) = c("permutation", "t-test");
colnames(results) = c("p");
results;

hist( permu$t, breaks=20, xlab="t-statistic", main="p-value and simulation" );
arrows( t.obs, n/10, t.obs, n/20, col="red");