Code for Monte Carlo simulation: generate samples of given size in R

Question

I started by generating a sample of 500 uniformly-distributed random numbers between 0 and 1 using the code below:

set.seed(1234)
X<-runif(500, min=0, max=1)

Now, I need to write a psuedocode that generates 10000 samples of N=500 for a MC simulation, compute the mean of my newly created X, and store the iteration number and mean value in a result object. I have never attempted this, and so far I have this:

n.iter <-(10000*500)
results <- matrix (0, n.iter, 4)

Finally, once this is accomplished, I'm to run it, then obtain median, mean, and min/max of the accrued sample means and save them to a data frame called MC.table. (Also note, above, I have no idea why there's a "4" in the matrix code --- I'm working off of previous examples). Any advice or help would be greatly appreciated.

EDIT: I have an example that may work, but I don't really understand what's going on with it, so please elaborate on its validity for this:

Ni <- 10000
n <- 500
c <- 0

for (i in n){
for (j in 1:Ni){
c <- c+ 1
d <- data.frame (x= , y= )
results [c,1] <- c
results [c,2] <- j
results [c,3] <- i
results [c,4] <- something( d$x, d$y)
rm (d) } }

If you could even take the time to explain what that means, that'd go a long way to helping me! Thanks!

Answer 1

You could try using data.table, a package that can be installed using install.packages("data.table"). With that installed, you would run something like...

> require(data.table)
> dt <- data.table(x=runif(500*10000),iter=rep(1:500,each=10000))
                  # x iter
      # 1: 0.48293196    1
      # 2: 0.61935416    1
      # 3: 0.99831614    1
      # 4: 0.26944687    1
      # 5: 0.38027524    1
     # ---                
# 4999996: 0.11314160  500
# 4999997: 0.07958396  500
# 4999998: 0.97690312  500
# 4999999: 0.81670765  500
# 5000000: 0.62934609  500
> summaries <- dt[,list(mean=mean(x),median=median(x)),by=iter]
     # iter      mean    median
  # 1:    1 0.5005310 0.5026592
  # 2:    2 0.4971551 0.4950034
  # 3:    3 0.4977677 0.4985360
  # 4:    4 0.5034727 0.5052344
  # 5:    5 0.4999848 0.4971214
 # ---                         
# 496:  496 0.5013314 0.5048186
# 497:  497 0.4955447 0.4941715
# 498:  498 0.4983971 0.4910115
# 499:  499 0.5000382 0.4997024
# 500:  500 0.5009614 0.4988237
> min_o_means <- min(summaries$mean)
# [1] 0.4914826

I think the syntax is fairly straightforward. You may want to look up some of the functions using ? (e.g., ?rep). The lines starting with # are just displaying the generated objects. In data.tables, the number to the left of the : is just the row number and --- indicates rows that are skipped in the display.

Answer 2

I guess the answer I would give would really depend on if you want to learn to pseudocode or if you want to learn the "R" ish way of doing it. This answer is what I would recommend for somebody who wanted to learn how to work with R.

First I would make a matrix with N columns and 10000 rows. R appreciates it when we make the space ahead of time for the numbers to go into.

X=matrix(NA,nrow=10000,ncol=500)

You know how to generate the 500 random variables for one row.

runif(500,0,1)

Now you need to figure out how to get that to happen 10000 times and assign each one to X. Maybe a for loop would work.

for(i in 1:10000) X[i,]=runif(500,0,1)

Then you need to figure out how to get the summaries of each row. One function that may help is rowMeans(). Look at its help page and then try to get the mean of each row of your table X

to get the means of each iteration

rowMeans(X)

then to get an idea of what these numbers are like I might be inclined to

plot(rowMeans(X))

Answer 3

I think you are describing a simple bootstrap. Eventually, you might want to use the function boot. But until you understand the mechanics, I feel that loops are the way to go. This should get you started:

test<-function(
    seed=1234,
    sample.size=500,
    resample.number=1000,
    alpha=0.05
    )
    {

        #initialize original sample
        original.sample<-runif(sample.size, min=0, max=1)   

        #initialize data.frame
        resample.results<-data.frame("Run.Number"=NULL,"mean"=NULL)
        for(counter in 1:resample.number){
            temp<-sample(original.sample, size=length(original.sample), replace = TRUE)
            temp.mean<-mean(temp)
            temp.table.row<-data.frame("Run.Number"=counter,"mean"=temp.mean)
            resample.results<-rbind(resample.results,temp.table.row)
        }
        resample.results<-resample.results[with(resample.results, order(mean)), ]

        #for the mean information
        lowerCI.row<-resample.number*alpha/2
        upplerCI.row<-resample.number*(1-(alpha/2))
        median.row<-resample.number/2

        #for the mean information
        median<-resample.results$mean[median.row]
        lowerCI<-resample.results$mean[lowerCI.row]
        upperCI<-resample.results$mean[upplerCI.row]

        #for the position information
        median.run<-resample.results$Run.Number[median.row]
        lowerCI.run<-resample.results$Run.Number[lowerCI.row]
        upperCI.run<-resample.results$Run.Number[upplerCI.row]

        mc.table<-data.frame("median"=NULL,"lowerCI"=NULL,"upperCI"=NULL)
        values<-data.frame(median,lowerCI,upperCI)
        #as.numeric because R doesn't like to mix data types
        runs<-as.numeric(data.frame(median.run,lowerCI.run,upperCI.run))
        mc.table<-rbind(mc.table,values)
        mc.table<-rbind(mc.table,runs)

        print(mc.table)
    }

After resampling your data, you find the mean. Then you order all of your resampled means. The middle of that list is the median. And, for example, with 10000 resamples, the 250th ordered resampled mean will be your lower 95% CI. Though I didn't do it here, the min value will just be at position 1, and the max value will be at position 10000. Be careful when you lower the resampling number: the way I calculated positions might become decimal values which will confuse R.

By the way, I put this in function form. If you like doing things line-by-line, just make sure to run all the lines between function() and the following main {}