Tips for Writing Computer Simulations of Business Processes

Written by Peter Rosenmai on 20 Dec 2019.

Here's a talk I gave in 2019 in Birmingham, Alabama, on creating computer simulations of business processes. It was based on work I had been doing writing Python simulations.

Some of the more interesting sections:
22:50: The sorts of graphs I produce to understand simulated processes over time.
37:54: A spreadsheet I used as the input to a simulation engine that I wrote in Python.
55:18: An example of the animations produced by one of my simulations.

Generating Random Survival Times From Any Hazard Function

Written by Peter Rosenmai on 14 Apr 2017.

Let's get 1,000 random survival times (for use, perhaps, in a simulation) from a constant hazard function (hazard = 0.001):

And let's check that the Kaplan-Meier curve for these survival times appproximates, as expected, the curve P(t) = exp(-0.001t):

require("survival")
plot(survfit(Surv(survival_times, rep(1, length(survival_times)))~1), xlab="Time", ylab="Survival Probability",
     main="Sampled and Expected Survival Curves for h(t) = 0.001")
neg_exp_fn = function(x){exp(-0.001 * x)}
curve(expr=neg_exp_fn, from=0, to=max(survival_times), add=TRUE, col="red", lwd=2)

Atlanta Crime Map

Written by Peter Rosenmai on 27 Dec 2016. Last revised 14 Apr 2017.

This interactive map is based on crime data provided by the Atlanta Police Department for 2015 and 2016. That dataset is subject to considerable change over time, as crimes are often reported months after they actually occurred. I last downloaded and incorporated that data into this map on 14 April 2017.

Note that the date shown here for a crime is the date on which it occurred, not the date on which it was reported.

I built this map using Leaflet, Leaflet.Markercluster and jQuery Date Range Picker Plugin.

Fitting an Exponential Curve to a Stepwise Survival Curve

Written by Peter Rosenmai on 27 Aug 2016. Last revised 13 Mar 2017.

Let's fit a function of the form f(t) = exp(λt) to a stepwise survival curve (e.g. a Kaplan Meier curve).

Here's the stepwise survival curve we'll be using in this demonstration:

# The data that describes the stepwise survival curve.
surv_time <- c(0.2,  0.5, 1.2,   2,   5,  6.5,    9,   15, 17.3,   35,   80,  200)
surv_prob <- c(  1, 0.94, 0.9, 0.8, 0.6, 0.57, 0.41,  0.2, 0.18, 0.07, 0.02, 0.01)

# Plot the curve.
plot_stepwise_surv_fn <- function(time, prob, col="black"){
   time <- c(0, rep(time, each=2))
   prob <- c(1, 1, rep(prob, each=2)[1:(2 * length(prob) - 1)])
   plot(x=time, y=prob, ylim=c(0, 1), type="l", col=col, xlab="Time", ylab="Probability")
}
plot_stepwise_surv_fn(surv_time, surv_prob)

Calculating Kaplan Meier Survival Curves and Their Confidence Intervals in SQL Server

Written by Peter Rosenmai on 1 Jan 2016.

I provide here a SQL Server script to calculate Kaplan Meier survival curves and their confidence intervals (plain, log and log-log) for time-to-event data.

Example 1: Customer Attrition, Ungrouped, Without Censoring

Suppose a web-application company has seen its ten customers cancel their subscriptions after 0.5, 1, 10, 10, 11, 13.5, 14, 19, 19.5 and 30 months (from the start of their respective subscriptions). Based on this data, we want to estimate the probability of a new customer remaining a customer more than, say, 12 months—and we want a confidence interval around that estimate.

We create the input table:

Calling SPSS Modeler from R

Written by Peter Rosenmai on 12 Dec 2015.

SPSS Modeler streams can be executed from R via input files and command-line calls. It's a hacky technique, but it works. And it's an easy way to make use of Modeler's excellent Expert Modeler functionality.

First, create an example of the data file that you want Modeler to read in. For example, if you want to run Modeler on a single time series, your data file will probably be a text file comprising a date column and a value column.

Then create and save a Modeler stream that reads in that file, fits the required model and produces an output file. This is fairly easy so I won't cover it here. But this is how it might look:

Kaplan Meier Survival Curve Grapher

