Miscellaneous Topics

---

This chapter deals with miscellaneous NumPy and SciPy functions and techniques. The goal is to present representative examples so you'll be able to search the SciPy documentation more efficiently. The following screenshot shows you where this chapter is headed.

Figure 25: Miscellaneous NumPy Functions Demo

In section 5.1, you'll learn how to use the NumPy searchsorted() binary search function and how to interpret its unusual return value.

In section 5.2, you'll learn how to use SciPy to perform LU decomposition on a square matrix and why decomposition is important.

In section 5.3, you'll learn about NumPy and SciPy statistics functions such as mean() and std().

In section 5.4, you'll learn how to generate random values from a specified distribution such as the Normal or Poisson, and how to bin data using the histogram() function.

In section 5.5, you'll learn about SciPy miscellaneous functions such as the double factorial.

In section 5.6, you'll learn how to use special SciPy functions such as bernoulli() and  gamma().

5.1 Array Binary Search

To search a sorted array, you can use the NumPy searchsorted() function. The searchsorted() function is quite different from the binary search functions in other languages. Also, you must be careful when dealing with arrays that have floating-point values.

Code Listing 22: Array Binary Search Demo

The demo program execution begins by setting up an array to search and a target value to search for:

arr = np.array([1.0, 3.0, 4.0, 6.0, 8.0, 11.0, 13.0])
print "Array arr is "
print arr

target = 11.0
print "Target value to find is " + str(target)

If you need to search a very large array and the array is already sorted, a binary search is often the best approach because it's much faster than a simple sequential search. For small arrays (typically those with less than 100 cells), the marginally faster performance of a binary search is often unimportant, and if your array is not already sorted, the time required to sort the array usually wipes out any time saved by a binary search.

Next, the demo calls the NumPy searchsorted() function like so:

print "Searching array using np.searchsorted() function "
idx = np.searchsorted(arr, target)
if idx < len(arr) and arr[idx] == target:
  print "Target found at cell " + str(idx)
else:
  print "Target not found " 

The binary search functions in most programming languages return a -1 if the target is not found, or return the cell index that holds the target value if the target is found. The searchsorted() function works a bit differently.

A call to searchsorted(arr, x) returns the cell index in sorted array arr where x would be inserted so that the array would remain sorted. For example, if arr = [2.0, 5.0, 6.0, 9.0] and x = 3.0, then searchsorted(arr, x) returns 1 because the 3.0 would be inserted at cell 1 in order to keep the array sorted. If x = 11.0, then searchsorted(arr, x) would return 4 because the 11.0 would have to be inserted beyond the end of the array.

If x is a value that is already in the array, then searchsorted(arr, x) will return the cell where the value is. Therefore, to determine if a value is in an array arr using the return value idx from searchsorted(arr, x), you must first check that idx is less than the length of arr and then check to see if the value at arr[idx] equals the target value.

If the search array holds floating-point values, using searchsorted() is somewhat risky. For example, if the target value is 11.0000000000000001 (there are 15 zeros), it would not be found by the demo program, but a slightly less precise target of 11.000000000000001 (there are 14 zeros) would be found.

The lesson is that, when searching a sorted array of floating-point values using the NumPy searchsorted(), you don't have control over how the function determines floating-point value equality, so you may want to write a program-defined binary search function like my_bin_search() in the demo program:

def my_bin_search(a, t, eps):
  lo = 0
  hi = len(a)-1
  while lo <= hi:
    mid = (lo + hi) / 2
    if np.isclose(a[mid], t, eps):
      return mid
    elif a[mid] < t:
      lo = mid + 1
    else:
      hi = mid - 1
  return -1   # not found

The program-defined function my_bin_search() uses a standard iterative (as opposed to recursive) binary search algorithm with early check for equality, combined with an epsilon parameter to control how close two floating-point values must be for them to be evaluated as equal.

Resources

For details about the NumPy searchsorted() function, see
http://docs.scipy.org/doc/numpy-1.10.0/reference/generated/numpy.searchsorted.html.

For information about the array binary search algorithm used by the demo, see
https://en.wikipedia.org/wiki/Binary_search_algorithm.

