In [1]:
require('ggplot2')
require('reshape')


Loading required package: ggplot2
Loading required package: reshape

In [2]:
trials <- 100

In [3]:
epsilons <- seq(0, 0.2, 0.01)
epsilons


Out[3]:
  1. 0
  2. 0.01
  3. 0.02
  4. 0.03
  5. 0.04
  6. 0.05
  7. 0.06
  8. 0.07
  9. 0.08
  10. 0.09
  11. 0.1
  12. 0.11
  13. 0.12
  14. 0.13
  15. 0.14
  16. 0.15
  17. 0.16
  18. 0.17
  19. 0.18
  20. 0.19
  21. 0.2

In [4]:
rcontnorm <- function(n=1, mean=0, sd=1, eps=0.1, sdc=10) {
    # Generate random numbers distributed according to a contaminated normal distribution.
    #
    # Arguments:
    #        n: number of observations.
    #     mean: vector of means.
    #       sd: vector of standard deviations for base.
    #      eps: factor of contamination.
    #      sdc: vector of standard deviations for contamination.
    t_choice <- rbinom(n, 1, 1-eps)
    t_norm <- rnorm(n)
    t_contam <- rnorm(n, sd=sdc)
    return (t_choice * t_norm + (1-t_choice) * t_contam)
}

In [15]:
df <- data.frame(trial=1:trials, sapply(epsilons, rcontnorm, n=trials, mean=0, sd=1, sdc=10))
colnames(df)[2:(length(epsilons)+1)] <- epsilons
melt_df <- melt(df, id.vars='trial', variable_name='eps')

In [16]:
ggplot(melt_df, aes(trial, value, color=eps)) + geom_point()



In [20]:
colMeans(df)


Out[20]:
trial
50.5
0
0.221357298830456
0.01
-0.132603125023709
0.02
-0.261128044449697
0.03
0.0478295471543567
0.04
0.139574306292092
0.05
0.103918395034697
0.06
0.226033750113427
0.07
-0.187696232658634
0.08
0.049426890917177
0.09
-0.56710234784771
0.1
0.183406088957881
0.11
0.246545905589974
0.12
0.110739121536215
0.13
-0.112260834657708
0.14
0.198349889213403
0.15
0.165620190075764
0.16
0.186791335660266
0.17
0.210885571603444
0.18
-0.419974027925565
0.19
0.164310283487101
0.2
-0.00697696986081378

In [ ]: