text:
# Simple, fast and accurate evaluation of the action of the exponential of a rate matrix on a probability vector

Given a time-homogeneous, finite-statespace Markov chain with a rate matrix, or infinitesimal generator of Q, an initial distribution of ν and a time interval of t, the distribution across states at time t, ν^ (t):=ν^ (Q t) is typically required for one or many time intervals, t, either as a contribution to the likelihood to enable statistical inference, or to understand the evolution of the process itself. When, as in many statistical applications in, for example, epidemics, systems biology and ecology, Q arises from a reaction network, or when Q corresponds to a negative graph Laplacian, it is usually sparse. Building on a relatively recent development for the action of the exponential of a general sparse matrix on a vector, we use the special properties of rate matrices to create the Single Positive Series method, which accurately and quickly evaluates the distribution across states at a single time point in O((r+1)ρ d) operations, where d is the dimension of the statespace, ρ=_i=1,...,d|Q_ii| and r is the average number of positive entries in the rows of Q. We also present the Multiple-Use Single Positive Series algorithm which simultaneously evaluates the distribution vector at k different time points in O((k+r)ρ d) operations. We demonstrate across a range of examples that the Single Positive Series method is both more accurate and more than an order of magnitude more efficient than the nearest generic and generally available competitor, and that the Multiple-Use Single Positive Series algorithm improves upon this still further.

## Authors

• 14 publications
• ### Convergence of Contrastive Divergence Algorithm in Exponential Family

This paper studies the convergence properties of contrastive divergence ...
03/17/2016 ∙ by Tung-yu Wu, et al. ∙ 0

• ### Efficient construction of an HSS preconditioner for symmetric positive definite ℋ^2 matrices

In an iterative approach for solving linear systems with ill-conditioned...
11/15/2020 ∙ by Xin Xing, et al. ∙ 0

• ### A Monte Carlo method for computing the action of a matrix exponential on a vector

A Monte Carlo method for computing the action of a matrix exponential fo...
04/29/2019 ∙ by Juan A. Acebron, et al. ∙ 0

• ### Variable and Fixed Interval Exponential Smoothing

Exponential smoothers are a simple and memory efficient way to compute r...
02/11/2015 ∙ by Javier R. Movellan, et al. ∙ 0

• ### On the Computational Complexity of Linear Discrepancy

Many problems in computer science and applied mathematics require roundi...
07/31/2020 ∙ by Lily Li, et al. ∙ 0

• ### Estimating entropy rate from censored symbolic time series: a test for time-irreversibility

In this work we introduce a method for estimating entropy rate and entro...
09/23/2020 ∙ by Raul Salgado-Garcia, et al. ∙ 0

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

For many applications, statistical methods require the exponentiation of large, sparse, matrices. Such calculations are needed, for example, to evaluate the likelihood for continuous-time Markov models

(e.g. Andersson and Britton, 2000) or to evaluate certain metrics for networks (e.g. Hammond et al., 2013). Standard algorithms for evaluating the matrix exponential (such as expm from the Matrix package in R) scale poorly with dimension, . However, in many applications interest lies only in the action of the matrix exponential on a vector; furthermore, when the matrix to be exponentiated is sparse, the calculation can be performed in operations. Recently, Al-Mohy and Higham (2011) introduced an algorithm which in the vast majority of examples considered was both faster and more accurate than its competitors, but this has largely been ignored in the statistics literature. Here we present further developments of the algorithm of Al-Mohy and Higham (2011) which take advantage of the special structure of matrix exponential calculations that are common in statistics. Our key algorithm is more straightforward to implement than that of Al-Mohy and Higham (2011), yet faster, and at least as accurate, as well as being an order of magnitude faster than its nearest widely-available competitor, making accurate likelihood-based inference feasible and straightforward for systems with statespace sizes in the tens of thousands.

The exponential of a square matrix, is defined via its infinite series: . As might be anticipated from the definition, algorithms for evaluating take operations see Moler and Van Loan (2003) for a review of many such methods. However, for a -vector, , the product is the solution to the initial value problem , , and is the key component of the solution to more complex differential equations such as . For this reason the numerical evaluation of the action of a matrix exponential on a vector has received considerable attention of itself (e.g. Gallopoulos and Saad, 1992; Saad, 1992; Sidje, 1998; Al-Mohy and Higham, 2011).

When is dense, evaluating

 eMv=∞∑i=01i!Miv (1)

takes operations. However, motivated by the examples we detail in Section 1.1 our interest lies in large sparse matrices, and the number of operations then reduces to , where is the average number of entries in each row of . In this case, as exemplified in the popular Expokit FORTRAN routines Sidje (1998)

, it is usual to estimate

via its projection on to the Krylov subspace , where .

With double-precision arithmetic, real numbers are stored to an accuracy of approximately . Thus, when a scalar, , is very large and negative, naive numerical evaluation of the series for can be innaccurate due to the cancellation of very large positive and negative terms. Evaluation of is, typically, similarly affected, and to control these errors, algorithms for evaluating (e.g. Sidje, 1998) use the identity

 eMv=(eM/s)sv=eM/seM/s…eM/sv, (2)

where the integer is chosen so that is small enough (for a suitable choice of norm) that cancellation of positive and negative terms which are ‘too large’ does not occur. The calculation on the right of (2) typically involves many more numerical operations than the direct calculation on the left and so should be the smallest integer that mitigates against the cancellation of large positive and negative terms, so as to minimise the accumulation of rounding errors.

The relatively recent, ground-breaking paper Al-Mohy and Higham (2011) presents a new algorithm for evaluating the action of the exponential of a sparse matrix on a vector. Rather than reducing the dimension using Krylov subspaces and then applying standard techniques, it uses (2) directly, together with straightforward truncation of the series (1) for (first with , then with etc.). The three key contributions which make this algorithm typically both faster and more accurate than that of Sidje (1998) are detailed in Section 2.1; however, to motivate our own contribution, we briefly note two of these now. Firstly, a detailed error analysis allows for the choice of sensible values for and the truncation point of each series, . The process by which these values are chosen is, however, relatively complex, requiring a look-up table and other bespoke code from Al-Mohy and Higham (2009), and may relate to the reason why the only implementation we were able to find was the original matlab code from the authors. Secondly, the identity

 eMv=eaeM−aIdv (3)

for any scalar, , is used to reduce the Frobenius norm of by setting its trace to .

Our interest lies specifically in the left action of the exponential of a scalar multiple of the rate matrix (or infinitesimal generator), , of a Markov chain, on a probability vector, :

 ν(t)⊤=ν⊤eQt=∞∑i=0tii!ν⊤Qi. (4)

Our key observation is that, since only the diagonal elements of are non-positive, it is possible to shift using (3) so that no elements are negative. Since probability vectors are also non-negative, all elements of all vectors in the series in (4) then become non-negative. Hence there is no cancellation of large positive and negative terms and no need to scale by a factor . This Single Positive Series algorithm (SPS) reduces the number of matrix-vector operations compared to that in the algorithm of Al-Mohy and Higham (2011), typically by a factor of , concomittantly reducing the opportunities for rounding errors and speeding up the evaluation. The SPS algorithm is straightforward to implement, since

is simply a high quantile of the Poisson distribution; we do, however, provide a method for evaluating this quantile which is typically faster than standard methods. Our second contribution lies in the evaluation of (

4) for multiple time points, , for example to chart the evolution of the Markov chain. These distributions are typically evaluated sequentially, since . Using the fact that, in our case, , we provide the Multiple-Use Single Positive Series algorithm (MUSPS) which evaluates all at once from a single evaluation of the series, rather than sequentially.

With our algorithm, or that of Al-Mohy and Higham (2011), for a Markov chain with thousands of states, evaluation of for values of that see a notable change in

takes only a small fraction of a second, enabling direct maximum likelihood estimation of the rate parameters when the network is observed precisely at a discrete set of time points, and making Bayesian inference (via latent variables for the states at observation times) feasible for partially observed networks. As well as being faster than the algorithm of

Al-Mohy and Higham (2011), the SPS algorithm is short, self-contained and straightforward to code should our C++ implementation not fit with the user’s preferred platform. For common epidemic models in particular, the alternative statespace formulation of Ho et al. (2018) (see Section 2.2) makes Bayesian inference feasible for relatively large populations; for example, Ho et al. (2018) perform inference for the SIR model using data from the Ebola outbreak in regions of Guinea. This reduced statespace is usable, and indeed improvable, by the SPS algorithm. In systems biology and elsewhere, reaction networks often have a countably infinite statespace. Georgoulas et al. (2017) provides a method for dealing with such cases using matrix exponentiation by applying the random truncation algorithm of Glynn and Rhee (2014) to the statespace itself; this method would also benefit from using the Single Positive Series algorithm, rather than the algorithm of Al-Mohy and Higham (2011).

In the remainder of this article, Section 1.1 provides several motivating examples, which will form the basis of the numerical comparisons of SPS and MUSPS against competitors in Section 4. Section 2 provides further background on the algorithm of Al-Mohy and Higham (2011) and aspects of the algorithm of Ho et al. (2018), whilst the SPS and MUSPS algorithms themselves are detailed in Section 3.

### 1.1 Examples and motivation

Both by way of motivation and because we shall use them later to illustrate our method, we now present three examples of continuous-time Markov processes, and one example motivated by networks, where a finite, sparse rate matrix contains all of the information about the dynamics.

For each Markov process, the set of possible states can be placed in one-to-one correspondance with a subset of the non-negative integers . The off-diagonal elements of the rate matrix, , are all non-negative, and the th diagonal element is . A chain that is currently in state leaves this state upon the first event of a Poisson process with a rate of ; the state to which it transitions is with a probability of . Whilst the rate matrix, , is a natural description of the process, the likelihood for typical observation regimes involves the transition matrix, , the th element of which is exactly .

Given a network, or graph with vertices (or nodes), the adjacency matrix, is the matrix with and, for , if nodes and are connected then is some positive weight representing the strength of the connection, else . The graph Laplacian is , where is the diagonal matrix with . The negative Laplacian, , has the same properties as a rate matrix and, as we shall see in Example 4, this connection can be related to certain key properties of a graph.

