Linear systems are systems that have predictable outputs when there are small changes in the inputs to the system. Nonlinear systems are those that produce disproportionate results for proportional changes in the inputs. Both linear and non-linear systems are common enough in nature and industrial processes, or more accurately, many industrial and natural processes can actually be modeled as standard linear or nonlinear systems.

## Nonlinear Dynamics and Chaos

Nonlinear dynamical systems are essentially systems that exhibit time-dependent behaviour and in a non-linear manner. A special class of such systems also exhibit chaos, which is defined as sensitive dependence upon initial conditions. There are great textbooks available on the subject, by researchers such as Steven Strogatz (Cornell University, Ithaca, New York).

While R is often used to run statistical analysis and studies of various kinds including advanced machine learning algorithms, it isn’t often used to study systems like this. However, they are easily modeled in R, and like any programming language, the surfeit of functions can help us understand statistical aspects of the behavior of such systems. There’s extensive material available on the internet, and Steven Strogatz’s lectures on Nonlinear Dynamics and Chaos provide a very deep treatment of the subject.

## Logistic Maps and Bifurcation Diagrams

A logistic map is a function that describes how the state of a particular dynamical system evolves from one point in time to the next. It allows us to understand bifurcations, and understand what kinds of conditions produce sensitive dependence on initial conditions. More on bifurcation diagrams here.

## Nonlinear Dynamical System in R

The system I’ll describe here is a probabilistic system that is based on the binomial distribution’s mechanics. This distribution is used to model events with two probabilities (success or failure), of some probability. A special case of this is the coin toss, the Bernoulli distribution.

In our example, we’re trying to understand the probability of success in a repetitive or sequential, identical game with two outcomes, provided we know the initial chance of success. In doing so, we’re exploring the impact of the number of games we’re playing in sequence, and the number of wins or losses we expect in each case. The end result from this, is a matrix, which we call * problemset * in this specific case.

## Animations in R

The R package “animation” has functions which can enable sequential graphics (such as that generated within a loop) to be saved as a GIF animation file. This is especially handy when we’re trying to understand the impact of parameters, or when we’re trying to illustrate the data, analysis and results in our work in some sequence. We’ll use the saveGIF() function here to do just such a thing – to save a sequence of images of logistic maps in succession, into a single GIF file.

library(animation) #Set delay between frames when replaying ani.options(interval=.05) #Do our plots within the saveGIF command parantheses, in order to capture the matrix plots we're generating saveGIF({ for (inval in seq(0,1,length.out = 200)){ pfirst <- inval #Defining a function to calculate event probability based on starting probability assumptions prob <- function(game){ n <-game[1]; k <-game[2]; p <-game[length(game)]; return( factorial(n) / (factorial(n-k) * factorial(k)) * p^k * (1-p)^(n-k) ); } iter <-100 k <- 2 games <- seq(2,100,1) victories <- rep(5,length(games)) problemset <- cbind(games, victories, rep(pfirst, length(games))) #Setting up a temporary variable to store the probability values per iteration out<-NULL for (i in seq(1,iter,1)){ for (i in seq(1,length(problemset[,1]), 1)){ out<-rbind(out,prob(problemset[i,])) } problemset <-cbind(problemset,out) out<-NULL } #Using the matrix plot function matplot() to plot the various columns in our result matrix, together matplot(problemset[,seq(3,length(problemset[1,]), 1)], type = "l", lwd = 1, lty = "solid", main = paste("Logistic Map with initial probability = ",round(pfirst,2)), ylab = "Probability", xlab = "Number of games", ylim = c(0,0.5) ) } })

The code above generates an animation (GIF file) in your default R working directory. It allows us to examine how different system parameters could affect the probability of events we’re evaluating in the sample space we have in this problem. Naturally, this changes depending on the number of games you indicate in the * games* variable in the code. The GIF files are shown at the end of this section – two different cases are shown.

Here’s a logistic map generated by a slightly modified version of the code above. This shows the calculated probabilities for different combinations of games, and won games, based on initial assumed win percentages. The initial assumed win percentage in this case is 0.1 or 10%, the horizontal line running through the graph.

A number of systems can only be described well when we see their performance or behaviour in the time domain. Visualizing one of the parameters of any model we construct along the time domain therefore becomes useful to try and spot patterns and trends. Good data visualization, for simple and complex problems alike, isn’t only about static images and infographics, but dynamic displays and dashboards of data that change and show us the changing nature of the system being modeled or data being collected. Naturally, this is of high value when we’re putting together a data science practice in an organization.

## Concluding Remarks

We’ve seen how animated logistic maps are generated in the example here, but equally, the animated plots could be other systems which we are trying to describe in the time domain, or systems we’re trying to describe in an elaborate way, especially when many parameters are involved. Linear and nonlinear systems aside, we also have to deal with nonlinear and dynamical systems in real life. Good examples are weather systems, stock markets, and the behaviours of many complex systems which have many variables, many interactions, although simple rules. Logistic maps can tell us, based on simplified representations of our data, about the regimes of chaos and order in the behaviour of the system, and what to expect. In this sense, they can be used in cases where it is known that we’re dealing with nonlinear and dynamical systems.

Pingback: 2 – Animated Logistic Maps of Chaotic Systems in R | Exploding Ads

Pingback: Found on WordPress.com: Animated Logistic Maps of Chaotic Systems in R | Simon Ouderkirk

Interesting article!Thank you for the writeup.

Logistic maps were unfamiliar to me, so I tried as best as I could to understand your code and use Wikipedia to get a handle on what was going on. If I understand correctly, I think that the label of your horizontal (x) axis is misleading. It isn’t the “number of games”, but the “number of recursive games”.

For example, (using x to represent “number of games” and y to represent “probability”) when x=1, then y represents the probability of exactly 5 heads in n flips, where each line represents a different n, from 5 to 100. However, when x=2, then y represents the probability that exactly 5 of n sets of n flips has 5 heads. I think it is the recursive nature of the problem that you are missing with the description “number of games”.

I found your code a little hard to follow. Your prob function is the dbinom function with the order of arguments changed. dbinom (and your own prob function) are vectorized, so you need not iterate through the values of problemset. Also, I believe there was a bug in the code where 2 to 4 flips were being considered, even though 5 victories were necessary. Below, I supply a more concise version of your code.

library(animation)

ani.options(interval=.05)

plot.set <- function(pfirst, iter = 100, victories = 5) {

games <- victories:100

problemset <- matrix(rep(pfirst, length(games)))

for (i in seq(iter))

problemset <-cbind(problemset,

dbinom(victories, games, problemset[,ncol(problemset)]))

matplot(problemset,

type = "l", lwd = 1, lty = "solid",

main = paste("Logistic Map with initial probability =",round(pfirst,2)),

ylab = "Probability", xlab = "Number of games", ylim = c(0,0.5))

}

saveGIF(for(pfirst in seq(0,1,length.out = 200)) plot.set(pfirst))

LikeLiked by 1 person

Aman, thanks for taking interest in this blog and for your informative comment, appreciate your time.

You’re right about “number of recursive games” being a better label for the x axis.

I was indeed using the binomial distribution, so the dbinom() function would have been a good fit. Thanks for pointing out the bug too.

As someone who is still picking up fluency in R (this blog is also helping me do that), I found your code far simpler to understand than my own! Thanks for posting it here, really appreciate it. I haven’t run it yet since my work laptop doesn’t have ImageMagick at the moment, but I will run the code when I have the time later on.

LikeLike