**Avulsos by Penz - Articles tagged as R**, and kindly contributed to R-bloggers)

On the second article, we saw how to use a Monte Carlo simulation generate sample of disk space delta for future dates and calculate the distribution probability of zeroing free space in the future.

In this article, we will see how we can plot the evolution of predicted distribution for the occupied disk space. Instead of answering que question "how likely is that my disk space will zero before date X?," we will answer "how much disk space will I need by date X, and with what probability?"

# The input data

This file has the dataset we will use as example. It's the same we used in the second part. The graph below shows it:

We now import this data into R:

duinfo <- read.table('duinfospike.dat', colClasses=c("Date","numeric"), col.names=c("day","usd")) attach(duinfo) totalspace <- 450000 today <- tail(day, 1)

We then build our simulations for the next 4 months:

# Number of Monte Carlo samples numsimulations <- 10000 # Number of days to simulate numdays <- 240 # Simulate: simulate <- function(data, ndays) { delta <- diff(data) dssimtmp0 <- replicate(numsimulations, tail(data, 1)) dssimtmp <- dssimtmp0 f <- function(i) dssimtmp <<- dssimtmp + replicate(numsimulations, sample(delta, 1, replace=TRUE)) cbind(dssimtmp0, mapply(f, seq(1, ndays))) } dssim <- simulate(usd, numdays) # Future days: fday <- seq(today, today+numdays, by='day')

# Visualizing the possible scenarios

What king of data have we built in our simulations? Each simulation is built by sampling from the delta samples and adding to the current disk space for each day in the simulated period. We can say that each individual simulation is a possible scenario for the next 4 months. The graph bellow shows the first 5 simulations:

plot(fday, dssim[1,], ylim=c(min(dssim[1:5,]), max(dssim[1:5,])), ylab='usd', xlab='day', xaxt='n', type='l') axis.Date(1, day, at=seq(min(fday), max(fday), 'week'), format='%F') lines(fday, dssim[2,]) lines(fday, dssim[3,]) lines(fday, dssim[4,]) lines(fday, dssim[5,]) abline(h=totalspace, col='gray')

From this graph we can clearly see that the range of possible values for the used disk space grows with time. All simulations start with the same value - the used disk space for today - and grow apart as we sample from the delta pool.

We can also plot all simulations in a single graph:

plot(fday, dssim[1,], ylim=c(min(dssim), max(dssim)), ylab='usd', xlab='', xaxt='n', type='l') axis.Date(1, day, at=seq(min(fday), max(fday), 'week'), format='%F') f <- function(i) lines(fday, dssim[i,]) mapply(f, seq(2, numdays)) abline(h=totalspace, col='gray')

This plot gives us an idea of the overall spread of the data, but it fails to show the density. There are 10000 black lines there, with many of them overlapping one another.

# Visualizing the distribution for specific days

There is another way to look at our data: we have created, for each day, a sample of the possible used disk spaces. We can take any day of the simulation and look at the density:

dssimchosen <- list(density(dssim[,5]), density(dssim[,15]), density(dssim[,45]), density(dssim[,120])) colors <- rainbow(length(dssimchosen)) xs <- c(mapply(function(d) d$x, dssimchosen)) ys <- c(mapply(function(d) d$y, dssimchosen)) plot(dssimchosen[[1]], xlab='usd', ylab='dens', xlim=c(min(xs),max(xs)), ylim=c(min(ys),max(ys)), col=colors[1], main='') lines(dssimchosen[[2]], col=colors[2]) lines(dssimchosen[[3]], col=colors[3]) lines(dssimchosen[[4]], col=colors[4]) abline(v=totalspace, col='gray') legend('top', c('5 day', '15 days', '45 days', '120 days'), fill=colors)

By looking at this graph we can see the trend:

- The curves are getting flatter: we are getting more possible values for occupied disk space.
- The curves are moving to the right: we have more simulations with higher occupied disk space values.

# Visualizing the evolution of the distribution

So far, we have seen how we can visualize some simulations along the 4 months and how we can visualize the distribution for some specific days.

We can also plot the distribution of the values for each day in the simulated 4 months. We can't use the kernel density plot or the histogram, as they use both axes, but there are other options, most of them involving some abuse of the built-in plot functions.

## Boxplot

We can use the *boxplot* function to create a boxplot for each day in R in a
very straightforward way:

boxplot(dssim, outline=F, names=seq(today, as.Date(today+numdays), by='day'), ylab='usd', xlab='day', xaxt='n') abline(h=totalspace, col='gray')

The boxplots glued together form a shape that shows us the distribution of our simulations at any day:

- The thick line in the middle of the graph is the median
- The darker area goes from the first quartile to the third - which means that 50% of the samples are in that range
- The lighter area has the maximum and minimum points, if they are within 1.5 IQR of the upper/lower quartile. Points out of this range are considered outliers and are not plotted.

## Quantile lines

We can use the *quantile* function to calculate the values of each quantile per day, and plot the lines:

q <- 6 f <- function(i) quantile(dssim[,i], seq(0, 1, 1.0/q)) qvals <- mapply(f, seq(1, numdays+1)) colors <- colorsDouble(rainbow, q+1) plot(fday, qvals[1,], ylab='usd', xlab='day', xaxt='n', type='l', col=colors[1], ylim=c(min(qvals), max(qvals))) mapply(function(i) lines(fday, qvals[i,], col=colors[i]), seq(2, q+1)) axis.Date(1, day, at=seq(min(fday), max(fday), 'week'), format='%F') abline(h=totalspace, col='gray')