###### Example 1.

The SIR model for epidemics. The SIR model for a disease epidemic has ‘species’: those who are susceptible to the epidemic, , those both infected and infectious, , and those who have recovered from the epidemic and play no further part in the dynamics, . The non-negative counts of each species are denoted by , , and . For relatively short epidemics the population, , is assumed to be fixed, and so the state of the Markov chain, represented by , is subject to the additional constraint of . The two possible reactions and their associated rates are:

 S+IβSI⟶2I,   and   IγI⟶R.
###### Example 2.

The SEIRS model for epidemics. The SEIRS epidemic model starts from the SIR model and inserts an additional state of infected but not yet infections, (‘exposed’). The number of these individuals is denoted by and the state is subject to the constraint that . Removed subjects also gradually lose their immunity so that the reactions are:

 S+IβSI⟶E+I,   EδE⟶I,   IγI⟶R   and   RηI⟶S.
###### Example 3.

The Moran model for allele frequency descibes the time evolution of the frequency of two alleles, and in a population with a fixed size of . Individuals with allele reproduce at a rate of , and those with reproduce at a rate of . When an individual dies it is replaced by the offspring of a parent chosen uniformly at random from the whole population (including the individual that dies). The allele that the parent passes to the offspring usually matches its own, however as it is passed down an allele may mutate; allele switching to with a probability of , and switching to with a probability of . Let be the number of individuals with allele . The two reactions are

 A1λN⟶A2   and   A2μN⟶A1.

Setting , the corresponding infinitesimal rates are

 λN=(1−fN)[αfN(1−u)+β(1−fN)v]   and   μN=fN[β(1−fN)(1−v)+αfNu],

where the unit of time is the expectation of the exponentially distributed time for an individual to die and be replaced. ∎

Our fourth example provides a characterisation of the discrepancy between two graphs in terms of the transmission of information. Given two graphs, and

with the same vertices but different edges (e.g. representing different observations of the connections, possible imputations given incomplete information, or observations of two different types of connections), and Laplacians

and , Hammond et al. (2013) defines the graph diffusion distance between and as:

 ξ(G1,G2;t):=||exp(−L1t)−exp(−L2t)||fr,

where is the Frobenius norm. The idea is that two graphs might be deemed similar if they enable similar patterns of information transmission. The negative Laplacian is a rate matrix, and the th row of represents the distribution over the network after time of a unit of information that was placed on node at time . The Frobenius norm, therefore, represents an average over all possible initial starting nodes of the discrepancies in the time- distributions.

Calculating and storing the full matrix exponential for a set of large graphs is prohibitive both in computational time and storage space; moreover, some initial conditions will be of much more interest than others.

###### Example 4.

Random graphs Given a particular, shared initial distribution across nodes,

, the evolution of the discrepancy between the probability distributions can provide a useful measure of the difference in the transmission properties of the graphs. We have:

 dν(G1,G2;t):=||νTexp(−L1t)−νTexp(−L2t)||1. (5)

Other examples of finite reaction networks include the dimerisation reactions and the Michaelis-Menten reaction kinetics (e.g. Wilkinson, 2012) and other epidemic models such as the SEIR and SIRS models (e.g. Andersson and Britton, 2000).

## 2 Further Background

Redefining , we require , where and are probability -vectors and is a rate matrix for a reaction network. We now detail the algorithm of Al-Mohy and Higham (2011) and summarise the key numerical results from that article. We then briefly describe the alternative statespace formulation for epidemic models, from Ho et al. (2018).

### 2.1 The algorithm of Al-Mohy and Higham

As stated in the introduction the algorithm presented in Al-Mohy and Higham (2011), which we, henceforth refer to as AMH, uses simple series truncation to evaluate using the decomposition in (2). AMH benefits from three key contributions. Firstly, an error analysis allows for sensible choices of and the value, , at which the Taylor series should be truncated; more detail on this is given below. Secondly, the identity (3) is used to reduce the Frobenius norm of by setting its trace to . Finally, the evaluation of each series for is truncated early, at the th term, if the sum of the sizes of the two most recent terms, (where ), is less than the product of some tolerance and the size of the partial sum, . This leads to a problem-specific speed up which can be substantial.

When is required for multiple time points, , the natural procedure is to calculate and then etc. However, when the set of time points is dense this can lead to a large accumulation of rounding errors. AMH reduces the build up of rounding errors by using, for example, whenever there would be no additional computational cost over using .

The complexity of the AMH algorithm, and perhaps the reason why we could only find it in the authors’ original matlab code, lies in the calculation of appropriate values for and , which, in particular, requires further, bespoke matlab code (by the same authors) for approximating the induced norm of a matrix (and its powers) using one or several matrix-vector multiplications for each power. With the hyper parameters for the algorithm given by the authors and hard-coded into the matlab code, except when the approximation of is found to be sufficiently small with ‘sufficient’ depending on the error tolerance chosen (double or single), is required for each between and , and each of these is approximated using a single matrix-vector multiplication.

In Al-Mohy and Higham (2011) the new algorithms are compared against many competitors across a very wide variety of examples. Whilst they do not dominate all other algorithms across all examples, against a given competitor the relevant new algorithm is both more accurate and faster in the vast majority of cases. In particular, the single-time algorithm is found to be both more accurate and typically between and times faster than expv. On the only occasion when the new algorithm is found to slower than expv the relative slow down is due to the matrix norm estimation, which is not required by SPS.

### 2.2 Birth formulation of the SIR and SEIR processes

Ho et al. (2018) reformulate the SIR and SEIR models in terms of births of the (for SEIR only), and species. A recursive formula for the Laplace transform of the transition probability to a given new state in terms of transition probabilities for old states then permits estimation of the transition vector from a known initial starting point in operations, where is the dimension of the statespace. When exact observations are available, the size of the statespace can then be reduced dramatically.

Consider, first, an SIR model with observations of and at times of and . The number of births is , and the number of births is . The total size of the statespace of possible events given the observations is therefore . For an SEIR model a similar argument leads to a product of the number of births of , and .

We may use the same statespace formulation as Ho et al. (2018), provided we include an additional coffin state, , with for all . Any births that would leave the statespace (and hence contradict the observation at time ) instead go to .

## 3 New algorithms

Given a vector , a matrix and a precision , the Single Positive Series (SPS) algorithm estimates , whilst the Multiple-Use Single Positive Series (MUSPS) algorithm estimates for a sequence of times . Unlike any previous algorithms, however, both the SPS and MUSPS algorithms operate by truncating a single series none of whose terms can be negative, rather than truncating multiple series where terms may change sign. Given an , our algorithms calculate such that the (guaranteed to be non-negative) true missing probability mass in each of the dimensions is controlled:

 0<1−d∑i=1^μ∗i<ϵ   and   0<1−d∑i=1^ν(tj)∗i<ϵ   (j=1,…,n),

where and are the probability vectors that would be calculated if there were no rounding errors, and the only errors were due to the truncation of the infinite series. Typically we aim for the error to be similar to the machine’s precision. We control the absolute truncation error since with any truncation of the power series it is impossible to obtain general control of the relative error in a given component of , . Consider, for example, a Moran process (Example 3), where the matrix is tridiagonal. Then is also banded, with a band width of . For any given , and , set . The truncated approximation to gives a transition probability of for all states above , yet, in truth there is a non-zero probability of such a transition.

### 3.1 The Single Positive Series algorithm

Let

 ρ:=maxi=1,…,d|Qii|   and   P=(1/ρ)Q+I. (6)

P is a Markov transition matrix. However, . Thus

 μ⊤=ν⊤eQ=e−ρν⊤∞∑i=0ρii!Pi≈e−ρm∑i=0ρii!ν⊤Pi=:^μ∗⊤.

Since no element of or is negative, there is no possibility for cancellation of terms of opposite sign. Now

 d∑i=1^μ∗i=^μ∗⊤1=e−ρm∑i=0ρii!ν⊤Pi1=e−ρm∑i=0ρii!.

So the absolute relative error, or missing probability mass, due to truncation is

 rm(ρ):=e−ρ∞∑i=m+1ρii!,

the tail probability of a random variable. Of direct interest to us is

 mϵ(ρ):=inf{m∈N:rm(ρ)≤ϵ},

the smallest required to achieve an error of at most , or, essentially, the quantile function for a random variable, evaluated at . Theorem 1, which is proved in Appendix A, provides relatively tight bounds on for any and that might be used.

###### Theorem 1.

If , , and if , . More generally: , where

 m+:=ρ−13logϵ{1+(1−18ρlogϵ)1/2}−1. (7)

Furthermore,

 ⌊m−⌋≤mϵ(ρ)≤⌈m++⌉,

where both inequalities require and the latter also requires and , where

 m− :=ρ+{2ρ}1/2{−log(ϵ√2π)−32logA+log(A−1)}1/2, (8) m++ :=ρ+B−logϵ3{1+(1+18ρB−logϵ)1/2}, (9)
 A:=2ρh(m++1ρ),   B:=−12log4πρh(m−ρ),

and .

The bound (7) arises from a standard argument, whereas those in (8) and (9) are derived from extremely sharp but intractable bounds on in Short (2013); our bounds use only elementary functions and so are much quicker to compute than the quantile upper bound in Short (2013), yet from Figure 1 they are still very sharp. The bounds in (8) and (9) together with the alternative form in (10) permit a simple but fast binary search for (see also Section 3.4). Equation (10) follows from the equivalence between at least events of a Poisson process with a unit rate occuring by time and the time until the th event being at most :

 rm(ρ)=1Γ(m+1)∫ρ0xme−xdx. (10)

As well as providing a tight strict upper bound on the number of sparse vector-matrix multiplications required by our algorithm, the bounds lead directly to the following asymptotic cost:

###### Corollary 1.

For any fixed error tolerance, , .