Written by Peter Rosenmai on 13 Jan 2015.

Beta Distribution PDF Grapher

Written by Peter Rosenmai on 1 Jan 2015.

Here's a D3-rendered graph of the probability density function (PDF) of the beta distribution. Move the sliders to change the shape parameters or the scale of the y-axis.

Installing Rpy2

Written by Peter Rosenmai on 16 Nov 2014.

Here's how I installed the rpy2 module for communicating between R and Python. I did this for R version 3.0.2, Python version 3.4.1, and rpy2 version 2.4.4 on a 64-bit machine running Windows 7.

Setting System Variables

First, I got the full pathname of my R executible by right-clicking the R icon in my Start menu and selecting Properties. It was: C:\Program Files\R\R-3.0.2\bin\x64\Rgui.exe. And it only took a moment of poking around to find the full pathname of my Python executible: C:\Anaconda3\python.exe.

Time to look at my system variables. I right-clicked on Computer in my Start menu and selected Properties; I then clicked on Advanced System Settings in the window that appeared. In the Advanced tab of the System Properties window, I clicked the Environment Variables button. The lower half of the resulting Environment Variables window showed my system variables.

The first system variable I had to deal with was Path. It didn't include the directory in which my R executable sits, so I added it: C:\Program Files\R\R-3.0.2\bin\x64\.

I then added an R_HOME system variable and set it to the top level directory of my R installation: C:\Program Files\R\R-3.0.2\.

And I added an R_USER system variable and set it to the directory that the rpy2 module would install into: C:\Anaconda3\Lib\site-packages\rpy2\.

Estimating the Distance Between GPS Points While Accounting for Circular Error Probable (CEP)

Written by Peter Rosenmai on 27 Sep 2014. Last revised 12 Oct 2014.

Given two GPS points recorded as being d metres apart with circular error probable (CEP) of c₁ and c₂ metres respectively, the true distance between the recorded points has the distribution

where:

• C is a non-central chi-square random variable having 2 degrees of freedom and non-centrality parameter d² / σ²
• 
• 

(I give a proof of this easy result below.)

Generating Random Numbers from Any Non-Negative Function

Written by Peter Rosenmai on 2 Aug 2014.

Here's some R code that generates random numbers from the probability distribution described by a given non-negative function. This can be useful when running simulations or generating datasets for testing purposes.

For example, suppose you want to generate random data from a distribution that looks something like this back-of-an-envelope sketch:

Dragging Outliers Towards the Mean or Median

Written by Peter Rosenmai on 21 Jun 2014.

Here's an example of how to use R to smoothly drag towards the mean outliers that are more than a given number of standard deviations from the mean—or median absolute deviations from the median, or whatever—so that the most extreme outliers are dragged in the most. I don't really agree with mangling data in this way and I think the task is a trivially simple one, but I've often been asked how to do it… so here's how you might go about it.

First, for demonstation purposes, I create a dataset with some obvious outliers:

I drag the outliers towards the mean using the standard deviation:

Graphing Survival and Hazard Functions

Written by Peter Rosenmai on 11 Apr 2014. Last revised 13 Jun 2015.

Here's some R code to graph the basic survival-analysis functions—s(t), S(t), f(t), F(t), h(t) or H(t)—derived from any of their definitions.

For example:

# Define time points at which to plot
t <- seq(0, 4, 0.01)

# Define h(t), the hazard function
hazard_fn <- function(t){2*t}

# Graph S(t), the survival function derived from h(t)
y <- apply_survival_function(t, hazard_fn, supplied_fn_type="h", fn_type_to_apply="S")
plot(x=t, y=y, xlim=c(0, max(t)), ylim=c(0, max(y)), main="S(t)", ylab="Survival Probability", type="l")

Calculating Influence Matrices

Written by Peter Rosenmai on 22 Feb 2014.

I present here R code to calculate the influence that a set of points have upon each other, where influence is a function of (1) the distance between the points and (2) the inherent influence of the points.

Example: Light Bulbs Shining on Light Bulbs

Consider, for example, five light bulbs with brightness given by this vector:

Now, suppose that the distance between the light bulbs (in metres) is given by this distance matrix:

