R snippets

Missing tikzDevice

2013-04-30T13:42:00.000-07:00

I love using tikzDevice. When preparing LaTeX documents I switched to prepare all graphs in GNU R and then port them to TeX using tikzDevice. Recently I have moved to GNU R 3.0.0 and was shocked to find that this package is no longer available on CRAN.
Fortunately I was able to find package on the web. Unfortunately - installing old package on GNU R version 3.0.0 simply generates an error when trying to load it.

So for now the only solution for me was to:

Install GNU R 2.15.3 in parallel to new one (fortunately this is not a problem)
downloading tikzDevice 0.6.3 which can be found at http://download.r-forge.r-project.org/bin/windows/contrib/2.15/
installing tikzDevice manually from local disk
installing filehash package that it requires

Fortunately this worked.

I really hope that tikzDevice will be back to CRAN in the future.

Estimating continuous piecewise linear regression

2013-04-02T08:02:00.002-07:00

When talking about smoothing splines a simple point to start with is a continuous piecewise linear regression with fixed knots. I did not find any simple example showing how to estimate the it in GNU R so I have created a little snippet that does the job.

Assume you are given continuous predictor x and continuous predicted variable y. We want to estimate continuous piecewise linear regression with fixed knots stored in variable knots using standard lm procedure.

The key to a solution is proper definition of regression formula. In order to introduce possibility of change of slope in knot k we have to add a so called hinge term to the model max(0, x-k).

In the code given below function piece.formula automatically generates a proper right hand side of the regression formula given variable name and list of required knots. It is next tested on a simple function.

N <- 40 # number of sampled points
K <- 5  # number of knots

piece.formula <- function(var.name, knots) {
  formula.sign <- rep(" - ", length(knots))
  formula.sign[knots < 0] <- " + "
  paste(var.name, "+",
        paste("I(pmax(", var.name, formula.sign, abs(knots), ", 0))",
              collapse = " + ", sep=""))
}

f <- function(x) {
    2 * sin(6 * x)
}

set.seed(1)
x <- seq(-1, 1, len = N)
y <- f(x) + rnorm(length(x))

knots <- seq(min(x), max(x), len = K + 2)[-c(1, K + 2)]
model <- lm(formula(paste("y ~", piece.formula("x", knots))))

par(mar = c(4, 4, 1, 1))
plot(x, y)
lines(x, f(x))
new.x <- seq(min(x), max(x) ,len = 10000)
points(new.x, predict(model, newdata = data.frame(x = new.x)),
      col = "red", pch = ".")
points(knots, predict(model, newdata = data.frame(x = knots)),
       col = "red", pch = 18)

Below we can see the graph of estimation result. Red line is the desired continuous piecewise linear regression with fixed knots given by red diamonds. Notice that the plot uses points procedure to plot the red line to highlight that the generated predictions have the required properties.

An additional value of the presented solution that we do not do any preprocessing of predictor variable if we want to make a prediction - all the calculations are made within the formula.

Of course this simple example can be easily extended to obtain a simple smoother. For example we can set K to be large and use some regularized regression like ridge or lasso.

GNU R loop speed comparison

2013-03-16T05:26:00.001-07:00

Recently I had several discussions about using for loops in GNU R and how they compare to *apply family in terms of speed. I have not seen a direct benchmark comparing them so I decided to execute one (warning: some of the code presented today takes long time to execute).

First I have started by comparing the speed of assignment operator for lists vs. numeric vectors in standard for loops. Here is the code:

speed.test <- function(n) {
  gc()
  x1 <- numeric(n)
  x2 <- vector(n, mode = "list")
  c(system.time(for (i in 1:n) { x1[i] <- i })[3],
    system.time(for (i in 1:n) { x2[[i]] <- i })[3])
}

n <- seq(10 ^ 4, 10 ^ 6, len = 5)
result <- t(sapply(n, speed.test))
par(mar=c(4.5, 4.5, 1, 1))
matplot(n / 1000, result, type = "l" , col = 1:2, lty = 1,
        xlab = "n ('000)", ylab = "time")
legend("topleft", legend = c("numeric", "list"),
       col = 1:2, lty = 1)

The picture showing the result of the comparison is the following:

As we can see - operation numeric vectors are significantly faster than list, especially for large vector sizes.

But how does this relate to *apply family of functions? The issue is that the workhorse function there is lapply and it works on lists. Other functions from this family call lapply internally.

So I have run the second test comparing: (a) lapply, (b) for loop working on lists and (c) for loop working on numeric vectors. Here is the code:

aworker <- function(n) {
  r <- lapply(1:n, identity)
  return(NULL)
}

lworker <- function(n) {
  r <- vector(n, mode = "list")
  for (i in 1:n) {
    r[[i]] <- identity(i)
  }
  return(NULL)
}

nworker <- function(n) {
  r <- numeric(n)
  for (i in 1:n) {
    r[i] <- identity(i)
  }
  return(NULL)
}

run <- function(n, worker) {
  gc()
  unname(system.time(worker(n))[3])
}

compare <- function(n) {
  c(lapply = run(n, aworker),
    list = run(n, lworker),
    numeric = run(n, nworker))
}

n <- rep(c(10 ^ 6, 10 ^ 7), 10)
result <- t(sapply(n, compare))
par(mfrow = c(1,2), mar = c(3,3,3,1))
for (i in n[1:2]) {
  boxplot(result[n == i,],
          main = format(i, scientific = F, big.mark = ","))
}

On the picture below we can see the result. For 1,000,000 elements in a vector lapply is the fastest. The reason it that it executes looping in compiled C code. However for 10,000,000 elements for loop using numeric vector is faster as it avoids conversion to list.