The Single Positive Series (SPS) algorithm is presented as Algorithm 1. For large values of , although there is no problem with large positive and negative terms cancelling, it is possible that the partial sum might exceed the largest floating point number storable on the machine. Our algorithm circumvents this problem by occasionally renormalising the vector partial sum when the most recent contribution is large, and compensating for this at the end; see lines 5, 12 and 14. One could, alternatively, renormalise when the partial sum itself becomes large, but the translation of the algorithm to multiple time points (Section 3.3) becomes both messier and computationally slower, with no benefits for the accuracy of the algorithm for any values for which the algorithm might feasibly be used.

### 3.2 Options

We now describe two optional extensions to the SPS algorithm: renormalisation and two-tailed truncation.

Since there is no need to keep track of the logarithmic offset ( in Algorithm 1). Instead the final vector ( in Algorithm 1) is renormalised at the end so that its components sum to . Renormalisation when a rate matrix is being exponentiated is included as an option in expoRkit Sidje (1998), for example. In experiments detailed in Section 4.1 and 4.2 we find that it improves the accuracy of the algorithm by between a factor of and several orders of magnitude at no additional computational cost.

In contrast to the general applicability of renormalisation, two-tailed truncation is unique to the SPS algorithm and reduces the computational cost with no loss of accuracy. When is moderate or large, the total mass of probability from the initial value of and the early values accumulated into (Steps 6 and 10 of Algorithm 1) is negligible (has a relative value smaller than , say) compared with the sum of the later values. In such cases may be initialised to and step 10 omitted for values of beneath some . Proposition 1 (see Appendix B for the proof) shows that if is chosen such that then setting ensures that the missing probability mass is no more than . For large , , so with two-tailed truncation the cumulative cost of Step 10 dwindles compared with the other costs, all of which are repeated times.

###### Proposition 1.

Given , let , and let . Then for , .

### 3.3 The Multiple-Use SPS algorithm

The SPS algorithm evaluates for a single time, . Suppose, now, that transitions for multiple time points, , are required. We use the identity

 ν⊤eQt=e−ρtν⊤eρPt=e−ρt∞∑i=0(tt∗)i×ρiti∗i!ν⊤Pi,

where is as defined in (6), and .

Then, with , and as in Algorithm 1, for each additional time interval, , since , we initialise . Within the iterations, we update and as in Algorithm 1 but for some relevant (rather than ) and, effectively, include the additional step for each : . We, therefore, avoid the potential build up of errors from partitioning the total time interval into sub-intervals. Furthermore, we reduce the cost since the vector-matrix multiplication, which involves multiplications and additions is only performed once for each ; for all times , the only update is to , and this involves a single addition and division. Moreover, although the loop over is repeated times, for any given , is only updated when the contribution to the eventual total would be non-negligible.

Both renormalisation and two-tailed truncation can be applied to the MUSPS algorithm if desired, and the algorithm with both of these options (MUSPS2r) typically provides the best combination of speed and accuracy. The bare bones of MUSPS2r is provided in Algorithm 2; for simplicity of presentation the version does not deal with numerical underflow or overflow, and the necessary additional steps are provided in Appendix C.

### 3.4 Implementation details

C++ code for SPS and MUSPS together with some working examples is (for now) besides the link to this paper at https://www.maths.lancs.ac.uk/sherlocc/Publications/publications.html

For a single time point, our binary search algorithm homes in on the required using the upper and lower bounds of Theorem 1 together with the identity (10), the right hand side of which can be evaluated quickly and accurately using the standard C++ toolbox, boost. This is quicker than the standard implementation of the Poisson quantile function (e.g. as implemented in boost), which uses the Cornish-Fisher expansion to approximate the quantile, hence needing an expensive evaluation of , and then conducts a local search. When evaluating for multiple times , presented in ascending order, our function can be speeded up still further for in by using as a lower bound for .

The speed of a vector multiplication by a sparse-matrix depends on the implementation of the sparse matrix algorithm. In R R Core Team (2018) and in C++ Armadillo Sanderson and Curtin (2016, 2018), sparse matrices are stored in column-major order. Hence pre-multiplication of the sparse matrix by a vector, , is much quicker than post multiplcation, . In other languages, such as Matlab, sparse matrices are stored in row-major order and post-multiplication is the quicker operation, so should be stored and used, rather than .

## 4 Numerical comparisons

In Al-Mohy and Higham (2011) the new algorithm (henceforth referred to as AMH) is compared across many examples against state-of-the-art competitors, including, in particular, the expokit function expv Sidje (1998). In most of the experiments, for single times and multiple time points, AMH is found to give comparable or superior accuracy together with superior computational speed. Given these existing comparisons and that comparing direct timings of algorithms coded in different languages confounds the efficiencies of the algorithms and of the compilers, we perform a relatively short comparison of accuracy and clock time between SPS, expv and a bespoke algorithm for epidemic processes, before comparing SPS and AMH in terms of accuracy and computational complexity across a range of examples. We then move on to multiple time points and exhibit additional gains in both speed in accuracy, via the MUSPS algororithm.

The highest accuracy available in C++ using sparse matrices and the armadillo linear algebra library is double precision, which we used throughout in our implementation of both our algorithm and that of Al-Mohy and Higham (2011). For the SPS algorithm we consider two settings: (for accuracy comparable with the algorithm of Ho et al. (2018); see Section 2.2) and , where the missing probability mass is just below the available precision. For the SPS algorithm, following the four options in Section 3.2; the label ‘r’ indicates that renormalisation was used and the label ‘2’ indicates that two-tailed truncation was used.

### 4.1 Comparison with state-of-the art algorithms other than AMH

We consider the collection (see the first three rows of Table 1) of (susceptible and infected) values for the Eyam plague that originated in Raggett (1982) and were used in Ho et al. (2018) to examine the speed and accuracy of the algorithm for the birth representation of the SIR model briefly described in Section 2.2. We set the parameters to their maximum-likelihood estimates, and consider the likelihood for the data in Table 1. In addition, to mimic the size of potential changes between observation times and the size of the elements of the rate matrix from a larger population, we also evaluated the likelihood for the jump directly from the data at time to the data at time . The final two rows of Table 1 refer to the rate matrix for the transition between consecutive observations and provide the dimension the matrix and the absolute value of is largest entry ().

For the SPS algorithm, with and , the algorithm of Ho et al. (2018) and the FORTRAN routine expv from the expoRkit package Sidje (1998) we found the CPU time for estimations of the likelihood ( estimates for the likelihood for the jump from to ); each experiment was performed three times, and the range of timings is given. We also recorded the error in the evaluation of the log likelihood. Given that the true likelihood is not known, the error using SPS with was approximately bounded by examining the discrepancy from the result produced by AMH with chosen for double-precision arithmetic.

The results are presented in Table 2. Since using renormalisation and two-tailed truncation together produced the fastest and most accurate evaluations, we only considered this combination when using . When using expv we did not include the (considerable) cost of creating the sparse rate matrix for each transition, since this was necessarily performed in R. The timings for the SPS algorithm do, however, include matrix-creation costs, which accounted for around of the total cost.

