Introduction to Data Science: A Comp-Math-Stat Approach

1MS041, 2021

©2021 Raazesh Sainudiin, Benny Avelin. Attribution 4.0 International (CC BY 4.0)

07. Modular Arithmetic, Linear Congruential Generators, and Pseudo-Random Numbers

What's this all about?

Question:

How can we produce realisations from Uniform$(0,1)$, the fundamental random variable?

i.e., how can we produce samples $(x_1, x_2, \ldots, x_n)$ from $X_1, X_2, \ldots, X_n$ $\overset{IID}{\thicksim}$ Uniform$(0,1)$?

What is SageMath doing when we ask for random()?

Answer:

Modular arithmetic and number theory gives us pseudo-random number generators.

Question:

What can we do with samples from a Uniform$(0,1)$ RV? Why bother?

Answer:

We can use them to sample or simulate from other, more complex, random variables. These simulations can be used to understand real-world phenomenon such as:

See how the modern history of Monte Carlo Methods starts in Los Alamos!

Remark: When mathematical models become too complicated to explain a complex phenomenon of interest (eg. traffic or supply-chain management, predictive maintenance of an expensive component in a national grid/system of power plant, etc.), you have little choice but to resort to modular arithmetic based simulations with pseudo-random numbers.

John von Neumann invented such a pseudo-random number generator called the middle-square method for solving complex problems in the Manhattan Project.

Finally, statistical inference methods you are about to see can be combined with simulation to conduct simulation-intensive statistical inference! These are called Monte Carlo methods and the naming has an interesting folk-tale from Los Alamos.

The starting point for Monte Carlo methods is modular arithmetic ...

Modular arithmetic

Modular arithmetic (arithmetic modulo $m$) is a central theme in number theory and is crucial for generating random numbers from a computer (i.e., machine-implementation of probabilistic objects). Being able to do this is essential for computational statistical experiments and methods that help do this are called Monte Carlo methods. Such computer-generated random numbers are technically called pseudo-random numbers.

