# simulate N trials of two die four sided die tosses
N = 10000;

print("dice 1 tosses");
d1s = sample(1:4,N,replace=TRUE)

print("dice 2 tosses");
d2s = sample(1:4,N,replace=TRUE)

scores <- list();
for(tn in 1:N) {
	scores[[tn]] = max(d1s[tn],d2s[tn]);
}

mean_score = mean(unlist(scores))

print("sample mean value:");
print(mean_score)

print("here is a histogram of the scores");
hist(unlist(scores));

# note for small N, mean_score is not necessarily close to expected value 
# of probability distribution of this problem. 
# for large N, it gets very close to this value.