As will be discussed further in Section 4.2, the choice of tolerance, , typically has only a small effect on the speed of SPS. For the full likelihood evaluation, SPS() is over times as fast as HCS and more accurate. The more accurate SPS() is roughly times faster than HCS as well as being more accurate even than expv and over twice as fast as expv; however this timing does not take into account that expv was using eight times the processor power of the other algorithms (expv is automatically optimised to the machine on which it is compiled, and attempting to restrict it to a single core caused it to slow down by a factor of , rather than ). This factor also does not account for the fact that FORTRAN is known to deal more efficiently with matrices than does C or C++ (e.g. https://modelingguru.nasa.gov/docs/DOC-2625).

For the single large jump between observations, we see the same pattern in terms of accuracy. The gain in efficiency by using two-tailed-truncation is larger because is larger ( and ), but despite this, HCS is now more efficient than SPS, although considerably less accurate than SPS(). However, in addition to the savings in state-space size made by Ho et al. (2018), further savings, are possible for SPS and expv. For example, for the SIR model, the constraint that: implies that . For the single big observation jump this gives that , almost halfing the number of states that need to be considered. With this saving, SPS would be of a similar speed to HCS. Analogous savings of up to a factor of would be possible for an SEIR model.

We consider a further example, Example 5, chosen for tractability so that the values can be calculated using greater precision and the true errors ascertained.

###### Example 5.

Immigration-death A single species has members in a confined space with available slots. The species dies out with a rate proportional to the current frequency, and immigration occurs with a rate proportional to the number of empty slots.

 XμX⟶∅   and   ∅γ(n−X)⟶X.

Transition probabilities from are derived in Appendix D. ∎

Setting , , and , we timed calculations of the probability vector and examined . Using 8 cores at , expv took seconds, with an error discrepancy of , whereas on a single core, SPS2r() took seconds with a discrepancy of ; the accuracy of SPSr was almost identical, whereas the discrepancies of both SPS and SPS2 where ; the timings were all within of that of SPS2r, which was the quickest. We then repeated the experiment but with and repetitions. Using cores at , expv took seconds, with an accuracy of , whereas using a single core, SPS2r() took seconds with an accuracy of (SPS2r/SPSr) and (SPS/SPS2).

Even allowing for a reduction in efficiency because of communication costs between cores, the results for the immigration-death model and those from Table 2 suggest that SPS is at least an order of magnitude faster than its nearest generic and easily available competitor, expv from expoRkit, as well as several orders of magnitude more accurate.

### 4.2 Comparison with the algorithm of Al-Mohy and Higham

We now compare the computational complexity and the accuracy of our method and that of Al-Mohy and Higham (2011) across four reaction networks. Complexity is measured in terms of , the number of (sparse) vector-matrix multiplications required in evaluating . The first example is chosen specifically because its transition probabilities are tractable and exact errors in the transition probability vectors can, therefore, be provided. Since an absolute truth is not available for the second and third examples, for a transition of time from the initial condition, accuracy for these examples is indicated by the absolute discrepancy in the probability vector: , where and represent two different estimates of the vector . We refer to the double and single precision error tolerance settings for the algorithm of Al-Mohy and Higham (2011) as AMHd and AMHs respectively (double and single precision arithmetic correspond to relative errors in the representation of any real number of approximately and ). AMHd and AMHs differ only in the choice of and , with AMHs being set up for single precision arithmetic and requiring fewer multiplications in all.

We compare the algorithms on four reaction networks, evaluating transitions of time from an initial condition to a set of logarithmically-spaced values up to a maximum of , which, in all examples, is a natural, large but reasonable value to consider. The first two examples lead to a rate matrix with a maximum (at ) and ; the third leads to a rate matrix with a maximum and the fourth to a matrix with a maximum .

Our first example is the tractable immigration-death model (Example 5), where we used the same parameters as in Section 4.1, with . Our second example is the Moran model (Example 3) with , and . The parameter values imply that the state is absorbing and, more importantly, that the rate matrix is poorly conditioned and exponentiation is prone to producing large errors (see Crawford and Suchard, 2012). Finally, ; i.e., by the final time most of the population has the new allele. The third example is an SIR model (Example 1) for a population of with . Setting , so that at the start of the epidemic an infective is expected to infect approximately four susceptibles, we find that ; there is a high probability that the epidemic is over by the final time considered. The final example is an SEIRS model with , and so that , and epidemics that have not died out have typically suffered from two ‘waves’ of the disease.

In Figure 2 we first examine the immigration-death model, where accuracies are, effectively, exact. For values of , AMHs is much less accurate than AMHd, SPS2 and SPS2r, which all have similar accuracy. However, for larger values of , AMHs, AMHd and SPS2 all have comparable accuracy, which degrades approximately linearly (on the log scale) as , and hence the number of vector-matrix multiplcations and, thus, the accumulated error, increases. The renormalisation step in SPS2r appears to mitigate against much of the degradation with . Accuracy results for SPS with single-tailed truncation are indistinguishable from those for SPS with two-tailed truncation and so are not shown. For the four combinations of SPS with (not shown), for accuracies are between and , with greater accuracy for . SPSr is the least accurate and SPS2 is the most accurate because arithmetic is still carried out to double precision so errors are entirely due to the premature truncation. SPS2 truncates with the upper tail at and the lower strictly below and so enjoys greater accuracy; the renormalisation in SPSr and SPS2r redistributes the missing mass over the existing entries rather than placing some of it on the entries which are currently zero yet should not be.

For small values of , SPS() and AMHd have similar complexities, as do AMHs and SPS(), this latter explaining the relative innacuracy of AMHs compared with AMHd for low values. For moderate and large the complexities of SPS() and SPS() are very similar (since for both) and are roughly half of an order of magnitude smaller than the complexities of both AMHd and AMHs. Indeed, at the largest values the ratio for AMHd to SPS() is . The main reason for the efficiency gain is that to preserve precision in the face of potential cancellation of large positive and negative terms, the algorithm of Al-Mohy and Higham (2011) increases so as to limit the highest series power to . For large this leads to a cost of approximately and for AMHd and AMHs respectively, whereas for our algorithm, from Corollary 1, the cost is approximately . The true cost of AMHs and AMHd, visible in Figure 2, is less than suggested above because, as described in Section 2.1, the algorithms include a stopping rule for early truncation of each series expansion.

Largely the same pattern as above is observed for the Moran, SIR and SEIRS models and so Figure 2 only includes the plots for the SIR model. The following minor differences are evident: (i) for low values, AMHs is (Moran), (SIR) and (SEIRS) rather than orders of magnitude less accurate than AMHd and SPS(); (ii) with no ground truth we are unable to ascertain the relative accuracies of AMHd and SPS() but can observe that their discrepancy increases with rho and that they agree to within approximately for ; (iii) for SPS() renormalisation improves the accuracy slightly for the SIR and SEIRS models but reduces it slightly for the Moran and the Immigration-Death models. Using two-tailed truncation always improves the accuracy slightly. At the largest values the ratio of complexities for AMHd to SPS() are (Moran), (SIR) and (SEIRS).

For the algorithm of Al-Mohy and Higham (2011), because of the further complication of including the norm-approximation algorithm of Al-Mohy and Higham (2009), we did not approximate the norms , but calculated them exactly. The cost of these evaluations was not included in the algorithm complexity; the vector-matrix multiplications that would have been used for these calculations in Al-Mohy and Higham (2009) were also not included in the complexity.

### 4.3 Simultaneous evaluation at multiple time points

Finally, we investigate the MUSPS algorithm of Section 3.3; its accuracy, speed and usefulness. We examine the immigration-death and SEIRS algorithms with the parameterisation as specified in Section 4.2 and four variations on Barábasi-Albert random network model that will be described, briefly, below. For each system we considered time points, equally spaced between and , where for the immigration-death system, for the SEIRS system, and for each of the four the random networks; in each case the upper bound, , was sensible given the system and (for SEIRS and random graphs) application of interest. Using time points for the intractable systems was not feasible as the cost of obtaining a set of appoximate truths using SPS (see below) is quadratic in .

We created a set of four random networks based upon a pair of independent realisations from the Barabasi-Albert model (see Appendix E), and , with parameters . Each of and has nodes, with an average degree of , a minimum degree of and the probability that the degree of a randomly chosen node is being proportional to (e.g. Durrett, 2007, Chapter 4). Our four graphs are obtained by connecting and through a single edge with a weight of , as follows. Let and be the nodes in and with the highest degree, and let and be the last node added to and respectively (this has the lowest possible degree). Joining to gives , joining to gives , joining to gives and joining to gives . We took and .

For all three applications we ran MUSPS and a sequential SPS which evolved a current vector by a time of to obtain a new current vector times. We compared the speeds of the two algorithms and their accuracies across all of the times from to ; for the immigration-death system, accuracy was relative to the known true transition probabilities; for the other two systems it was relative to an individual, direct calculation for each using SPS. We used MUSPS2r and SPS2r throughout.

Each application was run times; Table 3 shows the dimension and value for each system as well as the range of (i) times for MUSPS and (ii) ratios of the time for MUSPS to the time for sequential SPS. For the random networks model, for a given realisation of and , the variation in timings for was less than , so we only show the results for . Total time is roughly proportional to , as expected. For SEIRS and the Barabási-Albert model MUSPS is more efficient than SPS. For the immigration-death system MUSPS is more expensive than sequential SPS for two reasons. Firstly, the dimension is small so the relative cost of the -dependent overheads of the MUSPS algorithm is large; moreover, has only entries per column and so each vector-matrix multiplication has a similar cost to the vector-scalar multiplication and vector-vector addition that MUSPS must still perform. The additional cost may, however, be justified if accuracy is important, as we shall see.

Figure 3 plots the accuracy (base logarithm of the discrepancy from the ‘truth’, with removed as all algorithms produce the same probability vector) of MUSPS against that of SPS for the immigration-death system, the SEIRS system and the first set of realisations from the random-network model (the pattern was persistent across replications). Whenever there is a substantial difference in the accuracies, MUSPS is more accurate than sequential SPS.

The left and central panels of Figure 4 demonstrates an example use of MUSPS within an SEIRS model by plotting the probability of extinction against (left panel) and the expected number of people infected or exposed against amongst infections which have not died out (central panel). The kink in the first graph around corresponds to the peak in the second. However, the central panel shows a dip in the presence of infections before (conditional) equilibrium is reached. Interestingly, the solution to the initial value problem with and

 dSdt=−θ1SI+θ4(npop−S−E−I), dEdt=θ1SI−θ2E, dIdt=θ2E


text:
# Synchrocyclotron

Sketch of a synchrocyclotron from McMillan's patent.[1]

A synchrocyclotron is a special type of cyclotron, patented by Edwin McMillan, in which the frequency of the driving RF electric field is varied to compensate for relativistic effects as the particles' velocity begins to approach the speed of light. This is in contrast to the classical cyclotron, where this frequency is constant.[1]

There are two major differences between the synchrocyclotron and the classical cyclotron. In the synchrocyclotron, only one dee retains its classical shape, while the other pole is open (see patent sketch). Furthermore, the frequency of oscillating electric field in a synchrocyclotron is decreasing continuously instead of kept constant so as to maintain cyclotron resonance for relativistic velocities. One terminal of the oscillating electric potential varying periodically is applied to the dee and the other terminal is on ground potential. The proton or deuterons to be accelerated are made to move in circles of increasing radii. The acceleration of particles takes place as they enter or leave D. At the outer edge, the ion beam can be removed with the aid of electrostatic deflector. It was possible to produce 200 MeV deuterons and 400 MeV α-particle with the first synchrocyclotron.[citation needed]

## Differences to the classical Cyclotron

In a classical cyclotron, the angular frequency of the electric field is given by

${\displaystyle \omega ={\frac {qB}{m}}}$,

Where ${\displaystyle \omega }$ is the angular frequency of the electric field, ${\displaystyle q}$ is the charge on the particle, ${\displaystyle B}$ is the magnetic field, and ${\displaystyle m}$ is the mass of the particle. This makes the assumption that the particle is classical, and doesn't experience relativistic phenomena such as length contraction. These effects start to become significant when ${\displaystyle v}$, the velocity of the particle greater than ${\displaystyle \approx {\frac {c}{3}}}$. To correct for this, the relativistic mass is used instead of the rest mass, to do this, a factor of ${\displaystyle \gamma }$ multiplies the mass. Such

${\displaystyle \omega ={\frac {qB}{m\gamma }}}$,

Where,

${\displaystyle \gamma ={\frac {1}{\sqrt {1-{\frac {v^{2}}{c^{2}}}}}}}$,

This is then the applied frequency angular frequency to the particles as they are accelerated around the synchrocyclotron.

A part of the former Orsay synchrocyclotron

The chief advantage of the synchrocyclotron is that there is no need to restrict the number of revolutions executed by the ion before its exit. As such, the potential difference supplied between the dees can be much smaller.

The smaller potential difference needed across the gap has the following uses:

1. There is no need for a narrow gap between the dees as in the case of conventional cyclotron, because strong electric fields for producing large acceleration are not required. Thus only one dee can be used instead of two, the other end of the oscillating voltage supply being connected to earth.
2. The magnetic pole pieces can be brought closer, thus making it possible to increase greatly the magnetic flux density.
3. The frequency valve oscillator is able to function with much greater efficiency.

format: plain


text:
# Attempt to export EPS figures from TikZ fails

I'm attempting to output EPS figures when using MikTeX 2.8 under Windows XP, and am following Export eps figures from TikZ. My problem is that I can't get the bare example given there to run.

The code I'm using was copied from that thread verbatim:

\documentclass{article}
\usepackage{tikz}

% set up externalization
\usetikzlibrary{external}
\tikzset{external/system call={latex \tikzexternalcheckshellescape -halt-on-error
-interaction=batchmode -jobname "\image" "\texsource" &&
dvips -o "\image".ps "\image".dvi}}
\tikzexternalize[shell escape=-enable-write18] % MikTeX uses a -enable-write18 instead of --shell-escape.

\begin{document}
\begin{tikzpicture}
\draw (0,0) circle (1cm);
\end{tikzpicture}
\end{document}

When running I get the following error message:

This is pdfTeX, Version 3.1415926-1.40.10 (MiKTeX 2.8)
entering extended mode
(C:/junk/untitled-2.tex
LaTeX2e <2009/09/24>
. . .
)
===== 'mode=convert with system call': Invoking 'latex -halt-on-error -interaction=batchmode -jobname "untitled-2-figure0" "\def\tikzexternalrealjob{untitled-2}\input{untitled-2}" && dvips -o "untitled-2-figure0".ps "untitled-2-figure0".dvi' ========

! Package tikz Error: Sorry, the system call 'latex -halt-on-error -interaction=batchmode -jobname "untitled-2-figure0" "\def\tikzexternalrealjob{untitled-2}\input{untitled-2}" && dvips -o "untitled-2-figure0".ps "untitled-2-figure0".dvi' did NOT result in a usable output file 'untitled-2-figure0' (expected one of .pdf:.jpg:.jpeg:.png:). Please verify that you have enabled system calls. For pdflatex, this is 'pdflatex -shell-escape'. Sometimes it is also named 'write 18' or something like that. Or maybe the command simply failed? Error messages can be found in 'untitled-2-figure0.log'. If you continue now, I'll try to typeset the picture.
See the tikz package documentation for explanation.
Type  H <return>  for immediate help.
...

l.14 \end{tikzpicture}

?

Can someone help me over this hurdle?

My previous post was written on my Windows XP system using Chrome, where the answer box limited me to about 600 characters and did not accept markup or newlines :-(. I'm writing this on my Windows 7 system with Chrome, and the answer box has all markup bells and whistles I need -- very strange.

Following is the compilation output:

This is pdfTeX, Version 3.1415926-1.40.10 (MiKTeX 2.8) entering extended mode ("C:/Documents and Settings/Moshe/My Documents/test.tex" LaTeX2e <2009/09/24> . . .

===== 'mode=convert with system call': Invoking 'latex -enable-write18 -halt-on
-error -interaction=batchmode -jobname "test-figure0" "\def\tikzexternalrealjob {test}\input{test}" && dvips -o "test-figure0".ps "test-figure0".dvi' ========

! Package tikz Error: Sorry, the system call 'latex -enable-write18 -halt-on-er ror -interaction=batchmode -jobname "test-figure0" "\def\tikzexternalrealjob{te st}\input{test}" && dvips -o "test-figure0".ps "test-figure0".dvi' did NOT resu lt in a usable output file 'test-figure0' (expected one of .pdf:.jpg:.jpeg:.png :). Please verify that you have enabled system calls. For pdflatex, this is 'pd flatex -shell-escape'. Sometimes it is also named 'write 18' or something like  that. Or maybe the command simply failed? Error messages can be found in 'test- figure0.log'. If you continue now, I'll try to typeset the picture.

See the tikz package documentation for explanation. Type  H <return>  for immediate help.  ...

l.14 \end{tikzpicture}

I hope this is easier to read than my previous post :-)

-
You are compiling it with the -enable-write18 command line option, do you? It is required for this to work. – Martin Scharrer Mar 22 '11 at 8:38
Thank you for your prompt reply! Do I need to do anything other than include the line "\tikzexternalize[shell escape=-enable-write18]" in the script? I'm running the GUI version of MikTeX. Should I be setting some configuration option for -enable-write18? I did enable the "Allow scripts to running system commands" configuration option, but received the same error message. – Moshe Rubin Mar 22 '11 at 8:46
I'm not using MikTeX by myself much, so I can't tell you the exact way on Windows. Anyway, the -enable-write18 option must be passed to latex. The \tikzexternalize setting is AFAIK only for sub-calls of latex. Try to move this line before \tikzset and have a lock at the untitled-2-figure0.log file mentioned in the error message. It might contain more information. – Martin Scharrer Mar 22 '11 at 9:03
Thanks for pointing me to the log file (duh!). It seems to indicate a quotation error: runsystem(latex -halt-on-error -interaction=batchmode -jobname "untitled-2-figure0" "\def\tikzexternalrealjob{untitled-2}\input{untitled-2}" && dvips -o "untitled-2-figure0".ps "untitled-2-figure0".dvi)...quotation error in system command. – Moshe Rubin Mar 22 '11 at 9:26
My first instinct is also that the main latex process is not being allowed to run system commands. From the MiKTeX manual I read "The \write18 feature is only partially enabled by default to avoid security problems: only well known commands are allowed. You fully enable the feature by specifying --enable-write18 on the TeX command-line." So you could try doing a commandline compilation with pdflatex --enable-write18 and see if that works. – Loop Space Mar 22 '11 at 9:29

I tried to reproduce your example on my MikTeX 2.9 installation and found the same error message. Interestingly, the main .log file contains the hint "quotation error in system call". When I dropped all quotation marks, I got "forbidden, restricted" or something like that. Seems as if MikTeX has strange error checking here.

Anyway, as soon as I reconfigured TeXnicCenter to add the -enable-write18 switch, it worked - and it no longer complaint about 'quotation errors'.

So, I agree with the suggestions shown in comments above: I am pretty sure that you "only" need to figure out, how to add the -enable-write18 switch to your TeX Editor. If you tell us which editor you use, we may be able to assist you.

For TeXnicCenter, it works as follows: go to "Output >> Define output profiles (Alt+F7)". There, you should see the different output profiles on the left and on the right, you see the involved programs like "Path to LaTeX Compiler", "BibTeX", "MakeIndex". Select the profile LaTeX=>PS, then add -enable-write18 to the "Path to LaTeX Compiler". The "Path To LaTeX Compiler" should be something like C:Program Files\MikTeX 2.9\miktex\bin\latex.exe -enable-write18 afterwards.

I am using WinEdt, and having the same error. How ca i -enable-write18 in winEdt???. – Misaki Jul 17 '14 at 12:57


text:
# Schwartz space

For the Schwartz space of a semisimple Lie group, see Harish-Chandra's Schwartz space.

In mathematics, Schwartz space is the function space of functions all of whose derivatives are rapidly decreasing (defined rigorously below). This space has the important property that the Fourier transform is an automorphism on this space. This property enables one, by duality, to define the Fourier transform for elements in the dual space of S, that is, for tempered distributions. The Schwartz space was named in honour of Laurent Schwartz by Alexander Grothendieck.[1] A function in the Schwartz space is sometimes called a Schwartz function.

A two-dimensional Gaussian function is an example of a rapidly decreasing function.

## Definition

The Schwartz space or space of rapidly decreasing functions on Rn is the function space

${\displaystyle S\left(\mathbf {R} ^{n}\right)=\left\{f\in C^{\infty }(\mathbf {R} ^{n}):\|f\|_{\alpha ,\beta }<\infty \quad \forall \alpha ,\beta \in \mathbb {Z} _{+}^{n}\right\},}$

where α, β are multi-indices, C(Rn) is the set of smooth functions from Rn to C, and

${\displaystyle \|f\|_{\alpha ,\beta }=\sup _{x\in \mathbf {R} ^{n}}\left|x^{\alpha }D^{\beta }f(x)\right|.}$

Here, sup denotes the supremum, and we again use multi-index notation.

To put common language to this definition, we could note that a rapidly decreasing function is essentially a function f(x) such that f(x), f′(x), f′′(x), ... all exist everywhere on R and go to zero as x → ±∞ faster than any inverse power of x. In particular, S(Rn) is a subspace of the function space C(Rn) of infinitely differentiable functions.

## Examples of functions in the Schwartz space

• If i is a multi-index, and a is a positive real number, then
${\displaystyle x^{i}e^{-a|x|^{2}}\in S(\mathbf {R} ^{n}).}$
• Any smooth function f with compact support is in S(Rn). This is clear since any derivative of f is continuous and supported in the support of f, so (xαDβ) f has a maximum in Rn by the extreme value theorem.

## Properties

• S(Rn) is a Fréchet space over the complex numbers.
• From Leibniz' rule, it follows that S(Rn) is also closed under pointwise multiplication: if f, gS(Rn), then fgS(Rn).
• If 1 ≤ p ≤ ∞, then S(Rn) ⊂ Lp(Rn).
• The space of all bump functions, C
c
(Rn), is included in S(Rn).
• The Fourier transform is a linear isomorphism S(Rn) → S(Rn).
• If fS(R), then f is uniformly continuous on R.

## References

1. ^ TerzioĞglu, T. (1969). On Schwartz spaces. Mathematische Annalen, 182(3), 236–242.
• Hörmander, L. (1990). The Analysis of Linear Partial Differential Operators I, (Distribution theory and Fourier Analysis) (2nd ed.). Berlin: Springer-Verlag. ISBN 3-540-52343-X.
• Reed, M.; Simon, B. (1980). Methods of Modern Mathematical Physics: Functional Analysis I (Revised and enlarged ed.). San Diego: Academic Press. ISBN 0-12-585050-6.
• Stein, Elias M.; Shakarchi, Rami (2003). Fourier Analysis: An Introduction (Princeton Lectures in Analysis I<). Princeton: Princeton University Press. ISBN 0-691-11384-X.


text:
# Exercises

## Submission

To submit your solutions create a folder called lab3 and push one or more file to your repository with your answers (it can be plain text, markdown, pdf or whatever other format is reasonably easy to read)

### Individual

Please push the deliverables into your personal repository, for math-related questions LaTeX is preferred but handwritten is accepted too.

### Team

Create a new repository in VNAV2021-submission called TEAM_<N>, replacing <N> with your team number. For example team 2 will create the repository TEAM_2.

Please push the source code of the entire controller_pkg package in the foder lab3 of the team repository.

Deadline: the VNAV staff will clone your repository on September 29th, 2021 at midnight.

## Individual

### Deliverable 1 - Transformations in Practice (10 pts)

This exercise will help you to understand how to perform transformation between different rotation representations using ROS. In addition, you will gain experience of reading through a documentation in order to find what you need.

#### A. Message vs. tf

In ROS, there are two commonly used datatypes for quaternions. One is geometry_msgs::Quaternion and the other is tf2::Quaternion. Based on the documentation for tf2 and its package listing page, answer the following questions in a text file called deliverable_1.txt:

1. Assume we have an incoming geometry_msgs::Quaternion quat_msg that holds the pose of our robot. We need to save it in an already defined tf2::Quaternion quat_tf for further calculations. Write one line of C++ code to accomplish this task.

2. Assume we have just estimated our robot’s newest rotation and it’s saved in a variable called quat_tf of type tf2::Quaternion. Write one line of C++ code to convert it to a geometry_msgs::Quaternion type. Use quat_msg as the name of the new variable.

3. If you just want to know the scalar value of a tf2::Quaternion, what member function will you use?

#### B. Conversion

The following questions will prepare you for converting between different rotation representations in C++ and ROS.

1. Assume you have a tf2::Quaternion quat_t. How to extract the yaw component of the rotation with just one function call?

2. Assume you have a geometry_msgs::Quaternion quat_msg. How to you convert it to an Eigen 3-by-3 matrix? Refer to this and this for possible functions. You probably need two function calls for this.

### Deliverable 2 - Modelling and control of UAVs (20 pts)

The figure above depicts two quadrotors (a) and (b). Quadrotor (a) is a fully functional UAV, while for Quadrotor (b) someone changed propellers 3 and 4 and reversed their respective rotation directions.

Show mathematically that quadrotor (b) is not able to track a trajectory defined in position $[x,y,z]$ and yaw orientation $\Psi$.

• Hint: write down the $F$ matrix (see lecture notes eq. (6.9)) for the two cases (a) and (b) and compare the rank of the two matrices

Assume that empirical data suggest you can approximate the drag force (in the body frame) of a quadrotor body as:

$F^b = \begin{bmatrix} 0.1 & 0 & 0 \\ 0 & 0.1 & 0 \\ 0 & 0 & 0.2\end{bmatrix} (v^b)^2$

With $(v^b)^2= [-v^b_x \mid v^b_x \mid ,\, -v^b_y \mid v^b_y \mid ,\, -v^b_z \mid v^b_z \mid]^T$, and $v_x$, $v_y$, $v_z$ being the quadrotor velocities along the axes of the body frame.

With the controller discussed in class (see referenced paper 1), describe how you could use the information above to improve the tracking performance.

• Hint: the drag force introduced above is an additional term in the system’s dynamics, which the controller could compensate for explicitly…

## Team

### Trajectory tracking for UAVs

In this section, we are going to implement the geometric controller discussed in class on a simulated drone.

#### Getting the codebase

First, let us install some prerequisites. In a terminal, run:

sudo apt install ros-melodic-mav-msgs ros-melodic-ackermann-msgs

Consistently with the folder structure introduced in previous assignments, we will navigate to our existing clone of the VNAV Labs repository and pull the latest updates:

# change the folder name according to your setup
cd ~/labs
git pull

In ~/labs/lab3 we now have the tesse_ros_bridge and the controller_pkg folders, which are both ROS packages. Let us copy these folders in our VNAV workspace and build the workspace as follows:

cp -r ~/labs/lab3/. ~/vnav_ws/src
cd ~/vnav_ws

In your src folder, you should see three folders: catkin_simple, controller_pkg and tesse-ros-bridge. Then, we need install some Python dependencies:

cd ~/vnav_ws/src/tesse-ros-bridge/tesse-interface
pip install -r requirements.txt  # Dependencies
pip install .

At this point, you can both keep your lab 2 code or remove all the folders related to the previous labs. Either way, you need to compile the new code.

catkin build
source devel/setup.bash

After doing so, please proceed here to download the binary executable for the simulator we are going to use for lab 3. Unzip the file, put the folder in ~/vnav-builds, and run the following commands:

cd ~/vnav-builds/vnav-2021-lab3/
chmod +x vnav-2021-lab3.x86_64

This will make sure you can run the executable in the command line.

### Launching the TESSE simulator with ROS bridge

In a terminal window, run this command to start the simulator:

cd ~/vnav-builds/vnav-2021-lab3/
./vnav-2021-lab3.x86_64

One the simulator starts, you should see the unity simulator (fig. below).

In the simulator, you can:

• press W to turn on all propellers;
• press R to respawn the quadrotor at the initial location.

Note that you can resize the window to help with debugging later on.

Once the workspace has been compiled successfully, with your workspace sourced execute:

# make sure to source devel/setup.bash!

This launches two ROS nodes: tesse_ros_brige and quadrotor_control_interface. The former receives state messages from our simulator and publishes them in a few ROS topics, and the latter listens to our controller node and sends the propeller commands to our simulator. To visualize the state of our quadrotor, we will use rviz. Run the following command to load our rviz preset:

cd ~/vnav_ws/src/controller_pkg
rviz -d rviz/lab3.rviz

You should see a screen similar to this:

A few highlights:

• the purple upward arrow corresponds to the IMU message;
• the base_link_gt axes correspond to the body frame of the quadrotor;
• the world axes correspond to the world frame in ROS.

### Implement the controller

From the previous section, you must have noticed that the drone is stationary. Unsurprisingly, this happens because the UAV is not receiving any control inputs. In this section, you will be precisely asked to take care of it!

Given four propeller speeds, the simulator propagates forward the UAV dynamics and outputs the current state of the robot. It is your job now to provide the propeller speeds based on the current state of the robot (from the simulator) and a desired one (from a trajectory generator). This problem is extensively described in the geometric controller paper 1, where the authors introduce an almost-globally convergent geometric controller for a quadrotor, also outlined in class. We will ask you to read the paper thoroughly and implement this controller in ROS.

For the purposes of this lab, the trajectory generator is given by the traj_publisher node, which outputs a time-varying desired state along a hard-coded circular trajectory. The navigation stack (almost) implemented in our codebase has the architecture shown below:

where the ROS nodes are highlighted as vertices in the graph and the ROS topic names are reported on the edges.

In particular:

• /desired_state is of type trajectory_msgs/MultiDOFJointTrajectoryPoint and has the following fields:
• transforms - type geometry_msgs/Transform[]:
• translation: a 3D vector with the desired position of the c.o.m. in the world frame
• rotation: a quaternion rapresenting a rotation (in our case will be only a yaw component)
• velocities - type geometry_msgs/Twist[]:
• linear: a 3D vector with the desired velocity of c.o.m. in world frame
• angular: we will ignore this field as it is not used by the controller
• accelerations - type geometry_msgs/Twist[]:
• linear: a 3d vector with the desired acceleration of c.o.m. in world frame
• angular: we will ignore this field as it is not used by the controller
• /current_state is of type nav_msgs/Odometry with the relevant fields:
• pose.pose - type geometry_msgs/Pose:
• position: a 3d vector with the current position of the c.o.m. in the world frame
• orientation: a quaternion representing the current orientation of the UAV
• twist.twist - type geometry_msgs/Twist:
• linear: a 3d vector with the current velocity of c.o.m. in the world frame
• angular: a 3d vector with the current angular velocity in the world frame
• /rotor_speed_cmds is of type mav_msgs/Actuators and it contains the relevant field:
• angular_velocities - type float64[], containing the desired speeds of the propellers.

### Simulator conventions

In the picture below, we illustrate the conventions used in unity simulator for the propeller speeds and their positions in the body frame. Please look carefully as this will help you build the correct wrench-to-controls matrix (F2W) in the coding assignment

Note that

• as in most research papers, we assume that, for $\omega_i > 0$, all the propellers produce an upward force (lift). This is consistent with the first row of the matrix in eq. (6.9) on the lecture notes. While this might be confusing at first, note that on real quadrotors the motors are configured to rotate in different directions: two counter-clockwise (CCW), and two clockwise (CW). This configuration is necessary so that our quadrotor will be stable in its yaw direction. There are also two sets of propellers: one set for CCW motors, and another set for CW motors (both will produce upward force if mounted appropriately). See this video for an example.
• in the figure, we show the positive spinning directions for the propellers (blue) and the resulting direction of the yaw torques generated by aerodynamic drag. For example: propeller 1’s angular velocity $\omega_1$ is considered positive when propeller 1 spins clockwise (as seen from top), which in turn produces a counter-clockwise drag torque $\tau_{\text{drag}_1}$. The sign of the drag torque is determined by applying right-hand rule: positive for counter-clockwise, and negative for clockwise.
• the spinning directions are different from the paper. The expression $\tau_{\text{drag}}^i = (-1)^i c_d\,\omega_i | \omega_i |$ implies that the drag on propeller 1 produces a negative (clockwise) yaw torque instead, for $\omega_1$. However, the above convention is acceptable: propellers are mounted differently in the simulated quadrotor!
• the propeller speeds need to be provided in radians per second (rad/s), NOT in rpm. Initially, it will be tempting to go for the latter, as it will make the drone move more easily (control inputs get bigger), however, it will hurt you in the long run.

Before proceeding, keep in mind that there are a number of caveats that will make the math in your code slightly different from the one found in the paper. The main differences can be summarized as follows:

Reference frames. The paper and your code use different conventions for coordinate frames. Please refer to the figure and pay attention to the two main facts summarized below.

1. In the paper, the $z$-axes - both in the world and body frames - are downwards, while ROS has them pointing upwards. This affects the signs of $b_{3_d}$ in equation (12), as well as gravity and aerodynamic forces in equations (12), (15), (16). (equation numbers refer to the paper 1).
2. In the paper, the $x$-body axis is along a propeller arm, while we prefer having it at 45° between two propellers. This will require some thought on your end when converting the total desired force + torques (wrench) into the desired propeller speeds. In other words, you have to change equation (1) from the paper. For this step, refer to equation (6.9) in the lecture notes!

Aerodynamic coefficients. In the paper, equation (1), the quantity $c_{\tau f}$ relates directly the yaw torque with the lift produced by a propeller. In contrast, we prefer using the lift and drag coefficients $c_f$ and $c_d$, which have a clearer physical interpretation. Consider the relation $c_{\tau f}=c_d/c_f$ between these coefficients.

### Deliverable 3 - Geometric controller for the UAV (50 pts)

After reading the reference 1 - thoroughly - get together with your teammates and open the source file controller_node.cpp. In this file, you will find detailed instructions to fill in the missing parts and make your quadrotor fly! To have full credits for this part of the lab, your team needs to complete the following:

• Implement all missing parts in controller_node.cpp and ensure everything compiles.
• Tune parameters so that the quadrotor achieves stable circular trajectory.
• A video showing the quadrotor completing one round of the circular trajectory in rviz (can either be a screen capture, or a video shot using your phone). Please upload the video onto either Google drive or Dropbox, generate a publicly viewable link, and put the link in a text file called rviz_drone_flight.txt in your repo.

#### How to test

catkin build

from the folder ~/vnav_ws.

To try out your code, launch the simulator. Additionally, run your controller by typing, in a separate terminal:

roslaunch controller_pkg traj_tracking.launch

Depending on your focus - and degree of luck - it will take multiple attempts to have a stable and smooth controller and multiple iterations of parameter tuning. That’s normal!

As a suggestion, you could start by testing a stationary waypoint, before moving to the circular trajectory. This is possible by switching the STATIC_POSE  flag to 1 in controller_pkg/src/traj_publisher.cpp. Also, you should use rviz to visually compare the desired and the current poses.

#### Tuning tips

First of all, look closely at the equations and see what are the gain parameters directly affecting. Which one affects the yaw? Which one affects the speed? Which one affects the position? Doing so will give you a rough idea on what to do for tuning.

To provide a more graphical/qualitative understanding, here are some typical situations that you might encounter:

komega is too high:

kr is too low:

kv is too low:

NOTE: The specific gain parameters may differ across different computers.

#### What to expect

If your controller is working reasonably well we expect to see something like this:

Note: it is totally fine if the drone does a rapid maneuver to start tracking the circle as soon as you run the controller node.

Good luck!

## References

1. Lee, Taeyoung, Melvin Leoky, N. Harris McClamroch. “Geometric tracking control of a quadrotor UAV on SE (3).” Decision and Control (CDC), 49th IEEE Conference on. IEEE, 2010 Link  2 3 4


text:
Search packages:
 Sourcecode: c-icap version 1:0.1.4-11:0.1.5-11:0.1.6-11:0.1.6-1ubuntu120080706rc320080706rc3-120080706rc3-2

/*
*  Copyright (C) 2004-2008 Christos Tsantilas
*
*  This program is free software; you can redistribute it and/or
*  modify it under the terms of the GNU Lesser General Public
*  version 2.1 of the License, or (at your option) any later version.
*
*  This program is distributed in the hope that it will be useful,
*  but WITHOUT ANY WARRANTY; without even the implied warranty of
*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
*  Lesser General Public License for more details.
*
*  You should have received a copy of the GNU Lesser General Public
*  License along with this library; if not, write to the Free Software
*  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
*  MA  02110-1301  USA.
*/

#include "c-icap.h"

/**
\ingroup API
*/

ICAP_FROM, ICAP_HOST, ICAP_REFERER,
ICAP_USERAGENT,ICAP_PREVIEW
};

enum ci_encapsulated_entities {ICAP_REQ_HDR, ICAP_RES_HDR,
ICAP_REQ_BODY, ICAP_RES_BODY,
ICAP_NULL_BODY,ICAP_OPT_BODY };
CI_DECLARE_DATA extern const char *ci_encaps_entities[];

#ifdef __CYGWIN__

const char *ci_encaps_entity_string(int e);

#else

#define ci_encaps_entity_string(e) (e<=ICAP_OPT_BODY&&e>=ICAP_REQ_HDR?ci_encaps_entities[e]:"UNKNOWN")

#endif

/**
* This is a struct which can store a set of headers.
* The developers should not touch ci_headers_list_t objects directly but better use the
* documented macros and functions
*/
int size;
int used;
int bufsize;
int bufused;
char *buf;
int packed;

00077 typedef struct ci_encaps_entity{
int start;
int type;
void *entity;
} ci_encaps_entity_t;

#define BUFSIZE          4096

/**
* Allocate memory for a ci_headers_list_t object and initialize it.
\return the allocated object on success, NULL otherwise.
*
*/

/**
\param heads is a pointer to the ci_headers_list_t object to be destroyed
*
*/

/**
* Resets and initialize a ci_headers_list_t object
*
*/

/**
*
*example usage:
\code
\endcode
*
*/

/**
\return non zero on success zero otherwise
*/

/**
\return non zero on success, zero otherwise
*
*example usage:
\code
\endcode
*
*/

/**
\return a pointer to the start of the first occurrence of the header on success, NULL otherwise
*
*example usage:
\code
\endcode
* In this example on success the head pointer will point to a \em "Content-Lenght: 1025" string
*
*/

/**
* Search for a header in a header list and return the value of the first occurrence of this header
\return a pointer to the start of the header on success, NULL otherwise
*
*example usage:
\code
int content_length;
\endcode
*
*/

/*The following headers are only used internally */


text:
# Regulator constants of integral representations of finite groups

```@article{Torzewski2018RegulatorCO,
title={Regulator constants of integral representations of finite groups},
author={Alex Torzewski},
journal={Mathematical Proceedings of the Cambridge Philosophical Society},
year={2018},
volume={168},
pages={75 - 117}
}```
• A. Torzewski
• Published 30 March 2017
• Mathematics
• Mathematical Proceedings of the Cambridge Philosophical Society
Abstract Let G be a finite group and p be a prime. We investigate isomorphism invariants of \$\mathbb{Z}_p\$[G]-lattices whose extension of scalars to \$\mathbb{Q}_p\$ is self-dual, called regulator constants. These were originally introduced by Dokchitser–Dokchitser in the context of elliptic curves. Regulator constants canonically yield a pairing between the space of Brauer relations for G and the subspace of the representation ring for which regulator constants are defined. For all G, we show…
1 Citations

### Functoriality of motivic lifts of the canonical construction

Let \$\$(G,{\mathfrak {X}})\$\$ ( G , X ) be a Shimura datum and K a neat open compact subgroup of \$\$G(\mathbb {A}_f)\$\$ G ( A f ) . Under mild hypothesis on \$\$(G,{\mathfrak {X}})\$\$ ( G , X ) ,

## References

SHOWING 1-10 OF 28 REFERENCES

### Regulator constants and the parity conjecture

• Mathematics
• 2007
AbstractThe p-parity conjecture for twists of elliptic curves relates multiplicities of Artin representations in p∞-Selmer groups to root numbers. In this paper we prove this conjecture for a class

### Regulator constants of integral representations, together with relative motives over Shimura varieties

This thesis is split into three largely independent chapters. The first concerns the representation theory of Zp[G]-lattices. Specifically, we investigate regulator constants, due to

### Representations of Cyclic Groups in Rings of Integers, II

• Mathematics
• 1962
Let Gk be a cyclic group of order k, and let ZGk denote its group ring over the ring Z of rational integers. We denote by n(ZGk) the number of non-isomorphic indecomposable left ZGk-modules having

### On Brauer–Kuroda type relations of S-class numbers in dihedral extensions

Abstract Let F/k be a Galois extension of number fields with dihedral Galois group of order 2q, where q is an odd integer. We express a certain quotient of S-class numbers of intermediate fields,

### Integral Representations of Dihedral Groups of Order 2p

Introduction. Information about the integral representations of finite groups has been obtained to varying extents. For Z the ring of rational integers and G the cyclic group of prime order, the

### Brauer relations in finite groups

• Mathematics
• 2015
If GG is a non-cyclic finite group, non-isomorphic GG-sets X,YX,Y may give rise to isomorphic permutation representations C[X]≅C[Y]C[X]≅C[Y]. Equivalently, the map from the Burnside ring to the

### Brauer relations in finite groups II: Quasi-elementary groups of order paq

• Mathematics
• 2013
Abstract. This is the second in a series of papers investigating the space of Brauer relations of a finite group, the kernel of the natural map from its Burnside ring to the rational representation

### Failure of Krull-Schmidt for invertible lattices over a discrete valuation ring

Let p be a prime greater than 3, and let N be the semi-direct product of a group H of order p by a cyclic C group of order p - 1, which acts faithfully on H. Let R be the localization of Z at p. We

### Generalized GCD matrices

Let f be an arithmetical function. The matrix (f(i,j))n×n given by the value of f in greatest common divisor of (i,j), f (i,j) � as its i, j entry is called the greatest common divisor (GCD) matrix.


text: dc.contributor.advisor Wright, David W. en_US dc.contributor.author Lyon, Robyn Marie dc.date.accessioned 2009-06-29T13:12:05Z dc.date.available 2009-06-29T13:12:05Z dc.date.copyright 2008 en dc.date.issued 2008-05 dc.identifier.other t08026 dc.identifier.uri http://hdl.handle.net/10057/2053 dc.description Thesis [M.A.] - Wichita State University, College of Liberal Arts and Science, Dept. of Sociology en dc.description.abstract The household division of labor and pay inequality between men and women is en examined using the 2006 American Time Use Survey. There are three main theoretical perspectives to examine when discussing household division of labor and pay inequality between men and women. The individual model argues that an individual makes rational choices when investing in Human Capital, which directly affects their income. Structural theory states that the position that an individual occupies determines the income and reward that one receives. Gender theory states that gender is a process of devaluation and sorting that determines the type of jobs that one occupies and the income they receive. It is expected that and increase in household labor will result in a decrease in income. It is also expected that women will have lower income than men. An OLS regression analysis is performed. Women earn $60.40 less per week than men, net of other factors. Increased participation in household labor food preparation results in a slight decrease in income. However, this is only significant for women who receive a decrease in income of$26.62 per week. Policy implications are discussed, these include things such as encouraging women to further their education, reduce the inequality between jobs for men and women, enforcement of antidiscrimination laws based on marital status and an increase in family-friendly polices in the workplace. dc.format.extent vii, 37 leaves, ill. en dc.format.extent 157314 bytes dc.format.mimetype application/pdf dc.language.iso en_US en dc.publisher Wichita State University en dc.subject.lcsh Electronic dissertations en dc.title Household division of labor and pay inequaltiy between men and women en dc.type Thesis en


text:
# Transport and tumbling of polymers in viscoelastic shear flow

16 Mar 2020 Singh Sadhana Singh R. K. Kumar Sanjay

Polymers in shear flow are ubiquitous and we study their motion in a viscoelastic fluid under shear. Employing dumbbells as representative, we find that the center of mass motion follows: $\langle x^2_c(t) \rangle \sim \dot{\gamma}^2 t^{\alpha+2}, ~0< \alpha <1$, generalizing the earlier result: $\langle x^2_c(t) \rangle \sim \dot{\gamma}^2t^3 ~(\alpha = 1)$... (read more)

PDF Abstract

# Code Add Remove Mark official

No code implementations yet. Submit your code now

# Categories

• STATISTICAL MECHANICS
• SOFT CONDENSED MATTER


text:
Acta Mathematica Scientia (Series A)
Sponsored by　Wuhan Institute of Physics and
Mathematics, CAS, China
Edited by  Editorial Committee of Acta Mathematica
Scientia
Add: P. O. Box 71010, Wuhan 430071, China
Tel: 027-87199206(Series A & Series B)
086-27-87199087(Series B)
E-mail: actams@wipm.ac.cn
ISSN 1003-3998
CN 　42-1226/O
 The Solvability of Dual Minkowski Problem in $\mathbb{R}$2 Na Wei Acta mathematica scientia,Series A. 2019, 39 (6):  1314-1322.  Abstract ( 5 )   RICH HTML PDF(306KB) ( 8 )   In this paper, we study the existence of minimum of a constrained variational problem in the Sobolev space W1, 4($\mathbb{S}$). If ∫$_\mathbb{S}$g(θ)dθ>0, the minimum is a positive solution to the related Euler-Lagrange equation ${u^{\prime \prime }} + u = \frac{{g(\theta )}}{{u\left( {{u^2} + {u^{\prime 2}}} \right)}}, \theta \in {\rm{\mathbb{S}}}$ Based on this, we prove the solvability of the dual Minkowski problem in $\mathbb{R}$2 posed by Huang-Lutwak-Yang-Zhang[Acta Math, 2016, 216(2):325-338].
 Traveling Wave Solutions of the Generalized Hyperelastic-Rod Wave Equation Yongyi Gu,Wenjun Yuan,Yonghong Wu Acta mathematica scientia,Series A. 2019, 39 (6):  1342-1351.  Abstract ( 5 )   RICH HTML PDF(385KB) ( 21 )   In this paper, we study the generalized hyperelastic-rod wave equation. We changed the generalized hyperelastic-rod wave equation into a complex differential equation by using traveling wave transform and show that meromorphic solutions of the complex differential equation belong to the class W by the weak $\left\langle {h, k} \right\rangle$ condition and the Fuchs index. Furthermore, we find out all meromorphic solutions of the complex differential equation, then we obtain the traveling wave solutions of the generalized hyperelastic-rod wave equation. We can apply the idea of this study to some related mathematical physics equations.


text:
# A short(ish) introduction to the binomial distribution

Averages make for nice headline figures but are generally of limited use. Probability distributions can tell us far more. For example, if we think a coin is fair and we toss it 30 times and get only 12 tails should we start to question our assumption that the coin is unbiased? What are the chances of that happening? Knowing only that we should expect to see about 15 tails doesn’t tell us very much. Knowing the complete distribution of probabilities (if the coin were known to be fair), on the other hand, will help us to answer that question. In this case – and for many similar problems – we can use the binomial distribution.

The binomial distribution is (I think) one of the more intuitive distributions to get one’s head around. It’s also one of the most important for discrete data. The two parameters needed to describe an instance of the distribution are simple: the number of independent trials (n) and the probability of “success” (p) in each trial (which must be constant from trial to trial). A trial can be tossing a coin and recording whether it lands on tails, throwing a dice and recording whether it is a 1, or any number of other “experiments” with a binary outcome (i.e. “success” or “failure”). (The “success” terminology can be a little confusing since on occasion a success can represent something undesirable happening.)

We can get to the final formula for a binomially distributed discrete random variable (such as the number of tosses out of 30 that come up tails) using some basics of probability. Since the probability of success in any given trial is p and all trials are independent, the probability that the first r trials will lead to r successes is

$p^{r}.$

If we continue to perform trials thereafter until we have completed n trials in total we have to do another nr trials. The probability that any given one of these subsequent trials is a failure is 1-p so the probability they are all failures is

$(1-p)^{n-r}.$

Again, due to independence, the probability of both these things occurring – r successes followed by nr failures – is the product of the two:

$p^{r}(1-p)^{n-r}.$

But this is just one way of getting r successes from n events. nr failures followed by r successes is just as (un)likely. As is any other permutation of r successes and nr failures. So to find the total probability for getting r successes from n trials we need to add up the number of possible ways of getting r successes and nr failures from n trials and multiply by the above probability. The number of ways is given by the binomial coefficient

$\frac{n!}{r!(n-r)!},$

where “!” represents the factorial of the preceding number. If you want to know exactly where this came from then this article explains in detail. The complete formula, the probability that the number of successes X is exactly r, is then

$\textup{Pr}(X=r) = \frac{n!}{r!(n-r)!}p^{r}(1-p)^{n-r}$

for all integer values of r between 0 and n (and 0 otherwise).

The mean, μ, of the binomial distribution is simply

$\mu =np$

while the standard deviation, σ, is given by

$\sigma =\sqrt{np(1-p)}.$

The former may be obvious, the latter probably less so. Detailed derivations can be found here.

These calculations can be tedious and if you’re not careful you can end up with intermediate steps that your calculator cannot cope with thanks to the factorials giving rise to very large numbers. The interactive chart below can hopefully save you some effort (as long as you limit yourself to a thousand trials or less and p isn’t much less than 1 in 100). Drag the sliders to adjust n and p or click/tap on the associated arrows to change these values incrementally. You can also change the value of interest (for the table below the chart) using the slider (or associated buttons) or by clicking inside the rectangular chart region. The blue bars display the binomial distribution for the given parameters while the yellow vertical line marks the value of interest. The red curve will be discussed shortly. The middle column of the table provides details about the binomial distribution and the value of interest.

### Value of interest (r): ←→

Measure Binomial Poisson
Mean (μ)
Std. dev. (σ)
Pr(X = )
Pr(X < )
Pr(X)
Pr(X > )
Pr(X)

Often the trickier part of using the binomial distribution (and distributions in general) is making sure to ask the right question of it. Going back to the coin question – should we worry the coin is biased if we only observe 12 tails from 30 throws? – we might naively look for the probability of getting 12 successes from 30 trials with a fair coin. Using the chart and table we find this has an 8% probability. That’s small but clearly non-negligible. We’d be more suspicious if there were even fewer tails. 12 or fewer tails from 30 tosses is an 18% chance. Wouldn’t we be just as suspicious if we’d got 18 tails and only 12 heads? So there’s a 36% chance that the number of heads differs from the number of tails by at least 6. That’s more than 1 in 3. But we wouldn’t know that if all we looked at was the average.

As already noted, If you’re using a calculator to get values of the binomial distribution you may have a problem due to the inclusion of large factorials. You may be able to work around this using some clever maths. Alternatively, you may be able to use an approximation. One such approximation is the Poisson distribution.

The Poisson distribution concerns the probability of a certain number of independent events happening in a fixed interval when the average rate at which they happen is a known constant. With an average rate per interval of λ, the probability distribution is given by:

$\textup{Pr}(X=r) = \frac{\lambda^{r}e^{-\lambda }}{r!}$

for all integer values of r greater than or equal to 0 (and 0 otherwise). The mean and standard deviation of this distribution are λ and λ½, respectively.

The Poisson distribution deserves a whole article to itself but I’ll try to keep it short here. (You can find out more here and here, for example.) If n is large and p is small then substituting np for λ in the equation above will create a distribution that is a “good” approximation to the binomial distribution. How good an approximation depends on how large and how small n and p are, respectively. Because p is less than 1 (for all non-trivial cases), the factorial in the denominator is smaller than those seen in the binomial coefficient. This may ease use of a calculator.

I’ve added the corresponding Poisson approximation to the chart above (red line) and a third column to the table showing parameters calculated using the approximation. For the case of a coin toss, p is large and therefore the Poisson approximation is poor. On the other hand, if n = 400 and p = 0.04, the approximation is very good. One should still be aware that the Poisson approximation will give a finite probability that the number of successes exceeds the number of trials!

With the interactive tool above the binomial distribution is calculable for all values of n up to 1,000 so there is no particular reason to use the Poisson approximation. However, it does graphically illustrate when the latter may be appropriate next time you get stuck with only a calculator in your hand.