Of course probably on other machines than my notebook the difference in speed would manifest itself for other number of elements in a vector.
However one can draw a general conclusion: if you have large AND numeric vectors and need to do a lot of number crunching for loop will be faster than lapply.

Hexadecimal literals in GNU R

2013-03-11T05:04:00.002-07:00

Recently I have used hexadecimal numbers in GNU R. The way they are parsed surprised me and is inconsistent with Java. As R Language Definition pdf only briefly mentions hexadecimal numbers here is what I have found.

First I have checked the following code:

0x11

# 17

0x1.1

# 17

I seemed that GNU R drops the decimal point. So I have checked notation with exponent being power of 2 (as in Java):

0x1.1p0

# 1.0625

0x1.1p4

# 17

Now decimal point is correctly taken into account so I have checked further.

In the notation without exponent part GNU R just discards decimal points as in this example:

0x111

# 273

0x1.1.1

# 273

0x....1....1....1

# 273

However, in notation with exponent last decimal point is taken into consideration and all earlier decimal points are ignored:

0x11.1p0

# 17.0625

0x1.1.1p0

# 17.0625

0x....1....1....1p0

# 17.0625

In summary - handling of hexadecimal literals in GNU R is nonstandard and at least for me - unexpected.

Times per second benchmark

2013-03-05T15:41:00.004-08:00

In GNU R the simplest way to measure execution time of a piece code is to use system.time. However, sometimes I want to find out how many times some function can be executed in one second. This is especially useful when we want to compare functions that have significantly different execution speed.

Fortunately times per second benchmark for execution time can be simply evaluated using the following snippet:

