Calculating "generating functions" with numpy

Question

In mathematics, a "generating function" is defined from a sequence of numbers c0, c1, c2, ..., cn by c0+c1*x+c2*x^2 + ... + cn*x^n. These come as "moment generating functions", "probability generating functions" and various other types, depending on the source of the coefficient.

I have an array of the coefficients and I'd like a quick way to create the corresponding generating function.

I could do

import numpy as np
myArray = np.array([1,2,3,4])
x=0.2
sum([c*x**k for k,c in enumerate myArray])

or I could have an array having c[k] in the kth entry. It seems there should be a fast numpy way to do this.

Unfortunately attempts to look this up are complicated by the fact that "generate" and "function" are common words in programming, as is the combination "generating function" so I haven't had any luck with search engines.

You are only interested in evaluating the function for one value of ```x``` at a time? — wwii, Aug 20 '15 at 04:08
@JohnColeman Good question - suppose it's because I'm actually dealing with truncations of infinite sums, so it didn't occur to me. — Joel, Aug 20 '15 at 06:08

wwii · Accepted Answer · 2015-08-22T14:10:44.643

x = .2
coeffs = np.array([1,2,3,4])

Make an array of the degree of each term

degrees = np.arange(len(coeffs))

Raise x the each degree

terms = np.power(x, degrees)

Multiply the coefficients and sum

result = np.sum(coeffs*terms)

>>> coeffs
array([1, 2, 3, 4])
>>> degrees
array([0, 1, 2, 3])
>>> terms
array([ 1.   ,  0.2  ,  0.04 ,  0.008])
>>> result
1.552
>>>

As a function:

def f(coeffs, x):
    degrees = np.arange(len(coeffs))
    terms = np.power(x, degrees)
    return np.sum(coeffs*terms)

Or simply us the Numpy Polynomial Package

from numpy.polynomial import Polynomial as P
p = P(coeffs)
result = p(x)

score 2 · Answer 2 · answered Aug 20 '15 at 04:58

2

If you are looking for performance, using np.einsum could be suggested too -

np.einsum('i,i->',myArray,x**np.arange(myArray.size))

answered Aug 20 '15 at 04:58

Divakar

218,885
19
262
358

@Joel This [`Why is numpy's einsum faster than numpy's built in functions?`](http://stackoverflow.com/q/18365073/3293881) might give some insights on performance benefits with `np.einsum`. Though, in this case, since you are summing up elements only along one axis, the performance benefits over using `np.sum` might be marginal. – Divakar Aug 20 '15 at 06:37

score 1 · Answer 3 · answered Aug 20 '15 at 04:30

>>> coeffs = np.random.random(5)
>>> coeffs
array([ 0.70632473,  0.75266724,  0.70575037,  0.49293719,  0.66905641])
>>> x = np.random.random()
>>> x
0.7252944971757169
>>> powers = np.arange(0, coeffs.shape[0], 1)
>>> powers
array([0, 1, 2, 3, 4])
>>> result = coeffs * x ** powers
>>> result
array([ 0.70632473,  0.54590541,  0.37126147,  0.18807659,  0.18514853])
>>> np.sum(result)
1.9967167252487628

score 1 · Answer 4 · answered Aug 20 '15 at 04:44

Using numpys Polynomial class is probably the easiest way.

from numpy.polynomial import Polynomial

coefficients = [1,2,3,4]
f = Polynomial( coefficients )

You can then use the object like any other function.

import numpy as np
import matplotlib.pyplot as plt

print f( 0.2 )

x = np.linspace( -5, 5, 51 )
plt.plot( x , f(x) )

Calculating "generating functions" with numpy

4 Answers4