The advantage of this type of graph over the boxplot is that it is parameterized
by *q*. This variable tells us the number of parts that we should divide our
sample in. The lines above show us the division. If *q* is odd, the middle
line is exactly the median. If *q* is 4, the lines will draw a shape similar
to that of the boxplot, the only difference being the top and bottom line, that
will include outliers - the boxplot filters outliers by using the IQR as
explained above.

In the code above, we have used the *colorsDouble* function to generate a
sequence of colors that folds in the middle:

colorsDouble <- function(colorfunc, numcolors) { colors0 <- rev(colorfunc((1+numcolors)/2)) c(colors0, rev(if (numcolors %% 2 == 0) colors0 else head(colors0, -1))) }

## Quantile areas

We can also abuse the *barplots* function to create an area graph. We have to
eliminate the bar borders, zero the distance between them and plot a white bar
from the axis to the first quartile, if appropriate:

q <- 7 f <- function(i) { qa <- quantile(dssim[,i], seq(0, 1, 1.0/q)) c(qa[1], diff(qa)) } qvals <- mapply(f, seq(1, numdays+1)) colors <- c('white', colorsDouble(rainbow, q)) barplot(qvals, ylab='usd', xlab='day', col=colors, border=NA, space=0, names.arg=seq(min(fday), max(fday), 'day'), ylim=c(min(dssim), max(dssim))) abline(h=totalspace, col='gray')

In this case, using an odd *q* makes more sense, as we want to use the same
colors for the symmetric intervals. With an even *q*, there would either be a
larger middle interval with two quantiles or a broken symmetry. The code above
builds a larger middle interval when given an even *q*.

## Quantile heat map

If we increase *q* and use *heat.colors* in a quantile area plot, we get
something similar to a heat map:

q <- 25 f <- function(i) { qa <- quantile(dssim[,i], seq(0, 1, 1.0/q)) c(qa[1], mapply(function(j) qa[j] - qa[j-1], seq(2, q+1))) } qvals <- mapply(f, seq(1, numdays+1)) colors <- c('white', colorsDouble(heat.colors, q)) barplot(qvals, ylab='usd', xlab='day', col=colors, border=NA, space=0, names.arg=seq(min(fday), max(fday), 'day'), ylim=c(min(dssim), max(dssim))) abline(h=totalspace, col='gray')

# Visualizing past, present and future

We can also plot our data in the same graph as our simulations, by extending the
axis of the *barplot* and using the *points* function:

quantheatplot <- function(x, sim, ylim) { q <- 25 simstart <- length(x) - length(sim[1,]) f <- function(i) { if (i < simstart) replicate(q+1, 0) else { qa <- quantile(sim[,i-simstart], seq(0, 1, 1.0/q)) c(qa[1], diff(qa)) } } qvals <- mapply(f, seq(1, length(x))) colors <- c('white', colorsDouble(heat.colors, q)) barplot(qvals, ylab='usd', xlab='day', col=colors, border=NA, space=0, names.arg=x, ylim=ylim) abline(h=totalspace, col='gray') }

quantheatplot(c(day, seq(min(fday), max(fday), 'day')), dssim, ylim=c(min(c(usd, dssim)), max(dssim))) points(usd) abline(h=totalspace, col='gray')

# Training set, validation set

Cross-validation is a technique that we can use to validate the use of Monte Carlo on our data.

We first split our data in two sets: the training set and the validation set. We than use only the first in our simulations, and plot the second over. We can then see graphically if the data fits our simulation.

Let's use the first two months as the training set, and the other three months as the validation set:

# Number of days to use in the training set numdaysTrain <- 60 numdaysVal <- length(day) - numdaysTrain dssim2 <- simulate(usd[seq(1, numdaysTrain)], numdaysVal-1)

allvals <- c(usd, dssim2) quantheatplot(day, dssim2, c(min(allvals), max(allvals))) points(usd)

Looks like using only the first two months already gives us a fair simulation. What if we used only a single month, when no disk cleanup was performed?

# Number of days to use in the training set numdaysTrain <- 30 numdaysVal <- length(day) - numdaysTrain dssim3 <- simulate(usd[seq(1, numdaysTrain)], numdaysVal-1)

allvals <- c(usd, dssim3) quantheatplot(day, dssim3, c(min(allvals), max(allvals))) points(usd)

If we do regular disk cleanups, we must have at least one of them in our training set to get realistic results. Our training set is not representative without it.

This also tests our cross-validation code. A common mistake is using the whole data set as the training set and as the validation set. That is not cross-validation.

# Conclusions

We can use Monte Carlo simulations not only to generate a distribution probability of an event as we did in the previous article, but also to predict a possible range of future values. In this article, disk space occupation is not the most interesting example, as we are usually more interested in knowing when our used disk space will reach a certain value than in knowing the most probable values in time. But imagine that the data represents the number of miles traveled in a road trip or race. You can then not only see when you will arrive at your destination, but also the region where you will probably be at any day.

There are plenty of other uses for this kind of prediction. Collect the data, look at it and think if it would be useful to predict future ranges, and if it makes sense with the data you have. Predictions based on the evidence can be even used to support a decision or a point of view, just keep mind that you can only use the past if you honestly don't think anything different is going to happen.

**leave a comment**for the author, please follow the link and comment on his blog:

**Avulsos by Penz - Articles tagged as R**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...