This matrix tells us, for instance, that bulbs two and three are 12 metres apart. Note that the distance matrix is symmetrical about a zero diagonal. This corresponds with the ordinary notion of distance: Any point is a zero distance from itself, and the distance from point A to point B equals the distance from point B to point A. However, my code permits non-symmetric distances: If bulb two is "uphill" from bulb three, [2, 3] will be greater than [2, 3].

Calculating a Distance Matrix for Geographic Points Using R

Written by Peter Rosenmai on 30 Jan 2014.

Here's an example of how to calculate a distance matrix for geographic points (expressed as decimal latitudes and longitudes) using R:

For example, the above distance matrix shows that the straight-line distance—accounting for curvature of the earth—between Los Angeles and NYC is 3,945 km.

R's Crazy If-Else Parser Quirk

Written by Peter Rosenmai on 17 Jan 2014.

I present here what I consider to be a fiendishly weird quirk in R's code parser.

Okay, so what do you expect the following code to do?

Run it and you'll see that it prints the number 1, as you would expect.

Okay, now what happens when you remove the top-level if block?

Debugging with the Superassignment Operator

Written by Peter Rosenmai on 31 Dec 2013.

I show here how the <<- assignment operator may be used to debug R functions by writing local variables into the global environment.

Consider the following code:

Running that code produces a graph of sunspot activity since 1950 and an exponential smoother of those data. But there's a problem: The graph subtitle doesn't come out properly.

R's Flavours of Stacked Dot Plots

Written by Peter Rosenmai on 25 Nov 2013. Last revised 13 Jan 2014.

The humble stacked dot plot is, I think, often preferable to the histogram as a means of graphing distributions of small data sets. To gauge how closely a histogram approximates an underlying population distribution, one must take into account the number of points that the histogram is based on (the sample size). Many readers fail to do this—and all too often the sample size is not provided within the graph. However, a dot plot lets any reader make an immediate guess at how closely the graph follows the shape of the underlying distribution.

Several R functions implement stacked dot plots. But which one to use?

There's this one from the base graphics package:

stripchart(faithful$waiting, method="stack", offset=0.5, pch=1)

Creating Datasets Interactively

Written by Peter Rosenmai on 25 Nov 2013. Last revised 1 Jan 2014.

The following Create2DimData() R function allows two-dimensional datasets (e.g. of heights and weights) to be created by clicking with a mouse within a plot. This can be useful if you need to throw together a dataset for demonstration purposes.

For example, try calling Create2DimData() like this:

df.points <- Create2DimData(xlim=c(0,10), ylim=c(0,5))

Using Mahalanobis Distance to Find Outliers

Written by Peter Rosenmai on 25 Nov 2013. Last revised 30 Nov 2013.

R's mahalanobis() function provides a simple means of detecting outliers in multidimensional data.

For example, suppose you have a dataframe of heights and weights:

When plotting these data (generated for this example using an interactive plot), you could mark as outliers those points that are, for instance, more than two (sample) standard deviations from the mean height or mean weight:

Using the Median Absolute Deviation to Find Outliers

Written by Peter Rosenmai on 25 Nov 2013. Last revised 13 Jan 2013.

One of the commonest ways of finding outliers in one-dimensional data is to mark as a potential outlier any point that is more than two standard deviations, say, from the mean (I am referring to sample means and standard deviations here and in what follows). But the presence of outliers is likely to have a strong effect on the mean and the standard deviation, making this technique unreliable.

For an example of this well-known problem, try running the following R code:

Exploring the World Bank's Gini Index Data with R

Written by Peter Rosenmai on 17 Dec 2013. Last revised 18 Dec 2013.

Let's have a look at the Gini index data available from the World Bank through R's WDI package.

For those who haven't met it before, the Gini index is an elegantly constructed measure of, typically, income inequality. A Gini index of 0 represents a perfectly equal economy; a Gini index of 100 represents a perfectly unequal economy. (To find out more about the Gini index, have a look at my Gini index calculator.)

Let's search for the Gini index within the World Bank's datasets:

R Code to Remove Duplicates from a SQL Server Database

Written by Peter Rosenmai on 8 Dec 2013.

It's easy to remove duplicate rows from an R dataframe using the unique() function:

But this can be slow for large dataframes. Hideously slow, even.