5.2 Matrix Decomposition

Matrix decomposition is the process of breaking a matrix down into two smaller matrices that, when multiplied together, give the original matrix (or a slightly rearranged version of the original). An analogy in regular arithmetic is breaking down the number 21 into 3 and 7 because 3 * 7 = 21.

Matrix decomposition, also called matrix factorization, is rarely used by itself, but decomposition is the basis for efficient algorithms that find the inverse and the determinant of a matrix.

There are several kinds of matrix decomposition. The most common form is called lower-upper decomposition for reasons that will become clear shortly. The scipy.linalg.lu() function performs lower-upper matrix decomposition. It's sometimes useful to write a program-defined matrix decomposition function.

Code Listing 23: Matrix Decomposition Demo

The demo program begins by bringing the scipy.linalg submodule into scope:

import numpy as np
import scipy.linalg as spla

After creating the source matrix m and displaying its values, the matrix is decomposed using the linalg.lu() function like so:

print "\nDecomposing m using scipy.linalg.lu() "
(perm, low, upp) = spla.lu(m)

The return result is a tuple with three items. The first item, perm, will be explained shortly. The second and third items are the decomposed matrices. For the demo, return matrix low is:

[[ 1.          0.          0.          0.        ]
[ 0.57142857  1.          0.          0.        ]
[ 0.42857143  0.59090909  1.          0.        ]
[ 0.71428571  0.13636364  1.          1.        ]]

Notice that the relevant values are in the lower part of the matrix, and there are dummy 1.0 values on the main diagonal. The return matrix upp is:

[[ 7.          9.          8.          1.        ]
[ 0.         -3.14285714 -1.57142857 -0.57142857]
[ 0.          0.         -1.5         2.90909091]
[ 0.          0.          0.         -1.54545455]]

Here, all the relevant values are on the main diagonal and above. Next, the demo multiplies low and upp using the NumPy dot() function and displays the resulting matrix:

[[ 7.  9.  8.  1.]
[ 4.  2.  3.  0.]
[ 3.  2.  1.  3.]
[ 5.  6.  4.  2.]]

The original matrix is:

[[ 3.  2.  1.  3.]
[ 5.  6.  4.  2.]
[ 7.  9.  8.  1.]
[ 4.  2.  3.  0.]]

Notice that the product of matrices low and upp is almost the original matrix. Rows 0 and 1 have been swapped and rows 2 and 4 have been swapped. The swap information is contained in the perm matrix result:

[[ 0.  0.  1.  0.]
[ 0.  0.  0.  1.]
[ 1.  0.  0.  0.]
[ 0.  1.  0.  0.]]

This may be interesting, but what's the point? As it turns out, the lower and upper matrices of a decomposition can be used to easily calculate the determinant of the original matrix, and can also be used to compute the inverse of the original matrix.

The determinant of a matrix is the product of the parity of row swaps times the product of the diagonal elements of the upper matrix. The inverse of a matrix can be computed using a short helper function that performs what is called elimination on the lower and upper matrices.

This is exactly how SciPy calculates the determinant and inverse of a matrix. It may seem odd to use such an indirect approach, but decomposing a matrix and then finding the determinant or the inverse is much easier and faster than finding the determinant or inverse directly.

The LU decomposition functions in many other libraries return different values than the scipy.linalg.lu() function. The demo program implements a custom my_decomp() decomposition function that returns values in a different format. The call to my_decomp() is:

print "\nDecomposing m using my_decomp() "
(lum, perm, t) = my_decomp(m)

The program-defined function returns a tuple of three items. The first is a combined lower-upper matrix (instead of separate lower and upper matrices). The second item is a permutation array (instead of a matrix). And the third item is a toggle parity where +1 indicates an even number of row swaps and -1 indicates an odd number of row swaps. For the demo, the combined lower-upper matrix result from my_decomp() is:

 [[ 7.          9.          8.          1.        ]
  [ 0.57142857 -3.14285714 -1.57142857 -0.57142857]
  [ 0.42857143  0.59090909 -1.5         2.90909091]
  [ 0.71428571  0.13636364  1.         -1.54545455]]

These are the same values from linalg.lu() except combined into a single matrix to save space. The result perm array from my_decomp() is:

[2 3 0 1]

This contains essentially the same information as the perm matrix return from linalg.lu(), indicating that if the lower and upper matrices were extracted from the combined LU matrix, and then multiplied together, the result would be the original matrix with rows 0 and 2 swapped and rows 1 and 3 swapped.

Resources

For general information about matrix LU decomposition, see
https://en.wikipedia.org/wiki/LU_decomposition.

For details about the SciPy linalg.lu() decomposition function, see
http://docs.scipy.org/doc/scipy-0.16.0/reference/generated/scipy.linalg.lu.html.

5.3 Statistics

The NumPy and SciPy libraries have a wide range of statistics functions that work with arrays and matrices. Representative examples include the mean(), std(), median(), and corrcoef() functions.

Code Listing 24: Statistics Functions Demo

The demo program execution begins by setting up two parallel arrays. The first array represents the ability scores of six people. The second array represents the pay rates of the six people:

ability = np.array([0., 1., 3., 4., 4., 6.])
payrate = np.array([15., 15., 25., 20., 30., 33. ])

Next, after displaying the values in the two arrays, the demo illustrates the use of the NumPy median() and std() functions:

ma = np.median(ability)
print "The median ability score is "
print ma

s_sd = np.std(payrate, ddof=1)
print "The sample standard deviation of payrates is "
print s_sd

The median is the middle value in an array or, as in this example when there isn't a single middle value, the average of the two values closest to the middle.

By default, the NumPy std() function returns the population standard deviation of its array argument. If you want the sample standard deviation, you can use the ddof (delta degrees of freedom) parameter with value = 1.

Next, the demo computes and displays the Pearson r coefficient of correlation using the corrcoef() function:

pr = np.corrcoef(ability, payrate)
print "Pearson r calculated using np.corrcoef() = "
print pr

The correlation coefficient is a value between -1.0 and +1.0, the magnitude indicating the strength of the linear relation and the sign indicating the direction of the relationship. Notice the output is in the form of a matrix with the coefficient value (0.88700711) duplicated on the minor diagonal.

The demo concludes by calling a program-defined function my_corr() that also calculates the Pearson r coefficient of correlation:

pr = my_corr(ability, payrate)
print "Pearson r calculated using my_corr() = "
print "%1.8f" % pr

There's no advantage to using the program-defined correlation function. The point is that NumPy and SciPy have many built-in statistics functions, but in the rare situations when you need to implement a custom statistics function, NumPy and SciPy have all the tools you need.

Resources

For a list of the NumPy statistics functions, see
http://docs.scipy.org/doc/numpy-1.10.0/reference/routines.statistics.html.

For an explanation of the Pearson correlation coefficient that was used for my_corr(), see
https://en.wikipedia.org/wiki/Pearson_product-moment_correlation_coefficient.

For information about the NumPy corrcoef() function, see
http://docs.scipy.org/doc/numpy-1.10.0/reference/generated/numpy.corrcoef.html.

5.4 Random Numbers

The NumPy library has a wide range of functions that can generate pseudo-random values from a specified distribution type. Representative examples include the random.normal(), random.poisson(), random.exponential(), and random.logistic() functions.

Code Listing 25: Random Sampling Demo

The demo program begins by preparing to generate 100 random values that come from a Normal (also called Gaussian or bell-shaped) distribution with mean = 0.0 and standard deviation = 1.0.

np.random.seed(0)
mean = 0.0
std = 1.0
n = 100

Setting the global random seed, in this case to an arbitrary value of 0, means that the program results will be the same every time the program is run. For a Normal distribution with mean = 0.0, the vast majority of values will be between (-3 * std) and (+3 * std), so we expect all generated values to be in the range [-3.0, +3.0].

Next, the demo program creates an array with 100 cells and fills each cell with a Normal distributed random value:

print "Generating " + str(n) + " Normal values "
values = np.zeros(n)
for i in xrange(n):
  x = np.random.normal(mean, std)
  values[i] = x

An alternative approach is to create the array directly by supplying a value for the optional size parameter: values = np.normal(mean, std, 100). After displaying the 100 values, the demo program constructs histogram information from the values:

bins = 5
print "Constructing histogram data using " + str(bins) + " bins "
(histo, edges) = np.histogram(values, bins=5)

The NumPy histogram() function returns a tuple that has two arrays. The first array stores the count of values in each bin. The second array stores the boundary values for each bin. This is  clearer when you examine the output. The statements:

print histo
print edges

produce the following output:

Count of values in each bin:
[ 6 20 35 27 12]

The beginning and end values of each bin:
[-2.55298982 -1.58844093 -0.62389204  0.34065685  1.30520574  2.26975462]

This means there were 6 values in the interval [-2.55, -1.58), 20 values in [-1.58, -0.62), 35 values in [-0.62, 0.34), 27 values in [0.34, 1.30), and 12 values in [1.30, 2.26]. If you visually scan the 100 values, you can see the smallest value generated is -2.55298982 and the largest is 2.26975462.

The demo program concludes by showing you how to implement a Normal distribution value generator without using NumPy via a program-defined class named Gaussian. The class constructor accepts a mean, a standard deviation, and a seed:

class Gaussian:
  def __init__(self, mean, sd, seed):
    self.mean = mean
    self.sd = sd
    self.rnd = random.Random(seed)

The class uses a Random object from the Python random module. The next() function uses the clever Box-Muller algorithm to transform two uniform random values into one that is Normal.

def next(self):
  two_pi = 2.0*3.14159265358979323846
  u1 = self.rnd.random()  # [0.0 to 1.0)
  while u1 < 1.0e-10:
    u1 = self.rnd.random()
  u2 = self.rnd.random()
  z = math.sqrt(-2.0 * math.log(u1)) * math.cos(two_pi * u2)
  return z * self.sd + self.mean

The while loop in function next() guarantees that variable u1 is not a very small value so that log(u1) won't fail. This example illustrates that it's relatively easy to implement a custom generator in rare situations where NumPy doesn't have the generator you need.

Resources

For a list of NumPy random sampling functions, see
http://docs.scipy.org/doc/numpy-1.10.0/reference/routines.random.html.

For details about the NumPy histogram() function, see
http://docs.scipy.org/doc/numpy-1.10.1/reference/generated/numpy.histogram.html.

For information about the Box-Muller algorithm, see
https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform.

5.5 Double Factorial

The SciPy library has a collection of useful mathematical functions in the scipy.misc submodule. Examples include the misc.derivative(), misc.logsumexp(), and misc.factorial2() functions.

Code Listing 26: Double Factorial Demo

The demo program illustrates the double factorial function, which is best explained by example. The double factorial of n is often abbreviated as n!!, much like n! is an abbreviation for the regular factorial function.

7!! = 7 * 5 * 3 * 1 = 105
6!! = 6 * 4 * 2 = 48

In words, the double factorial is like the regular factorial function except every other term in the product is skipped in the product. The double factorial function is used as a helper in many important mathematical functions such as the specialized gamma function. The demo program begins by importing the scipy.misc submodule:

import scipy.misc as sm

Note that the factorial2() function is also in the scipy.special submodule. After import, the factorial2() function can be called like so:

n = 3
dfact = sm.factorial2(n)
print "Double factorial of " + str(n) + " using misc.factorial2() = "
print str(dfact)

The factorial2() function has an optional parameter exact that, if set to False, allows the function to do a fast approximation rather than a slower exact calculation.

The demo implements a program-defined version of the double factorial function named my_double_fact(). There's no advantage to a program-defined version unless you need some sort of specialized behavior, or wish to avoid importing a module for some reason.

Resources

For details about the misc.factorial2() function, see
http://docs.scipy.org/doc/scipy-0.16.1/reference/generated/scipy.misc.factorial2.html.

For information about the double factorial function, including alternate definitions, see
https://en.wikipedia.org/wiki/Double_factorial.

5.6 The Gamma Function

The SciPy library has a large collection of mathematical functions in the scipy.special submodule. Examples include elliptic functions, Bessel functions, advanced statistical functions, and gamma functions.

Code Listing 27: Gamma and Special Gamma Demo