In this notebook we are going to learn to add and multiply modulo $m$ (this part of our notebook is adapted from William Stein's SageMath worksheet on Modular Arithmetic for the purposes of linear congruential generators). If you want a more thorough treatment see Modular Arithmetic as displayed below.

Remember when we talked about the modulus operator %? The modulus operator gives the remainder after division:

Arithmetic modulo $m$ is like usual arithmetic, except every time you add or multiply, you also divide by $m$ and return the remainder. For example, working modulo $m=12$, we have:

$$8 + 6 = 14 = 2$$

since $2$ is the remainder of the division of $14$ by $12$.

Think of this as like the hours on a regular analog clock. We already do modular addition on a regular basis when we think about time, for example when answering the question:

Question: If it is 8pm now then what time will it be after I have spent the 6 hours that I am supposed to dedicate this week toward this course?

Answer: 2am.

Modular addition and multiplication with numbers modulo $m$ is well defined, and has the following properties:

Let us make a matrix of results from addition and multiplication modulo 4.

Before that realise that arithmetic over a set $\mathbb{M}$ is merely a function from the Cartesian product $\mathbb{M} \times \mathbb{M}$ of all pairs in $\mathbb{M}$ to $\mathbb{M}$:

$$ \star : \mathbb{M} \times \mathbb{M} \to \mathbb{M} $$

where $\star$ is usually $+$ for "addition" or $*$ for "multiplication".

Thus, we can think of arithmetic as a set of three-tuples: $$(x_1, x_2, x_3) \in \mathbb{M}^3 := \mathbb{M} \times \mathbb{M} \times \mathbb{M}, \quad \text{where} \quad x_3 = x_1 \star x_2. $$

Let us list all the three-tuples for addition modulo 4 next.

Since the arguments to addition modulo m are fixed in a square array $[0,1,2\ldots,m-1]\times[0,1,2\ldots,m-1]$ we can simply focus on the image of the arithmetic modulo $m$ operation and place it over the square array using the following matrix call:

Visualize Modular Arithmetic

In the following interactive image (created by William Stein) we make an addition and multiplication table modulo $m$, where you can control $m$ with the slider. Try changing $m$ and seeing what effect it has:

Look at this for a while. Make $m$ close to 200 slowly and see how the outputs of the two arithmetic operations change.

Answer the following questions (which refer to the default colour map) to be sure you understand the table:

You can change the colour map if you want to.

YouTry

You should be able to understand a bit of what is happening here. See that there are two list comprehensions in there. Take the line matrix(m,m,[(i+j)%m for i in range(m) for j in range(m)]). The list comprehension part is [(i+j)%m for i in range(m) for j in range(m)]. Let's pull it out and have a look at it. We have set the default modulo m to 4 but you could change it if you want to:

This list comprehension doing double duty: remember that a list comprehension is like a short form for a for loop to create a list? This one is like a short form for one for-loop nested within another for-loop.

We could re-create what is going on here by making the same list with two for-loops. This one uses modulo m = 4 again.

Notice that the last statement in the list comprehension, for j in range(m), is the inner loop in the nested for-loop.

The next step that Stein's interactive image is to make a matrix out of the list. We won't be doing matrices in detail in this course (we decided to concentrate on numpy.arrays and multi-dimensional numpy.ndarrays with floating-point numbers for numerical computing instead), but if you are interested, this statement uses the matrix function to make a matrix out of the list named listFormAdd. The dimensions of the matrix are given by the first two arguments to the matrix function: (m, m,...).

Optionally, you can find out more about matrices from the documentation.

YouTry

Try recreating the matrix for multiplication, just as we have just recreated the one for addition.

Modular arithmetic in SageMath

The simplest way to create a number modulo $m$ in Sage is to use the Mod(a,m) command. We illustrate this below.

Let's assign it to a variable so that we can explore it further:

We will compare myModInt to a "normal" SageMath integer:

We can see that myModInt and myInt are different types, but what does this mean? How do they behave?

Try addition:

Was this what you already expected?

What about multiplication?

What's going on here? As we said above, arithmetic modulo mm is like usual arithmetic, except every time you add or multiply, you also divide by mm and return the remainder. 8 x 6 is 48, and the remainder of 48 divided by 12 (12 is our modulo) is 0.

What about this one? What's happening here?

YouTry

Can you create the number giving your year of birth (±1 year) in a similar way. For example, if you are 19 years old now then find the number -19 modulo 2018.

Let us assign 10 modulo 12 to a variable named now:

YouTrys

Add -1 to now:

Put in the expression to do this in SageMath into the cell below:

And subtract 13 from the previous expression.

Put in the expression to do this in SageMath into the cell below:

Also try adding 14 to the previous expression -- the new postition is not shown on this clock but you should be able to see what it should be.

Put in the expression to do this in SageMath into the cell below:

And multiplying 2 by now (or now by 2)

Try making and using some modular integers of your own. Make sure that you can predict the results of simple addition and multiplication operations on them: this will confirm that you understand what's happening.

YouTry Later!

What happens if you try arithmetic with two modular integers? Does it matter if the moduli of both operands a and b in say a+b are the same if a and b are both modular integers? Does this make sense to you?

(end of YouTrys)



Linear Congruential Generators

A linear congruential generator (LCG) is a simple pseudo-random number generator (PRNG) - a simple way of imitating samples or realizations from the Uniform$(0,1)$.

"Pseudo-random" means that the numbers are not really random [From Ancient Greek pseudo means “false, lying”)].

We will look at what we mean by that as we find out about linear congruential generators.

The theory behind LCGs is easy to understand, and they are easily implemented and fast.

To make a LCG we need:

Using these inputs, the LCG generates numbers $x_1, x_2, \ldots x_{n-1} $ where $x_{i+1}$ is calculated from $x_i$ as defined by the following recurrence relation:

$$x_{i+1} \gets mod \left(a x_i + c , m \right)$$

In the above expression $\gets$ denotes "gets" as in variable assignment.

$x_0,x_1,\ldots,x_{n-1}$ is the sequence of pseudo-random numbers called the linear congruential sequence.

We can define a function parameterised by $(m,a,c,x_0,n)$ to give us a linear congruential sequence in the form of a list.