tps <- function(f, time) {

gc()

i <- 0

start <- proc.time()[3]

repeat {

i <- i + 1

f(i)

stop <- proc.time()[3]

if (stop - start > time) {

return (i / (stop - start))

}

}

This function takes two parameters: a function to be benchmarked (f) and how much time is to be used for evaluation (time). It returns an estimate how many times per second function f can be executed.

As a simple application of tps function consider calculating relative speed of standard, lattice and ggplot2 graphics. The following function compares them by plotting histograms:

library(ggplot2)

library(lattice)

test <- function(n, time) {

x <- runif(n)

b <- c(tps(function(i) {

hist(x, 10, main = i)

}, time),

tps(function(i) {

print(histogram(x, nbin = 10, main = format(i)))

}, time),

tps(function(i) {

print(qplot(x, binwidth=0.1, main = i))

}, time))

names(b) <- c("hist", "histogram", "qplot")

return(b)

}

The function takes two arguments. First is number of points to sample for the histogram and second is time passed to tps function. On my computer the test gave the following result for 10000 size of the sample and 5 seconds for each function each :

> test(10000, 5)

hist histogram qplot

192.614770 14.285714 5.544933

We can see that standard hist is over 10 times faster than from histogram from lattice and almost 40 times faster than qplot from ggplot2.

Solving 9-puzzle with GNU R

2013-01-04T14:23:00.001-08:00

During holiday break I have decided to solve 9-puzzle, which is 3x3 variant of a well known 15-puzzle. The solution has proven to be a nice application of igraph package. Warning: this time the code takes a bit more tame than usual in my posts to execute.
9-puzzle is a problem where numbers from 1 to 8 are placed on a board that has size 3x3 tiles. One spot is left empty. We can slide numbers to empty tile. Starting from an arbitrary setup we want to put the numbers in order: 1, 2, 3 in the first row, 4, 5, 6 in the second row and 7,8 and empty in third row.

To solve it I represent each positioning of the numbers as vertex in the graph. Two vertexes are connected by an edge (undirected) if one of them can be reached from the other using only one move.

In order to encode the number setups I use 9 for empty space. Then each position can be encoded as permutation of numbers from 1 to 9. The encoding is row wise, so target setup has 123456789 encoding. Here is the code that writes to a file all edges in the graph:

library(gtools)

moves <- list(c(2, 4), c(3, 5), 6, c(5, 7), c(6, 8), 9, 8, 9)

cat("Initializing permutations\n\n")

flush.console()

positions <- permutations(9, 9)

pos9 <- apply(positions, 1, function(x) { which(x == 9) } )

max.val <- length(which(pos9 != 9))

edge.write <- function(i) {

percentage <- i / max.val

setTxtProgressBar(progress.bar, percentage)

pos9i <- pos9[i]

if (pos9i == 9) {

stop("Should not happen")

}

from <- positions[i,]

to <- matrix(from, nrow = length(moves[[pos9i]]),

ncol = 9, byrow = T)

for (j in 1:nrow(to)) {

to[j, pos9i] <- to[j, moves[[pos9i]][j]]

to[j, moves[[pos9i]][j]] <- 9

cat(paste(from, collapse=""), ",",

paste(to[j,],collapse=""), "\n", file = out.file)

}

file.name <- tempfile()

cat("Writing edge file to: ", file.name, "\n",

rep(".", 50), "\n", sep="")

progress.bar <- txtProgressBar(0, 1, 0, width = 50)

out.file <- file(file.name, "w")

cat("from, to\n", file = out.file)

for (i in which(pos9 != 9)) {

edge.write(i)

}

close(progress.bar)

close(out.file)

Notice that the graph is undirected so we save only moves of empty tile "right" or "down". This is encoded in moves variable. For example if position 2 is empty it can be moved to position 1, 3 or 5, but only 3 and 5 are encoded as move to 1 is "left". In order to optimize the algorithm we remove all setups that have empty tile at position 9 as there are no moves "right" or "down" from such a setup.

Once the file is written it can be processed using igraph package as follows:

library(igraph)

cat("\nReading edge file to iGraph\n")

flush.console()

edgelist <- read.csv(file.name)

unlink(file.name)

puzzle9.graph <- graph.data.frame(edgelist, directed = TRUE)

cat("Number of clusers in the graph:",

clusters(puzzle9.graph)$no, "\n\n")

v.names <-get.vertex.attribute(puzzle9.graph,"name")

v.start <- which(v.names=="123456789")

sp <- shortest.paths(puzzle9.graph,v.start)

v.tough <- which(sp==max(sp[sp<Inf]))

formatted.tough <- character(2)

for (i in 1:length(v.tough)) {

split.tough <- strsplit(v.names[v.tough[i]],split="")[[1]]

split.tough[which(split.tough=="9")] <- "."

split.tough <- matrix(split.tough, nrow = 3)

formatted.tough[i] <- paste(apply(split.tough, 2, paste,

collapse=""),

collapse="\n")

}

par(mar=c(4.5, 4.5, 0.5, 0.5))

plot(prop.table(table(sp, exclude=Inf)),

xlab = "Steps needed", ylab = "Probability")

text(2, 0.12,"Tough cases:", family = "mono", pos = 4)

text(2, 0.1, formatted.tough[1], family = "mono", pos = 4)

text(9, 0.1, formatted.tough[2], family = "mono", pos = 4)

cat("Done!\n")

What we learn is that the graph has 2 clusters and only half of the possible setups can lead us to desired orderings (in fact we know that it is directly related to parity of the permutation). The code produces the graph given below. We see that most of the times 22-24 moves are needed to solve the puzzle if it is possible and in worst case 31 moves are needed. We find that there exactly 2 such tough cases and they are printed out on a graph:

Wrapper functions in GNU R

2012-12-26T06:40:00.002-08:00

Recently I have been working with GNU R optimization routines a lot. Function optim has a nice trace option that allows to monitor optimization progress. Another standard function optimize has no such feature but it is possible to add outside tracking to it.

Let us start with an example from documentation:

f <- function (x, a) { (x - a) ^ 2 }

print(optimize(f, c(0, 1), tol = 0.0001, a = 1/3))

which produces:

$minimum

[1] 0.3333333

$objective

[1] 0

We only get the final result of the optimization. The simple way to add tracing is by inserting print statement inside optimized function like this:

f <- function (x, a) {

print(x)

(x - a) ^ 2

}

print(optimize(f, c(0, 1), tol = 0.0001, a = 1/3))

The call produces the following output:

[1] 0.381966

[1] 0.618034

[1] 0.236068

[1] 0.3333333

[1] 0.3333

[1] 0.3333667

[1] 0.3333333

$minimum

[1] 0.3333333

$objective

[1] 0

However, GNU R allows for a much fancier method that at the same time shows some nice features of the language itself. We use original function but add a wrapper around it:

wrap.f <- function(f, minimum = TRUE) {

calls <- NULL

return(list(f = function(x, ...) {

calls <<- c(calls, x)

ifelse(minimum,f(x, ...), -f(x, ...))

getter = function() { return(calls) },

reset = function() { calls <<- NULL }))

}

wf <- wrap.f(function (x, a) { (x - a) ^ 2 })

print(optimize(wf$f, c(0, 1), tol = 0.0001, a = 1/3))

print(wf$getter())

wf$reset()

print(wf$getter())

The wrap.f function allows for tracking of calls of its first argument and to change optimization direction (this is possible through optimize function, but I add this as an example what can be achieved using this method). It holds internal variable calls which keeps track of calls to f function. It can be accessed using getter function and reset using reset function. The last three lines of code present the effect of resetting calls variable. Here is the output from the code:

$objective

[1] 0

[1] 0.3819660 0.6180340 0.2360680 0.3333333 0.3333000 0.3333667 0.3333333

NULL

The wrapping method presented here can be used as a replacement for S4 classes in simple situations. We use the environment created by wrap.f function to store attributes and return methods that can access and modify them.

Predictive models in R: a new book in Polish

2012-12-14T14:29:00.000-08:00

Together with Mateusz Zawisza I have just published a new book in Polish on building predictive models in GNU R. It can be bought at Oficyna Wydawnicza SGH. The book presents complete examples of basic data mining processes.
Although the book is in Polish, sources of all procedures used in it, which are available on my website, can be used without the book. Here is a simplified code from exercise 4.5 presenting neural network bagging:

library(nnet)

set.seed(1)

SAMPLE_SIZE <- 256

X <- seq(-2, 2, length.out = SAMPLE_SIZE)

TRUE_Y <- X ^ 2 / 2 + sin(4 * X)

y <- TRUE_Y + 2 * rnorm(SAMPLE_SIZE)

GetBootstrapPrediction <- function() {

bootstrap.indices <- sample(SAMPLE_SIZE, replace = T)

bootstrap.sample.y <- y[bootstrap.indices]

bootstrap.sample.x <- X[bootstrap.indices]

bootstrap.model <- nnet(bootstrap.sample.y ~bootstrap.sample.x,

lin = T, size = 4, trace = FALSE, maxit = 10 ^ 6)

return(predict(bootstrap.model, data.frame(bootstrap.sample.x = X)))

}

progress.bar <- winProgressBar("Progress in %", "0% done", 0, 1, 0)

BOOTSTRAP_REPLICATIONS <- 1024

bootstrap.predictions <- rep(0, SAMPLE_SIZE)

for (i in 1:BOOTSTRAP_REPLICATIONS) {

bootstrap.predictions <-bootstrap.predictions +

GetBootstrapPrediction()

percentage <- i / BOOTSTRAP_REPLICATIONS

setWinProgressBar(progress.bar, percentage, "Progress in %",

sprintf("%d%% done", round(100 * percentage)))

}

close(progress.bar)

plot(X, y,xlim = c(-2, 2), ylim = c(-5, 6))

lines(X, TRUE_Y, lwd = 4)

lines(X, bootstrap.predictions /BOOTSTRAP_REPLICATIONS,
lwd = 3, col = 3)

It produces the following graph. Circles represent training data, black line is true relationship and green line is prediction from bagging procedure:

Bridge hand distribution: simulation vs exact calculation

2012-12-08T02:42:00.001-08:00

Recently I played bridge with my friends. Being frustrated with several consecutive poor hand distributions we asked ourselves a question what is the probability of having a hand good enough for a small slam. A well known rule of thumb is that you need 33+ HCP for 6NT. But we could not find information about the probability of such an event. So we decided to calculate it.
And the question was whether to calculate it exactly or simulate the result. With GNU R you can do both. However - not surprisingly - simulation is much easier and faster to perform and exact calculation is coding error prone. Have a look at the codes.

First start with simulated result. The code is clean and simple:

set.seed(1)

deck <- rep(c(1:4, rep(0, 9)), 4)

points <- replicate(1000000, sum(sample(deck, 26)))

print(2* mean(points > 32))

and we get the following result:

[1] 0.00687

On the other hand exact result requires careful counting of all possible card combinations:

library(caTools)

result <- data.frame(HCP=0:40, exact = 0)

result[1, 2] <- choose(16, 0) * choose(36, 26)

for (i in 1:16) {

HCP.sums <- table(rowSums(combs(rep(1:4,4), i))) *

choose(36, 26 - i)

levels <- as.integer(names(HCP.sums)) + 1

for (i in 1:length(HCP.sums))

result[levels[i], 2] <- result[levels[i], 2] +

HCP.sums[i]

}

result[, 2] <- result[, 2] / sum(result[, 2])

print(2* sum(result$exact[result$HCP > 32]))

It took: (a) some thinking before coding, (b) 10x more time to code and (c) one mistake in the process to get the desired results. The output of the procedure is the following:

[1] 0.0069593

So the conclusion is that some pair gets a hand good enough for 6NT on the average once in 150 games.

I wanted check sure whether the simulation and exact results are similar so I decided to plot the resulting HCP distributions using the code:

par(mar=c(4, 4, 1, 1))

result$sim <- sapply(0:40, function(x) { mean(points == x) })

plot(result$HCP, result$exact, pch = 4,

xlab = "HCP", ylab = "Probability")

points(result$HCP, result$sim, pch = 3, col ="red")

Tthe result given below shows that simulation approximates exact result pretty well.

ESSA2013 Conference

2012-11-24T03:10:00.000-08:00

It has been just announced that during ESSA2013 conference I am planning to organize a special track on "Statistical analysis of simulation models". I hope to get some presentations using GNU R to promote it in social simulation community.

It is obvious that GNU R excels in analysis of simulation data. However, very often it can be neatly used to implement simulations themselves.

For instance I have recently implemented a simulation model proposed in Section 4 of Volatility Clustering in Financial Markets: Empirical Facts and Agent–Based Models paper by Rama Cont. The model is formulated as follows (I give only its brief description; please refer to the paper for more details).

Consider market with n trading agents and one asset. We simulate market for times periods. In each period each agent can buy asset, sell it or do nothing.
Asset return r[i] in period i equals to number of buy orders minus number of sell orders divided by number of agents n and multiplied by normalizing constant max.r. Thus it will always lie in the interval [-max.r,max.r].
Agents make buy and sell decisions based on random public information about an asset. The stream of signals are IID normal random variables with mean 0 and standard deviation signal.sd. Each investor holds an internal non negative decision making threshold. If signal is higher than threshold level buy decision is made. If it is lower than minus threshold level asset is sold. If signal is not strong enough investor does nothing.
After return r[i] is determined each investor with probability p.update performs threshold update to abs(r[i]).
As you can see the description is quite lengthily. However, the implementation of the model in GNU R is a genuine snippet as can be seen below:

cont <- function(times, n, signal.sd, max.r ,p.update) {

threshold <- vector("numeric", n)

signal <- rnorm(times, 0, signal.sd)

r <- vector("numeric", times)

for (i in 1:times) {

r[i] <- max.r * (sum(signal[i] > threshold)

- sum(signal[i] < (-threshold))) / n

threshold[runif(n) < p.update] <- abs(r[i])

}

return(r)

}

And an additional benefit is that one can analyze the simulation results in GNU R also. Here is a very simple example showing the relationship between signal.sd and standard deviation of simulated returns (the initial burn in period in the simulation is discarded):

cont.sd <- function(signal.sd) {

sd(cont(10000, 1000, signal.sd, 0.1, 0.05)[1000:10000])

}

sd.in <- runif(100, 0.01, 0.1)

sd.out <- sapply(sd.in, cont.sd)

plot(sd.in,sd.out)

and here is the resulting plot:

Possible error in Bayesian bootstrap

2012-11-08T05:51:00.003-08:00

After my last post on Bayesian bootstrap I got a question why the sample from Dirichlet distribution is taken as weights for calculating mean in the procedure and not as weights used for sampling from the original data set. Actually this mistake is subtle and occurs even in textbooks, see example Chernick (2008), page 122. In this post I want to clarify the issue.

In the example I give correct bootstrap and Bayesian bootstrap procedures and wrong ones. The wrong Bayesian bootstrap follows description from Chernick (2008), page 122 (that is equivalent to the comment to my last post).

Here is the code that I used:

library(gtools)

ok.mean.bb <- function(x, n) {

apply(rdirichlet(n, rep(1,length(x))), 1, weighted.mean, x = x)

}

ok.mean.fb <- function(x, n) {

replicate(n, mean(sample(x, length(x), TRUE)))

}

wrong.mean.bb <- function(x, n) {

replicate(n, mean(sample(x, length(x), TRUE,

diff(c(0, sort(runif(length(x) - 1)), 1)))))

}

wrong.mean.fb <- function(x, n) {

replicate(n, mean(sample(sample(x, length(x), TRUE),

length(x), TRUE)))

}

set.seed(1)

reps <- 10000

x <- cars$dist

par(mar=c(5,4,1,2))

plot(density(ok.mean.fb(x, reps)), main = "", xlab = "Bootstrap mean")

lines(density(ok.mean.bb(x, reps)), col = "red")

lines(density(wrong.mean.fb(x, reps)), col = "blue")

lines(density(wrong.mean.bb(x, reps)), col = "green")

The figure it produces is:

Black curve is standard bootstrap density, red is Bayesian bootstrap and blue and green are generated by wrong bootstrapping procedures (respectively frequentist and Bayesian).

We can see that wrong Bayesian bootstrap has an equivalent in standard bootstrap approach that is generated by repeating the sampling twice (sampling from a sample) and it clearly increases dispersion of the results.

Simple Bayesian bootstrap

2012-11-02T16:01:00.001-07:00

Bootstrapping is a very popular statistical technique. However, its Bayesian analogue proposed by Rubin (1981) is not very common. I was looking for an example of its implementation in GNU R and could not find one so I decided to write a snippet presenting it.

In standard bootstrapping observations are sampled with replacement. This implies that observation weights follow multinomial distribution. In Bayesian bootstrap multinomial distribution is replaced by Dirichlet distribution.

This observation leads to very simple implementation of Bayesian bootstrap using gtools package. Here is the code presenting it with a simple application giving frequentist and Bayesian 95% confidence interval for mean in fbq and bbq variables:

library(gtools)

# Bayesian bootstrap
mean.bb <- function(x, n) {

apply(rdirichlet(n, rep(1, length(x))), 1, weighted.mean, x = x)

}

# standard bootstrap
mean.fb <- function(x, n) {

replicate(n, mean(sample(x, length(x), TRUE)))

}

set.seed(1)

reps <- 100000

x <- cars$dist

system.time(fbq <- quantile((mean.fb(x, reps)), c(0.025, 0.075)))

system.time(bbq <- quantile((mean.bb(x, reps)), c(0.025, 0.075)))

As it can be seen implementation of Bayesian bootstrap is fairly simple.

On my computer Bayesian bootstrap is approximately 80% slower than standard bootstrap, but its performance probably could be improved.

Plotting Watts-Strogatz model

2012-09-21T07:29:00.001-07:00

Recently I wanted to reproduce Figure 2 from Watts and Strogatz (1998). The task using igraph is simple but an interesting task was annotation of the resulting plot.

Watts-Strogatz model generates graphs that have so called small-world network property. Such networks should have low average path length and high clustering coefficient. The algorithm has three parameters: number of nodes in the graph, initial number of neighbors of each node distributed on a ring and rewiring probability.

Interestingly in Watts-Strogatz model having small but positive values of rewiring probability generates graphs having desired properties - and this is exactly depicted on Figure 2 in their article.

I decided to replicate it. To enhance it I wanted to plot median and 5 and 95 percentile of distribution of average path length and clustering coefficient as a function of rewiring probability.

Here you have the code that generates the graph (warning: it takes about 1 minute to run):

library(igraph)

set.seed(1)

avg.stat <- function(nei, p) {

result <- replicate(1000, {

wsg <- watts.strogatz.game(1, 100, nei, p)

c(average.path.length(wsg),

transitivity(wsg))

})

apply(result, 1, quantile, probs = c(0.5, 0.05, 0.95))

}

nei <- 6

p <- 2 ^ -seq(0, 10, len = 21)

result <- sapply(p, avg.stat, nei = nei)

result <- t(result / rep(avg.stat(nei, 0)[1,], each = 3))

par(mar=c(3.2, 2, 0.2, 0.2), mgp=c(2, 1, 0))

matplot(p, result, type = "l", log = "x", xaxt = "n", ylab = "",

lty = rep(c(1,2,2),2), col=rep(c(1,2), each=3))

axis(1, at = 2 ^ -(0:10),

labels = c(1, parse(text = paste(2, 1:10, sep = "^-",

collapse = ";"))))

legend("bottomleft", c("average path length", "clustering coefficient"),

lty = 1, col = c(1, 2))

The result of the procedure is the following picture:

It looks very similar to what is shown in the article (apart from adding lines depicting 5 and 95 percentile of distributions of both graph characteristics).

However, the interesting part was to properly annotate X-axis on the plot. Of course you can use expression function to get it but then the problem is that you have to do it ten times. Interestingly parsing a string containing those ten expressions separated by semicolons works just as needed.

Simulation metamodeling with constraints

2012-09-14T07:45:00.002-07:00

Last week I have posted about using simulation metamodeling to verify results of analytical solution of the model. After posting it I realized that the solution presented there can be improved by using knowledge of simulation model structure.

The story is about calculation of the variance of the random variable given by formula:

q(p-2+q/N)

where p and N are model parameters and q is a random variable that has binomial distribution with N Bernoulli trials and probability of success 2-p. As I have written the properly specified formula for the metamodel has the following form:

v ~ (p + I(p ^ 2) + I(p ^ 3) + I(p ^ 4)) * (N + I(1 /N))

where v is an estimate of the variance. Such a model was discussed in my last post (it is described in more detail there so I omit the repetition in this post).

However one can easily notice that for p=1 and p=2 the variance of the random variable is equal to 0. This is because in such case q has no variance (it is equal to 1 and 0 respectively).

In estimation of the metamodel one can take this into account in two ways:

the weight of cases where p equals 1 or 2 might be increased (so we use Weighted Least Squares)
a restriction on values of prediction in points where p equals 1 or 2 can be made to ensure that it is exactly 0 (so we use Constrained Least Squares)

The code given below estimates: original model and two versions of corrected models:

library(plyr)

library(limSolve)

sim <- function(param) {

q <- rbinom(1000, param$N, 2 - param$p)

c(v = var((param$p - 2 + q / param$N) * q))

}

set.seed(1)

data.set <- expand.grid(p = seq(1, 2, len = 101), N = 1:10)

data.set <- ddply(data.set, .(p,N), sim,.progress = "text")

lm.coef <- coef(lm(v ~ (p + I(p ^ 2) + I(p ^ 3) + I(p ^ 4))

* (N + I(1 / N)), data = data.set))

weights <- ifelse(data.set$v==0, 100000, 1)

wlm.coef <- coef(lm(v ~ (p + I(p ^ 2) + I(p ^ 3) + I(p ^ 4))

* (N + I(1 / N)), data = data.set,

weights = weights))

A <- with(data.set, cbind(1, p, p^2, p^3, p^4, N, 1 / N,

p * N, p / N, p ^ 2 * N, p ^ 2 / N,

p ^ 3 * N, p ^ 3 / N, p ^ 4 * N, p ^ 4 / N))

colnames(A) <- names(lm.coef)

B <- data.set$v

E <- A[B == 0,]

F <- B[B == 0]

lsei.coef <- lsei(A, B, E, F)$X

Now we can compare the coefficients from these models to their analytical values given in the last post:

true.coef <- c(32, -88, 88, -38, 6, -8, -26, 20,

75, -18, -79, 7, 36, -1, -6)

coef.dev <- cbind(lm.coef, wlm.coef, lsei.coef) - true.coef

print(coef.dev, digits = 2)

dotchart(coef.dev[,1], pch = 19)

points(coef.dev[,2], 1:nrow(coef.dev), col = "red", pch = 19)

abline(v = 0, col = "gray", lty = 2)

The code first tabulates the deviation of estimates from true values:

lm.coef wlm.coef lsei.coef

(Intercept) -0.6056 0.060 0.060

p 1.5061 -0.235 -0.235

I(p^2) -1.3576 0.323 0.323

I(p^3) 0.5243 -0.185 -0.185

I(p^4) -0.0731 0.038 0.038

N 0.0088 -0.060 -0.060

I(1/N) 1.1510 0.354 0.354

p:N -0.0065 0.173 0.173

p:I(1/N) -3.0646 -0.965 -0.965

I(p^2):N -0.0102 -0.182 -0.182

I(p^2):I(1/N) 2.9942 0.953 0.953

I(p^3):N 0.0115 0.084 0.084

I(p^3):I(1/N) -1.2730 -0.404 -0.404

I(p^4):N -0.0029 -0.014 -0.014

I(p^4):I(1/N) 0.1991 0.062 0.062

and next plots the visualization of the differences (on the plot black points represent OLS estimates and red - WLS):

It can be seen that:

the precision of estimation is improved by using additional information on problem structure;
the estimates from WLS and CLS are identical to the third decimal place.

Simulation metamodeling with GNU R

2012-09-07T13:42:00.001-07:00

I am one of the organizers of ESSA2013 conference that will take place in September 2013 in Warsaw, Poland. The conference scope is social simulation and in particular methods of statistical analysis of simulation output (metamodeling). As we have just issued Call for Papers for the conference so I decided to post a simple example of a metamodel.

Recently I had to calculate variance of a random variable given as:

q(p-2+q/N),

where p and N are model parameters and q is a random variable that has binomial distribution with N Bernoulli trials and probability of success 2-p. After some lengthily calculations I came up with the following formula for the variance:

However the result is so long that I was not sure that it was correct, so I decided to verify it using metamodeling.

First I have generated the data for the simulation with the following code and visualized it. The data covers a grid where p is from interval [1;2] and N spans from 1 to 10. In each grid point 100000 simulations are performed and variance of q(p-2+q/N) is calculated:

library(plyr)

library(lattice)

sim <- function(param) {

q <- rbinom(100000, param$N, 2 - param$p)

c(v = var((param$p - 2 + q / param$N) * q))

}

set.seed(1)

data.set <- expand.grid(p = seq(1, 2, len = 201), N = 1:10)

data.set <- ddply(data.set, .(p,N), sim, .progress = "text")

levelplot(v ~ p + N, data = data.set,

col.regions = terrain.colors, xlab = "p", ylab = "N")

Here is the resulting plot:

Next I have estimated the linear regression model with parameters reflecting analytical specification given above:

summary(lm(v~(p + I(p ^ 2) + I(p ^ 3) + I(p ^ 4))

* (N + I(1 / N)), data = data.set))

The output looks as follows:

Call:

lm(formula = v ~ (p + I(p^2) + I(p^3) + I(p^4)) * (N + I(1/N)),

data = data.set)

Residuals:

Min 1Q Median 3Q Max

-0.0116801 -0.0007195 0.0000063 0.0007608 0.0098659

Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) 32.089136 0.235380 136.33 <2e-16 ***

p -88.232888 0.655760 -134.55 <2e-16 ***

I(p^2) 88.225183 0.674590 130.78 <2e-16 ***

I(p^3) -38.095417 0.303854 -125.37 <2e-16 ***

I(p^4) 6.014949 0.050597 118.88 <2e-16 ***

N -8.012673 0.027781 -288.42 <2e-16 ***

I(1/N) -26.082527 0.303365 -85.98 <2e-16 ***

p:N 20.034053 0.077398 258.84 <2e-16 ***

p:I(1/N) 75.213337 0.845166 88.99 <2e-16 ***

I(p^2):N -18.033947 0.079621 -226.50 <2e-16 ***

I(p^2):I(1/N) -79.203708 0.869434 -91.10 <2e-16 ***

I(p^3):N 7.014852 0.035863 195.60 <2e-16 ***

I(p^3):I(1/N) 36.085047 0.391618 92.14 <2e-16 ***

I(p^4):N -1.002405 0.005972 -167.85 <2e-16 ***

I(p^4):I(1/N) -6.013089 0.065212 -92.21 <2e-16 ***

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.002197 on 1995 degrees of freedom

Multiple R-squared: 0.9999, Adjusted R-squared: 0.9999

F-statistic: 1.976e+06 on 14 and 1995 DF, p-value: < 2.2e-16

As we can see the model fit is almost perfect and the estimates are very close to analytical results. This reassured me that my calculation was correct.

In this way you could use simulation to verify analytical results. Of course simulation metamodeling has much broader applications - especially when we do not have analytical results.

I hope to organize a special session during ESSA2013 dedicated to simulation metamodeling using GNU R. So if you are interested in such topics and willing to promote GNU R in simulation society you are welcome to come to Warsaw.

Exporting ctree object to Asymptote

2012-08-25T02:27:00.001-07:00

When producing regression or classification trees (standard rpart or ctree from party package) in GNU R I am often unsatisfied with the default plots they produce. One of many possible solutions is to export a tree plot to Asymptote.
The code I have prepared generates an Asymptote file based on generated ctree object. Here is the procedure that does the conversion.

treeAsy <- function(tree, # ctree to be plotted

off.f, # tree plot fixed shift

off.v, # tree plot variable shift

file.name, # output file name

preamble) { # preamble for asy

minv <- +Inf

maxv <- -Inf

response <- names(tree@responses@variables)

plot.node <- function(root,

nest = 0, # level in a tree

pOffset = 0, # plotting offset

condition = "root", # split condition text

id = "root") { # block name in asy

if (length(root$prediction) > 1) {

stop("Only single prediction value supported")

}

if (root$prediction < minv) minv <<- root$prediction

if (root$prediction > maxv) maxv <<- root$prediction

child.l <- ""

child.r <- ""

if (!root$terminal) {

if (class(root$psplit) == "orderedSplit") {

varN <- root$psplit$variableName

point <- root$psplit$splitpoint

left <- paste(varN, "$\\leq$", point, sep="")

right <- paste(varN, "$>$", point, sep="")

} else {

stop("Only orderedSplit supported")

}

add <- "add(new void(picture pic, transform t) {

blockconnector operator --=blockconnector(pic,t);\n "

child.l <- paste(plot.node(root$left, nest + 1,

pOffset - off.f - 1 / off.v ^ nest,

left, paste(id,"l",sep="")),

add,id,"--Down--Left--Down--",id,"l;\n});\n\n", sep="")

child.r <- paste(plot.node(root$right, nest + 1,

pOffset + off.f + 1 / off.v ^ nest,

right, paste(id,"r",sep="")),

add,id,"--Down--Right--Down--",id,"r;\n});\n\n", sep="")

}

paste("block ", id, " = rectangle(Label(\"",

condition, "\"),

pack(Label(\"n=", sum(root$weights), "\"),

Label(\"", response, "=",

format(root$prediction), "\")),

(", pOffset, ",", -nest, "), lightgray, col(",

root$prediction, "));",

"\ndraw(", id,");\n\n",

child.l, child.r, sep="")

}

treestruct <- plot.node(tree@tree)

cat(file=file.name,

preamble,

"\nimport flowchart;\n",

"pen col(real x) {

real minv = ", minv, ";

real maxv = ", maxv, ";

real ratio = 1 - (x - minv) / (maxv - minv);

return rgb(1, ratio, ratio);

}\n\n",

treestruct, "\n", sep="")

shell(paste("asy -f png", file.name))

}

Each node on the plot contains: the condition leading to it, number of observations and response variable prediction (also intensity of red indicates its relative value).

In order to keep the example simple it is very simplified. Currently handles only regression trees with continuous predictors and will generate errors if variable names contain TeX special characters (like &). Additionally you can control the tree layout only manually by setting variables off.f and off.v or by manipulating picture size in the preamble (and one could write a code to layout the plot automatically).

The code produces png output as I needed this format to show the picture on blog, but of course you can generate eps or pdf file which is probably a more suitable option.

And there is the example of the code use based on standard ctree example:

library(party)

airq <- subset(airquality, !is.na(Ozone))

airct <- ctree(Ozone ~ ., data = airq)

treeAsy(airct, -0.25, 1.4, "tree.asy",

"size(22cm,12cm, keepAspect=false);")

It gives the following output:

Which is much nicer for me in comparison to default plot generated by plot(airct):

An example of OOP in GNU R using S4 Classes

2012-08-18T01:45:00.002-07:00

Recently I have discussed with my friend from WLOG Solutions an implementation of banking cash management engine in GNU R. The code made a nice use of S4 classes so I thought it would be worth showing as an example.

The problem

Every commercial bank needs to provide its customers with access to cash via hundreds of cash access points like ATMs or branches. Bank managers facing this problem have to handle three conflicting objectives: (1) they have to ensure that there is enough cash in cash access points to maintain liquidity, (2) they want to minimize the amount of cash frozen because it is not working for bank and (3) they want to minimize transportation costs from central vault to access points. This is a complex optimization problem which in particular involves the need to predict cash balance in access point every day using historic data.
When designing a solution providing forecasts of cash balances in access points a case for typical application of object oriented approach arises. For each cash point we supply historical data having the same structure and want to obtain a balance prediction. However different access points have different customer usage characteristics and will require different statistical forecasting models. For example in ATM one can only withdraw the money but in branch you can as well make a deposit. Therefore we model the bank using S4 classes.

The implementation

There is a Bank class that can have many CashPoints associated with it. CashPoint is a virtual class that has two implementing classes ATM and Branch. This structure is shown on figure below. Each CashPoint holds its historical balances and has a givePrediction() method that provides a forecast. This method will be implemented differently in ATM and Branch classes.

The example code implementing this structure is given on listing below. First using setGeneric() function we create a generic function givePrediction() that will dispatch appropriate methods following the class of its arguments. Next we create definitions of Bank, ATM and Branch classes of S4 type using setClass() function and create formal method for givePrediction() function for those classes. In our example for ATMs we use linear regression and for Branches simple mean as balance predictors. Notice this method defined for CashPoint class will be invoked if it will not be overridden by appropriate methods in subclasses (it is not possible to create an object of class CashPoint as it is defined virtual).
The code is run by invoking givePrediction() function on a new Bank class instance. The constructor of Bank class reads bank structure data from bank_model.csv file that contains the list of cash points with their ids and types (ATM or Branch). Next it invokes creation of CashPoint-s. Each cash point is initialized with data from branch_balances_data.csv which contains three columns: BranchId, Date, Balance. An appropriate subset of data is first selected using BranchId column. Date and Balance are retained in balances field and contain historical data for this cash point. After creation of an object of type Bank its givePrediction() method is invoked which calls automatically via S4 class system either ATM or Branch givePrediction() method according to the run-time type of cash point.
And here is the code. I hope you will find this simple example of S4 classes application useful.

setGeneric("givePrediction", function(object) {

standardGeneric("givePrediction")

})

setClass("Bank", representation(cashPoints = "list"))

setMethod("initialize", "Bank", function(.Object){

BankModel <- read.table(file = "bank_model.csv",

sep = ";", header = TRUE, stringsAsFactors = FALSE)

.Object@cashPoints <- apply(BankModel, 1, function(cp) {

new(cp[2], cp[1])

})

names(.Object@cashPoints) <- apply(BankModel, 1,

paste, collapse = "_")

return(.Object)

})

setMethod("givePrediction", "Bank", function(object){

return(sapply(object@cashPoints, "givePrediction"))

})

setClass("CashPoint", representation(id = "character",

balances = "data.frame", "VIRTUAL"))

setMethod("initialize", "CashPoint", function(.Object, cashPointId){

.Object@id <-cashPointId

balances <- read.table(file = "branch_balances_data.csv",

sep = ";", header = TRUE)

.Object@balances <- subset(balances,

balances$BranchId == .Object@id, -BranchId)

.Object@balances$Date <- as.Date(.Object@balances$Date)

return(.Object)

})

setMethod("givePrediction", "CashPoint", function(object){

stop("no givePrediction method for this class")

})

setClass("Branch", contains = "CashPoint")

setMethod("givePrediction", "Branch", function(object){

return(mean(object@balances$Balance))

})

setClass("ATM", contains = "CashPoint")

setMethod("givePrediction", "ATM", function(object) {

LM <- lm(Balance ~ as.numeric(Date), data = object@balances)

prediction <- predict(LM,

data.frame(Date = 1 + max(object@balances$Date)))

return(unname(prediction))

})

print(givePrediction(new("Bank")))

To run the code you need to create bank_model.csv and branch_balances_data.csv files. Here you have sample truncated contents (you need to copy-paste the text given below and save in GNU R working directory under appropriate names).

bank_model.csv
CashPointId;CashPointType
CashPoint_1;ATM
CashPoint_2;Branch
CashPoint_3;ATM
CashPoint_4;Branch
CashPoint_5;ATM

branch_balances_data.csv
BranchId;Date;Balance
CashPoint_1;2012-12-01;423000
CashPoint_1;2012-12-02;312000
CashPoint_1;2012-12-03;220000
CashPoint_1;2012-12-04;123000
CashPoint_2;2012-12-01;223000
CashPoint_2;2012-12-02;212000
CashPoint_2;2012-12-03;320000
CashPoint_2;2012-12-04;223000
CashPoint_3;2012-12-01;323000
CashPoint_3;2012-12-02;312000
CashPoint_3;2012-12-03;270000
CashPoint_3;2012-12-04;223000
CashPoint_4;2012-12-01;323000
CashPoint_4;2012-12-02;412000
CashPoint_4;2012-12-03;320000
CashPoint_4;2012-12-04;373000
CashPoint_5;2012-12-01;223000
CashPoint_5;2012-12-02;192000
CashPoint_5;2012-12-03;150000
CashPoint_5;2012-12-04;133000

Animation basics for a vacation

2012-08-05T03:22:00.001-07:00

Since I have a vacation this time I decided to implement some entertaining graphics. I have chosen to animate a Cassini oval.
The task is can be accomplished using polar equation:

The implementation of the animation is given by the following code:

library(animation)

xy <- function(angle, a, b) {

ca <- cos(2 * angle)

r2 <- a ^ 2 * ca + sqrt(a ^ 4 * ca ^ 2 - (a ^ 4 - b ^ 4))

c(sqrt(r2) * cos(angle), sqrt(r2) * sin(angle))

}

go <- function() {

angle <- seq(0, 2 * pi, len = 1000)

par(mar = rep(0,4))

b <- 1 + 0.5 * seq(0, 1, len = 31) ^ 3

b <- c(b, 1.5, rev(b))

for (i in b) {

coord <- t(sapply(angle, xy, a = 1, b = i))

plot(coord, type = "l",

xlim = c(-2, 2), ylim = c(-1.25, 1.25))

polygon(coord, col = "gray")

}

ani.options(interval = 0.1)

saveGIF(go())

Parameter a is fixed to 1 and b changes from 1 to 1.5. In order to achieve smooth animation the sequence defining changes of b is not uniform but is more dense near 1.

And here is the result: