library(ggplot2)
set.seed(1) #set the seed so results are replicable
start = proc.time()[1] #start a clock so later we can determine how long this took

#define the ex-Gaussian parameter space to explore
mu = rep(c(300,350,400,450,500,550),each=2)
sigma = rep(c(20,50),times=6)
tau = rep(c(300,250,200,150,100,50),each=2)

#Obtain the expected values for the median via brute force data generation
a = data.frame(mu,sigma,tau)
a = ddply(
	.data = a
	,.variables = .(mu,sigma,tau)
	, function(x){
		to_return =  median(rnorm(1e7,x$mu,x$sigma)+rexp(1e7,1/x$tau))
		names(to_return) = 'expected_median'
		return(to_return)
	}
	,.progress = 'text'
)

#define skew as the difference between the mean (mu+tau) and the expected median
a$skew = a$mu+a$tau-a$expected_median

#add k (sample size) to the parameter space data frame
k = c(4,6,8,10,15,20,25,35,50,100)
temp = NULL
for(this_k in k){
	temp = rbind(temp,cbind(a,k=this_k))
}
a=temp
rm(temp)

#run the main simulation
a=ddply(
	.data = a
	, .variables = .(skew,k)
	, function(x){
		rt = matrix(rnorm(1e5*x$k,x$mu,x$sigma)+rexp(1e5*x$k,1/x$tau),1e5,x$k)
		obs_meds = apply(rt,1,median)
		bias = mean(obs_meds) - x$expected_median
		ev = sd(obs_meds)
		to_return = c(bias,ev)
		names(to_return) = c('bias','ev')
		return(to_return)
	}
	,.progress = 'text'
)
proc.time()[1]-start #show how long this took

#get ready for some plotting
a$bias_printable = as.character(round(a$bias))
a$ev_printable = as.character(round(a$ev))

#plot the bias
bias_plot = ggplot()
bias_plot = bias_plot + layer(
	geom = 'tile'
	, data = a
	, mapping = aes(
		x = k
		, y = skew
		, fill = bias
	)
)
bias_plot = bias_plot + layer(
	geom = 'text'
	, data = a
	, mapping = aes(
		x = k
		, y = skew
		, label = bias_printable
	)
	, size = 1
	, colour = 'white'
	, angle = 90
)
ggsave(bias_plot , file = 'median_bias.pdf' , width = 4 , height=  4)

#plot the ev
ev_plot = ggplot()
ev_plot = ev_plot + layer(
	geom = 'tile'
	, data = a
	, mapping = aes(
		x = k
		, y = skew
		, fill = ev
	)
)
ev_plot = ev_plot + layer(
	geom = 'text'
	, data = a
	, mapping = aes(
		x = k
		, y = skew
		, label = ev_printable
	)
	, size = 1
	, colour = 'white'
	, angle = 90
)
ggsave(ev_plot , file = 'median_ev.pdf' , width = 4 , height=  4)