(Remember about function parameters? The function parameters here are m, a, c, x0, and n. Note that the return value of the function is a list.

Implementing the LCG

Let's see how we can define a function for our LCG next.'

You know enough SageMath/Python to understand what every line in the above function linConGen is doing!

The function is merely implementing the pseudocode or algorithm of the linear congruential generator using a for-loop and modular arithmetic: note that the generator produces integers modulo $m$.

Linear Congruential Generators: the Good, the Bad, and the Ugly

Are all linear congruential generators as good as each other? What makes a good LCG?

One desirable property of a LCG is to have the longest possible period.

The period is the length of sequence we can get before we get a repeat. The longest possible period a LCG can have is $m$. Lets look at an example.

You should be able to see the repeating pattern 1, 9, 8, 6, 2, 11, 12, 14. If you can't you can see that the sequence actually contains 8 unique numbers by making a set out of it:

Changing the seed $x_0$ will, at best, make the sequence repeat itself over other numbers but with the same period:

At worst, a different seed specified by x0 might make the sequence get stuck immediately:

An LCG is just a discrete determinstic dynamical system and each seed acts as the initial condition with different behaviours. The behaviour above can be seen as a fixed point at 10 for this LCG initialized at 10.

Thus, the three sequences with different initial conditions or seeds cover all the points in our Ring of integers modulo 17.

What about changing the multiplier $a$?

We want an LCG with a full period of $m$ so that we can use it with any seed and not get stuck at fixed points or short periodic sequences. This is a minimal requirement for simulation purposes that we will employ such sequences for. Is there anyway to know what makes a LCG with a full period?

It turns out that an LCG will have a full period if and only if:

(Different conditions apply when $c=0$. The Proposition and Proof for this are in Knuth, The Art of Computer Programming, vol. 2, section 3.3).

SageMath has a function gcd which can calculate the greatest common divisor of two numbers:

SageMath can also help us to calculate prime factors with the prime_factors function:

$(m, a, c, x0, n) = (256, 137, 123, 13, 256)$ gives a linear congruential sequence with the longest possible period of 256. Let us see how these parameters satisfy the above three requirements while those earlier with $m=17$, $a=2$, and $c=7$ do not.

Thus, the parameters $(m, a, c, x_0, n) = (256, 137, 123, 13, 256)$ do indeed satisfy the three conditions to guarantee the longest possible period of 256.

In contrast, for the LCG example earlier with $m=17$, $a=2$, and $c=7$, although $m$ and $c$ are relatively prime, i.e., $gcd(m,c)=1$, we have the violation that $a-1=2-1=1$ is not divisible by the only prime factor $17$ of $m=17$. Thus, we cannot get a period of maximal length $17$ in that example.

Let us see if the parameters $(m, a, c, x_0, n) = (256, 137, 123, 13, 256)$ that satisfy the three conditions to guarantee the longest possible period of 256 do indeed produce such a sequence:

ourLcg is an example of a good LCG.

The next example will demonstrate what can go wrong.

Consider $m= 256$, $a = 136$, $c = 3$, $x_0 = 0$, $n = 15$.

But, since $a-1=$135 is not divisible by 2, the only prime factor of $m=$256, we get into the fixed point 91 no matter where we start from.

YouTry

See if changing the seed $x_0$ makes a difference?

We can look at the linear congruential sequence generated by $(m,a,c,x_0,n)=(256,137,0,123,256)$ from Knuth's classic [The Art of Computer Programming, vol. 2, Sec. 3.3.4, Table 1, Line 5] and compare it with ourLcg.

Note that although ourLcg has maximal period of 256, lcgKnuth334T1L5 has period of 32. Let us look at them as points.

We can plot each number in the sequence, say ourLcg, against its index in the sequence -- i.e. plot the first number 13 against 0 as the tuple (0, 13), the second number against 1 as the tuple (1, 112), etc.

To do this, we need to make a list of the index values, which is simple using range(256) which as you know will give us a list of numbers from 0 to 255 going up in steps of 1. Then we can zip this list with the sequence itself to make the list of the desired tuples.

Then we plot using points for ourLcg and for lcgKnuth334T1L5.

We can see that in section 3.3.4, Table 1, line 5 in The Art of Computer Programming, Knuth is giving an example of a particularly bad LCG.

When we introducted LCGs, we said that using an LCG was a simple way to imitate independent samples from the Uniform$(0, 1)$ RV, but clearly so far we have been generating sequences of integers. How does that help?

To get a simple pseudo-random Uniform$(0,1)$ generator, we scale the linear congruential sequence over [0, 1]. We can do this by dividing each element by the largest number in the sequence (256 in the case of ourLcg).

Important note: The numbers in the list returned by our linConGen function are integers modulo $m$.

You will need to convert these to Sage integers or Sage multi-precision floating point numbers using int() or RR() command before dividing by $m$. Otherwise you will be doing modular arithmetic when you really want the usual division.

Having sorted that out, we can use a list comprehension as a nice neat way to do our scaling:

This is more like it! We could have a look at this on a plot. Again again want tuples (index, element in scaled sequence at index position), which we can get using range(256) (to get the indexes 0, .., 255) and zip:

Now we have points on the real line. We could use a histogram to look at their distribution. If we are hoping that our LCG, once we have scaled the results, is imitating the Uniform$(0,1)$, what kind of shape would we want our histogram to be?

Is this roughly what you expected?

(Please note: we are using the histogram plots to help you to see the data, but you are not expected to be able to make one yourself.)

We could repeat this for the Knuth bad LCG example:

And show it as a histogram. Given the pattern above, what would you expect the histogram to look like?

Larger LCGs

The above generators are cute but not useful for simulating Lotto draws with 40 outcomes. Minimally, we need to increase the period length with a larger modulus $m$.

But remember that the quality of the pseudo-random numbers obtained from a LCG is extremely sensitive to the choice of $m$, $a$, and $c$.

To illustrate that having a large $m$ alone is not enough we next look at RANDU, an infamous LCG which generates sequences with strong correlations between 3 consecutive points, which can been seen if we manipulate the sequence to make 3-dimensional tuples out of groups of 3 consecutive points.

In the cell below we make our scaled sequence in one step, using a list comprehension which contains the expression to generate the LCG and the scaling.

Have a look at the results as a histogram:

Now we are going to use some of the array techniques we have learned about to resize the sequence from the RANDU LCG to an array with two columns. We can then zip the two columns together to make tuples (just pairs or two-tuples).

Let us resize the LCG to an array with three columns. We can then zip the three columns together to make tuples. The effect will be that if our original sequence was $x_0, x_1, x_2, x_3, \ldots, x_{n-3}, x_{n-2}, x_{n-1}$, we will get a list of triplets, tuples of length three, $(x_0, x_1, x_2), (x_3, x_4, x_5), \ldots, (x_{n-1-2}, x_{n-1-1}, x_{n-1})$. Unlike the pairs in 2D which seem well-scattered and random, triplets from the RANDU LCG are not very random at all! They all lie on parallel planes in 3D.

You can alter your perspective on this image using the mouse. From a particular perspective you can see that something has gone horribly wrong ... RANDU is a really ugly LCG.

The above generators are of low quality for producing pseudo-random numbers to drive statistical simulations. We end with a positive note with a LCG that is in use in the Gnu Compiler Collection. It does not have obvious problems as in small periods or as high a correlation as RANDU.

Even good LCG are not suited for realistic statistical simulation problems.

This is because of the strong correlation between successive numbers in the sequence. For instance, if an LCG is used to choose points in an n-dimensional space, the points will lie on, at most, $m^{1/n}$ hyper-planes. There are various statistical tests that one can use to test the quality of a pseudo-random number generator. For example, the spectral test checks if the points are not on a few hyper-planes. Of course, the Sample Mean, Sample Variance, etc. should be as expected. Let us check those quickly:

Recall that the population mean for a Uniform$(0,1)$ RV is $\frac{1}{2}$ and the population variance is $\frac{1}{12}$.

To go into this topic of how good a pseudo-random number generaror is and what "good" means in many different mathematical senses (in detail) is clearly beyond the scope of this course. You should just remember that using computers to generate pseudo-random numbers is not a trivial problem and use care when employing them especially in higher dimensional or less smooth problems. The mathematics behind this has one of the most beautiful roots of the subject called number theory.

More Sophisticated Pseudo-Random Number Generators

We will use a pseudo-random number generator (PRNG) called the Mersenne Twister for simulation purposes in this course. It is based on more sophisticated theory than that of LCG but the basic principles of recurrence relations are the same.

(The Mersenne Twister is a variant of the recursive relation known as a twisted generalised feedback register. See Makato Matsumoto and Takuji Nishimura, "Mersenne Twister: A 623-dimensionally equidistributed uniform pseudo-random number generator, ACM Transactions on Modelling and Computer Simulation, vol. 8, no. 1, Jan. 1998, pp. 3-20. if you have ACM access or via links at Mersenne Twister.

The Mersenne Twister has a period of $2^{19937}-1 \approx 10^{6000}$ (which is essentially a Very Big Number) and is currently widely used by researchers interested in statistical simulation.

Simulating a Drunkard's walk

We start with some warm up

Some Warmup before Simulating a Drunkard's Walk

Accumulating sequences with numpy.cumsum

Our example using prngs is going to use a very useful function from the numpy module called cumsum. This calculates a cumulative sum from a sequence of numeric values. What do we mean by a cumulative sum? If we have a sequence $x_0, x_1, x_2, x_3, \ldots, x_n$, then we can derive a sequence that is the cumulative sum,

$$x_0, \,x_0 + x_1,\, x_0+ x_1+x_2, \,\dots, \sum_{i=0}^nx_i$$

Try evaluating the next cell to get the cumulative sum of the sequence 0, 1, 2, ..., 9

Since pylab is a superset of numpy, we can also do the following imports instead:

You will see that the result is in a pylab.array or numpy.array. This can be useful, but all we need is a list, we can convert this back to a list with the array's tolist method:

Using the list function to make a list will do the same thing:

Now, say we have some probabilities in a list and we want to get the cumulative sum of them.

We can use cumsum as we did above. We are going to start with just two probabilties which (for reasons that will shortly become clear), we will assign as values to variables pLeft and pRight.

Note that cumsum works on a tuple as well as a list, providing that the elements of the tuple are some type of number.

Simulating a Drunkard's Walk

You are now going to use some simulated Uniform$(0,1)$ samples to simulate a Drunkard's Walk.

The idea of the Drunkard's Walk is that the Drunkard has no idea where s/he is going: at each decision point s/he makes a random decision about which way to go. We simulate this random decison using our Uniform$(0,1)$ pseudo-random samples.

We are going to have quite a limited version of this. At each decision point the Drunkard can either go one unit left, right, up or down (pretend we are in the middle of Manhattan with its streets and avenues).

Effectively, s/he is moving around on a $(x, y)$ coordinate grid. The points on the grid will be tuples. Each tuple will have two elements and will represent a point in 2-d space. For example, $(0,0)$ is a tuple which we could use as a starting point.

First, we recall useful feature of tuples: we can unpack a tuple and assign its elements as values to specified variables:

We use this useful feature in the functions we have defined below.

You now know what is happening in each step of these functions:

and should be able to understand what is happening from the detailed comment and doc-string (read!).

You should be thinking "Hey that is a lot of duplicated code: I thought that functions should help us not to duplicate code!". You are completely right, but this example will be easier to understand if we just live with some bad programming for the sake of clarity.

Now, experience the flexibility of Python: We can have lists of functions!

To demonstrate how we can use our list of functions, take some path which is represented by a list of points (tuples):

If you want to, see if you can find the tuple which is the last point in the path so far using len(somePath) and the indexing operator [ ].

We put the functions into the list in order, left to right, so we know that the first function in the list (index 0) is makeLeftTurnPoint, and the second function in the list (index 1) is makeRightTurnPoint.

We can now call makeLeftTurnPoint by calling movementFunctions[0], and passing it the arguement somePath which is our example path so far. We should get back a new point (tuple) which is the point to the immediate left of the last point in the list somePath:

This has not added the new point to the list, but we can do this as well:

YouTry

Now try doing a similar thing to find the next right point.

That's all very well, but what about some random moves? What we are now going to do is to use random() to make our decisions for us. We know that the number generated by random will be in the interval [0,1). If all directions (up, down, left, right) are equally probable, then each has probability 0.25. All directions are independent. The cumulative probabiltiies can be thought of as representing the probabilities of (up, up or down, up or down or left, up or down or left or right).

Using these accumulated probabilities we can simulate a random decision: if a realisation of a Uniform$(0,1)$, $u$, is such that $0 \le u < 0.25$ then we go left. If $0.25\le u < 0.50$ we go right, if $0.50 \le u < 0.75$, we go up, and if $0.75 \le u < 1$ we go down.

We can demonstrate this:

You will see that we have only dealt with the cases for left and right. You may well not get $n$ lines of output. You can have a go at improving this in very soon. First, one more thing ...

We can tidy up the if statement a bit by using another for loop and the break statement. The break statement can be used to break out a loop. In this case we want to compare u to each cumulative probability until we find a 'match' (u < cumulative probability) and then break out of the for loop.

YouTry

Now, try adding the cases to deal with up and down.

Now we can combine all this together to make a simulated Drunkard's Walk: First, a little helper function to plot a list of points as lines.

Now the real stuff by putting it all together:

A bit boring? A bit one-dimensional? See if you can add up and down to the Drunkard's repetoire. You will need to:

Try it and see!

This is a warmup for an assignment - you can do and show me that you comprehend it.

Simple Random Walks

Simple random walks on $\mathbb{Z}^d$ when scaled appropriately leads to Browninan motion that can be formalized as the Weiner process. This is a primary ingredient in stochastic calculus.

Simulating An Independent Rademacher Random Walk in $\mathbb{Z}$ over Discrete Time in $\mathbb{Z}_+$

Here discrete time is indexed by the non-negative intereger denoted by $\mathbb{Z}_+ := \{0,1,2, 3,\ldots\} = \{0\} \cup \mathbb{N}$.

The following @interact-ive function Independent1DRademacherRandomWalkInDiscreteTime has been adapted from a similar function by William Stein, the main founder of SageMath. William has left academia (Math Dept. at Univ. of Washington) and started a company called COCALC: Collaborative Calculation in the Cloud at https://cocalc.com around 2017 (personal communication).

Here we can develop intuition for an independent Rademacher random walk in $\mathbb{Z}^1$ (also known as simple symmetric random walk) over multiple dsicrete time-steps.

This can be seen to be the Weiner process in the limit (as the number of discrete steps approach infinity under appropriate scaling). Such processes are at the heart of stochastic calculus - an important toll in financial mathematics, stochastic control, statistical time series mdoels, etc.

But first we need to learn about an important Random Variable known as the Rademacher Random Variable, say $K$, whose probability mass function is: $$f(k) = \mathbf{1}_{\{-1,+1\}}(k)\frac{1}{2} = \begin{cases} \frac{1}{2} \quad \text{if } k=-1\\ \frac{1}{2} \quad \text{if } k=+1\\ 0 \quad \text{otherwise.}\end{cases}$$

Simulating An Independent Rademacher Random Walk in $\mathbb{Z}^3$ over Discrete Time in $\mathbb{Z}_+$

The following random walk should be natural now. Such simulations in 3D are the building blocks for modelling diffusion of airborne particles in the space (their extensions with drift from weather models are often used to predict pollutant flows in the atmosphere across state or national boundaries, for example).

Simulating Multiple Independent Random Walkers in 2D-space over Time

Such simulations are the starting point for traffic engineering simulations used to optimise city-planning of transportation networks. These random walkers can be made dependent and also made to move on a road network with right speeds, etc. Here we look at the vanilla example.

A Comment-it YouTry!

Try to comprehend this code and add comments to each line. This is mostly an opportunity for you to see for yourself what you should be able to already do, in terms of simulation.

Example of Simulation-Intensive Inference

Suppose you want to know the probability of a scaled random walk visiting a particular region after having started from an initial position, or say the expected time to visit, etc., then we can resort to simulations to get the estimate of interest.

Soon, we will see why this approach is justified mathematically via the Law of Large Numbers and the Central Limit Theorem. For now, let's build more intuition.

Animations of your Independent Simulation Trials

In summary, we can see from the animations of simulations above, we can multiply the independent estimates of the probabilities of going:

Think! Why can we multiply these probability estimates? What property of the random walk justifies this?

We know that the estimates themselves are simply simulations of very definition of probability as the expectation of the event, i.e., the sample mean of the indicator function of the event - in other words the fraction of simulation trails that actually resulted in the event of interest.

Imagine how more complex rare events (these are just events with very low probability) can be estimated by simulating various independent componenets of them, if they exist.

simulations based estimate of p = probrobability of our random walk going from first circle to second to third (the walks are shown up to few steps before they reach the destination circle)
1 of 10 trials go from first circle to second 1 of 10 trials go from second circle to third simulation-based estimate of p from 10 X 10 IID trials
10 simulations from first circle to second 10 simulations from second circle to third prob estimates from 10 X 10 IID simulations to estimate p

How to do it yourself if you want to - Some First Steps

Sometimes you may want to simulate a phenomenon and make an animation and save it as an animated GIF file, either to build intuition ro to make a presentation to a broader audience. Here the following code snippets can get you started.

This is advanced content and you need to know more about administering your system. It is here for those who want to take it further on their own. To understand the html code in the animation above go to w3 schools for table - this is how people keep learning things they don't know about coding certain things.

Also, your mathematical intuition can be strengthened in this way when done systematically using mathematical libraries in environments like SageMath.

You Try

If you have a SageMath installation, done carefully on Mac OS X or Linux OS or in Ubuntu 20.10 within WSL2 (Windows subsystem for Linux)by installing the needed additional packages at appropriate levels in your system, then you can make your own animated gif in your system, independent of the three OSs. This will be an Assignment.