Simulation experiments with "memisc"

For example, simulation results can be collected using the function Simulate and conveniently further processed to produce graphics such as in the following screenshot:

Output of *Simulate*

This shows how the law of large numbers can 'fail': Student's t-distribution with one degree of freedom (a Cauchy distribution) does not have a finite mean. Consequently averages of iid-random numbers with this distribution do not approach a finite value even if the sample size increases ad infinitum. (Actually, this example shows the relevance of the assumptions of finite variance and finite mean as premises of the law of large numbers.)

The code to produce the data and the plot is as follows:

library(memisc)

rtnorm <- function(n,df) if(df==Inf) rnorm(n) else rt(n,df)

# Run the simulation: t.normal.simres <- Simulate( mean(rtnorm(n,df)), # Conditions of the simulation experiment # that is, under what conditions the # replication series are run conditions = expand.grid( n=c(10,20,50,100,200,500,1000), df=c(Inf,5,3,2,1) ), nsim = 100, seed=4 # Fix the seed to make results reproducable )

# Turn the result into a data frame t.normal.simres <- as.data.frame(t.normal.simres)

# Some preparations for aggregation t.normal.simres <- within(t.normal.simres,{ df <- factor(df,levels=unique(df)) n <- factor(n) })

# Aggregate the results t.normal.sumry <- aggregate( c( mean=mean(result), sd=sd(result) )~df+n, data=t.normal.simres,sort=TRUE)

# Some preparations for the plot t.normal.sumry <- within(t.normal.sumry, type <- relabel(df, "Inf"="Normal", "5"="Student's t, df=5", "3"="Student's t, df=3", "2"="Student's t, df=2", "1"="Student's t, df=1" ) )

xyplot(cbind(mean,mean+2sd,mean-2sd)~n|type, data=t.normal.sumry, main="Mean", ylab="Mean of 100 replications +/- 2 standard deviations", xlab="Sample size", panel=panel.errbars, scales=list(y="free"), as.table=TRUE, type="o", pch=19)