Mathematics and Computing
Getting Started with Programming </h1>Dr. Jeff Knisley
In the first notebook, we introduced the Python language and used it to explore how mathematics can be implemented and advanced via algorithms. In this notebook, we begin to explore the relationship between algorithms and programs, as most people tend to think of "programming a computer" rather than "equipping a computer with algorithms."
The difference is both subtle and obvious. An algorithm is a sequence of steps designed to solve a problem, or equivalently, there are mathematical problems (e.g., estimating the square root of a positive number) whose solution is in the form of a sequence of steps. What is obvious is that a program is an implementation of an algorithm -- usually on a machine, most often a computer. What is subtle is that a program includes both the algorithm and the types of objects the algorithm is applied to.
This subtle yet obvious relationship is analogous to a cooking recipe. Although we think of the recipe as the sequence of steps that leads to the desired food dish, the types (and amounts) of ingredients at the top of the recipe is as necessary as the sequence of steps. Of course, substitutions are possible, and some ingredients ( a teaspoon of nutmeg ) are more specific than others (season to taste with some type of salt).
Similarly, we must convert algorithms into recipes implemented on a computer, which we tend to call programs. Let's first configure the notebook and then we'll explore in more detail the types of ingredients (data structures) included in Python.
Note: We are using Python 3 for this assignment, so no need to import from the __future__ (Python 3 is the future!)
Strings¶
As we saw in the first notebook, a string is a special type of list, one that is made solely of alphanumeric characters (and punctuation). (Please review first notebook if you do not remember what a string is). Strings are defined by enclosure in either single quotes ' ', double quotes " ... ", or triple quotes, ''' ... ''' or """ ... """.
Characters in a string are accessed using an index inside square brackets, beginning at 0. Thus, if
MyString = ' Hello World '
then MyString[0] returns the string ' ', which is the whitespace character that begins the string. To further explore this example, let's notice that len(MyString) is 22, and that we can index each character as follows:
1 2
0123456789012345678901
MyString = ' Hello World! '
Thus, MyString[5] returns the string 'H', and MyString[5:10] returns the string "Hello".
MyString = ' Hello World! '
MyString
MyString[0]
MyString[5]
MyString[5:10]
There are also a number of methods associated with strings. If a string is bound to a name such as MyString, then
dir(MyString)
produces a list of the methods that can be applied to the string.
dir(MyString)
For example, we can strip the spaces off of either end of MyString via MyString.strip().
MyString.strip() # removes leading and trailing spaces
The result does not change MyString, but instead creates a copy of the string with the leading and trailing spaces removed. Likewise, the methods attached to the MyString object produce a new string object, while the original MyString object remains unchanged.
MyString.upper()
MyString.find('lo')
MyString.split(' ') # ' ' is a space character
MyString.replace("World","Universe")
Casting and Comprehension¶
Obviously, a string is not an integer, and a float (decimal) is not a string. Likewise, a list is not a tuple and a tuple is not a string and so on. That is, each object has its own type, which can be revealed using the type command.
type("Hello World")
type( 1.2345 )
type( [1,2,3,4] )
Casting an object is the process of creating a new object from the old but with a different type. For example, we may want to use a tuple as a list, or we might want to use a decimal as a string. To do so, we simply apply the command for that type.
str(1.2345)
list( (1,2,3,4) )
It is important to note casting creates a new object. For example, below we bind the number 1.2345 to a and then cast it to a string, which we call b. This does not change the object bound to a!
a = 1.2345
b = str( a )
a
type(a)
b
type(b)
In contrast, comprehension is the concept of creating a new object using the content of a given object. That is, whereas casting an object is a new object with a different, comprehension is a new object with modified contents.
For example, list comprehension is the act of making a new list from any iterable object in the form
[ f(item) for item in IterableObject ]
where f(item) is some function or manipulation that is applied sequentially to each item in the object.
For example, the next few cells change a number into a list of its digits.
number = 1235
[ digit for digit in str(number) ]
That is, we cast a number into a string using str(number), and then we use comprehension to create a list from the individual letters in the string.
Moreover, if we so desired, we could convert each digit string into a numerical digit by casting via the command int.
[ int(digit) for digit in str(number) ]
To see how useful list comprehension is, let's recall that dir(MyString) produced a huge list of properties and methods attached to MyString, but in general, we are not interested in the ones that begin with an underscore _ because they tend to be used internally in Python rather than by programmers.
dir(MyString)
We can use list comprehension, however, to obtain only the names we are interested in! To do so, we add a conditional statement which says that we want names that do not begin with a _, which in Python is
name[0] != '_'
where != is Python for "not the same as".
[ name for name in dir(MyString) if name[0] != '_' ]
In fact, we can include only those names that begin with a given first letter. For example,
name[0] == 's'
means we only want names associated with MyString that begin with the letter 's'.
IMPORTANT NOTE!!!!!
The logical concept of equality as "the same as"
is denoted by 2 equal signs, ==.
This is because one equal sign alone is used to assign values to a name.
[ name for name in dir(MyString) if name[0] == 's' ]
Casting and comprehension are extremely useful, so much so that new features of Python and even entirely new langages (e.g., Julia) are based on the importance of types, casting, and comprehension.
Dictionaries¶
In addition to lists and tuples, a very useful type of object is the dictionary object. In particular, the command dict creates a Python dictionary, which is a set (order does matter) of key-value pairs of the form
{ key1:value1, key2:value2, … }
For example, dict( a = 5, b=7, c=3 ) creates the dictionary object
{ 'a': 5, 'b': 7, 'c':3 }
Below are additional examples of dictionary uses in Python.
MyDictionary = dict( a = 5, b=7, c=3 ) # bind to the name MyDictionary
MyDictionary
MyDictionary['a']
MyDictionary['c']
MyDictionary['b'] = 11
MyDictionary
We can also get the keys and values from a dictionary separately.
MyDictionary.keys() #returns the keys for the dictionary
MyDictionary.values() #returns the keys for the dictionary
Alternately, a dictionary can be defined as a set of key-value pairs. However, notice also that like a set, the order in which the values are entered is not preserved. Instead, Python chooses the order, but all the key-value pairs are maintained !!
{ 'a':2, 'c':7, 11:14 }
_['a']
Dictionaries are also quite useful, as we are about to discover.
Putting Python to work¶
We could continue in this fashion indefinitely, but ultimately, it distracts from the reason that computers languages exist -- namely, because people want computers to do things that they and others will find useful and exciting!
Instead, as the experiences of programmers and scientists the world over has shown repeatedly, one of the better ways to learn a new language is to
have a problem to solve and
hack a similar approach until it works for us (where "hack" means modifying until suitable for the task at hand).
In fact, our entire usage of Python in this course is motivated by
A Mathematical problem to solve that requires a computational approach.
An approach to the problem, which you the student may be asked to "hack" in some straightforward fashion.
Thus, we will motivate our initial foray into programming with Python via a problem to be solved. Our problem is to use computers to analyze strings. Sometimes, this is called text mining, which broadly is the problem of
__ _counting the number of occurences of a given pattern in a string._ __
For example, this is one of the most important computational problems in genomics, because a DNA molecule is nothing more than a (typically very very long) string which only includes the letters A, C, G, and T, which represent base pairs.
DNAsnippet = "AGTCGAATCGACAAAATCG"
Although DNA molecules tend to have millions or billions of base pairs, we should always remember one of the fundamental principles of programming:
Use a full range of examples to TEST the implementation of an algorithm as a program .
Thus, we start small with a string that will test our program.
Patterns are Regular Expressions¶
In the text mining mindset, a pattern is a special type of string known as a regular expression. Just as $\LaTeX$ has over many decades of development become the lingua franca of textual scientific communication, so also many decades of development have produced the regular expressions format for finding patterns in (often extremely large) strings. In Python, this means importing the re package.
import re
Regular expressions can themselves solve most simple pattern matching problems. We will see this using the findall method of re, which does exactly what the name implies.
re.findall('pattern','String with one pattern or many patterns')
Thus, we can now "find all" occurences of a given pattern in a given string, such as our DNAsnippet.
re.findall("CG", DNAsnippet)
re.findall("TA", DNAsnippet)
So clearly, the pattern CG occurs twice in DNAsnippet, whereas the pattern TA does not occur at all in
DNAsnippet = "AGTCGAATCGACAAAATCG"
However, we can do much more. The reason patterns are called regular expressions is that they can find a variety of patterns, such as all the 3 letter substrings that begin with the letter A. To do so, we place a period (".") in a place where any letter can occur (that is, the period is a wild card character -- whereas to search for a period, we use "\.").
re.findall("A..", DNAsnippet)
Other characters also have special meaning. For example, the + sign means "one or more" of whatever it follows. Thus, the pattern 'A+' means one or more repetitions of the letter A.
re.findall("A+", DNAsnippet)
To find 2 or more repetitions of the letter A, we simply use "AA+", which means the letter A followed by one or more repetitions of the letter A.
re.findall("AA+",DNAsnippet)
More information on using regular expressions can be found at
https://docs.python.org/2/howto/regex.html and also https://en.wikibooks.org/wiki/Python_Programming/Regular_Expression
Indeed, there are many tutorials available online.
Exercise 1¶
The shortest viral genomes are those in the virus family </i>Circoviridae</i>. The string below is the complete genome of the Porcine Cirocovirus and is 1,767 basepairs (i.e., letters) in length.
PorcineCirocovirusGenome = "ACCAGCGCACTTCGGCAGCGGCAGCACCTCGGCAGCACCTCAGCAGCAACATGCCCAGCAAGAAGAGTGGA"+\
"AGAAGCGGACCCCAACCACATAAAAGGTGGGTGTTCACGCTGAATAATCCTTCCGAAGACGAGCGCAAGAAAATACGGGAGCTCCCAATCTCCCTATT"+\
"TGATTATTTTATTGTTGGCGAGGAAGGTAATGAGGAGGGCCGAACACCTCACCTACAGGGGTTCGCTAATTTTGTGAAGAAGCAAACTTTTAATAAAG"+\
"TGAAGTGGTATTTTGGTGCCCGCTGCCACATCGAGAAAGCGAAAGGAACAGATCAGCAGAATAAAGAATATTGCAGTAAAGAAGGCAACTTACTGATA"+\
"GAATGTGGAGCTCCTAGATCTCAAGGACAACGGAGTGACCTCTCTACTGCTGTGAGTACCTTGTTGGAGAGCGGGAGTCTGGTGACCGTTGCAGAGCA"+\
"GCACCCTGTAACGTTTGTCAGAAATTTCCGCGGGCTGGCTGAACTTTTGAAAGTGAGCGGGAAAATGCAGAAGCGTGATTGGAAGACGAATGTACACG"+\
"TCATTGTGGGGCCACCTGGGTGTGGCAAAAGCAAATGGGCTGCTAATTTTGCAGACCCGGAAACCACATACTGGAAACCACCTAGAAACAAGTGGTGG"+\
"GATGGTTACCATGGTGAAGAAGTGGTTGTTATTGATGACTTTTATGGCTGGCTGCCGTGGGATGATCTACTGAGACTCTGTGATCGATATCCTTTGAC"+\
"TGTTGAGACTAAAGGTGGAACTGTACCTTTTTTGGCCCGCAGTATTCTGATTACCAGCAATCAGACCCCGTTGGAATGGTACTCCTCAACTGCTGTCC"+\
"CAGCTGTAGAAGCTCTCTATCGGAGGATTACTTCCTTGGTATTTTGGAAGAATGCTACAGAACAATCCACGGAGGAAGGGGGCCAGTTCGTCACCCTT"+\
"TCCCCCCCATGCCCTGAATTTCCATATGAAATAAATTACTGAGTCTTTTTTATCACTTCGTAATGGTTTTTATTATTCACTTAGGGTTAAGTGGGGGG"+\
"TCTTTAAGATTAAATTCTCTGAATTGTACATACATGGTTATACGGATATTGTAGTCCTGGTCGTATATACTGTTTTCGAATGCAGTGCCGAGGCCTAC"+\
"ATGGTCTACATTTCCAGTAGTTTGTAGTCTCAGCCAGAGTTGATTTCTTTTGTTATTGGGTTGGAAGTAATCGATTGTCCTATCAAGGACAGGTTTCG"+\
"GGGTAAAGTACCGGGAGTGGTAGGAGAAGGGCTGGGTTATGGTATGGCGGGAGGAGTAGTTTACATAGGGGTCATAGGTTAGGGCATTGGCCTTTGTT"+\
"ACAAAGTTATCATCTAGAATAACAGCAGTGGAGCCCACTCCCCTGTCACCCTGGGTGATTGGGGAGCAGGGCCAGAATTCAACCTTAACCTTTCTTAT"+\
"TCTGTAGTATTCAAAGGGCACAGTGAGGGGGTTTGAGCCCCCTCCTGGGGGAAGAAAATCATTAATATTAAATCTCATCATGTCCACATTCCAGGAGG"+\
"GCGTTCTGACTGTGGTTTTCTTGACAGTATAACCGATGGTGCGGGAGAGGCGGGTGTTGAAAATGCCATTTTTCCTTCTCCAGCGGTAACGGTGGCGG"+\
"GGGTGGACGAGCCAGGGGCGGCGGCGGAGGATCTGGCTAAGATGGCTGCGGGGGCGGTGTCTTCGTCTGCGGAAACGCCTCCTTGGATACGTCATCGC"+\
"TGAAAACGAAAGAAGTGCGCTGTAAGTATT"
Enter the pattern which would allow you to find all occurences in the PorcineCirocovirusGenome in which there is a single letter between two occurences of the letter A.
re.findall("A.A",PorcineCirocovirusGenome)
The shortest non-viral genome is more than 500,000 letters (basepairs) long (for a certain bacterium), so we won't go in that direction. Instead, we will explore some simple text mining applications.
A Simple Text Mining Application¶
In some sense, regular expressions make all basic problems rather easy. Thus, let's move on to a slightly more general problem -- one that is performed daily on millions of web pages by hundreds of different companies and organizations (e.g., Google).
The Problem:
Given a document as a list of words, count the number
of words in the document that satisfy a given pattern.
We will solve this problem using Python dictionaries, because key-value pairs are ideal for counting the number of times a pattern occurs in a string. For example, execute the cell belows:
TargetString = "There is a tree near the street."
PatternCounts = dict()
re.findall('ee',TargetString) #all occurences of the pattern 'ee'
PatternCounts['ee'] = len(re.findall('ee',TargetString))
PatternCounts
re.findall('aa', TargetString) #all occurences of the pattern 'aa'
PatternCounts['aa'] = len(re.findall('aa',TargetString))
PatternCounts
However, we would not want to use findall to determine how many times the word 'tree' occurs in the TargetString, because it would not only find the word 'tree', but it would also find the occurence of 'tree' inside of the word street.
re.findall('tree',"There is a tree near the street.")
Instead, we split the string at the spaces to create a list of words, after which we use a regular expression match object.
If, Else, and Match¶
The re.match command is used in the form
re.match('pattern',TargetString)
If the pattern matches the beginning of the TargetString, then a match object is returned. Otherwise, if there is no match, then the value returned is None.
re.match('tree','trees')
re.match('tree','street')
We ofen use re.match in conjunction with a Python if statement, which is of the form
if condition_to_be_tested :
First Command to execute if condition is True
Second Commmand to execute if condition is True
Again, whitespace (i.e., indentation) is used to show what is to be done if the condition is True. Try it out below.
REMINDER Python is case sensitive -- True must have a capital T
condition = True # Change to True to see what happens differently
if( condition) :
print("The condition variable now has the value True")
We can also include an else statement to be executed when the condition is False. This leads to if statements in the form
if condition_to_be_tested :
First Command to execute if condition is True
Second Commmand to execute if condition is True
else:
First Command to execute if condition is False
Second Commmand to execute if condition is False
condition = False # Change to True to see what happens differently
if( condition) :
print("The condition variable now has the value True")
else:
print("The condition has the value False")
When combined with re.match, it is possible to determine when two "words" begin with the same pattern.
if( re.match('tree','tree.') ):
print("True")
else:
print("False")
A None value produces no output. However, in logical terms, None has a value of False
if( re.match('tree','street') ):
print("True")
else:
print("False")
Thus, we need only split a document into individual words and count the number of matches of a given pattern.
Pattern Counting in Strings¶
With all these data types, the program (i.e., 'recipe') for counting patterns becomes rather straightforward. If given a Document and a list of strings as Patterns, do the following:
- Convert the Document into a WordList.
- Create a dictionary called PatternCounts whose keys are the Patterns and whose initial values are 0.
- For each word in WordList
- Compare to each pattern in PatternCounts .
- If a match, then increase the value of that pattern by 1
The resulting dictionary will then be the counts of the occurences of the patterns in the document.
For example, suppose we have the poem "Trees" by Joyce Kilmer (an American poet and soldier who died in battle in WW1), which we load as a triple quote in the next cell.
TreesByJoyceKilmer = """\
I THINK that I shall never see
A poem lovely as a tree.
A tree whose hungry mouth is prest
Against the sweet earth's flowing breast;
A tree that looks at God all day,
And lifts her leafy arms to pray;
A tree that may in summer wear
A nest of robins in her hair;
Upon whose bosom snow has lain;
Who intimately lives with rain.
Poems are made by fools like me,
But only God can make a tree."""
The Joyce Kilmer Memorial Forest in Western North Carolina is an old growth forest containing many, many trees, but how many time does the word "tree" occur in the poem it memorializes? To answer that, we want to count the number of times the pattern 'tree' occurs as a word in the poem, once we have converted the poem to all lower case.
treesbyjoycekilmer = TreesByJoyceKilmer.lower()
treesbyjoycekilmer
Now we must convert the poem into a list of words, which we do using the string method split.
WordList = treesbyjoycekilmer.split()
WordList[:10] #The first 10 words in the list
Matching is not as easy as we may think, however, as revealed by the last 10 words in WordList.
WordList[-10:] #The last 10 words in the list
The problem is that the "tree." has a period in it, but re can handle it easily. The match method defines a match to be when a string is the same as a pattern (regular expression) for the length of the regular expression.
Thus, we can continue, trusting that re.match will handle any subtleties of matching. Next, we define a dictionary with the pattern "tree" (lower case since we converted the poem to lower case).
PatternCounts = { 'tree':0 }
Now it is simply a matter of checking each word in the WordList to see if it matches the pattern. To do so, we do the following:
- for each word in the WordList
- compare the word to each of the keys in PatternCounts
- if a word matches a pattern, then increase the count for that pattern by 1
That is, the matching part of the algorithm requires the use of a Python if statement. Specifically, the if statement executes its commands if there is a match, but does not if None (i.e., no match) occurs.
Note: We repeat the PatternCounts statement since initialization is an important part of the algorithm.
PatternCounts = { 'tree':0 }
for word in WordList:
for pattern in PatternCounts.keys():
if( re.match(pattern,word) ):
PatternCounts[pattern] += 1
PatternCounts
The += operator changes the value of the object by increasing it by the increment, which in this case means adding 1 and thus increasing the count of the number of matches.
Exercise 2¶
The algorithm above is repeated below. Modify it so that it counts the number of times the word that occurs in the poem.
PatternCounts = { 'tree':0 }
for word in WordList:
for pattern in PatternCounts.keys():
if( re.match(pattern,word) ):
PatternCounts[pattern] += 1
PatternCounts
PatternCounts = { 'that':0 }
for word in WordList:
for pattern in PatternCounts.keys():
if( re.match(pattern,word) ):
PatternCounts[pattern] += 1
PatternCounts
A Simple Text Mining Application¶
The algorithm we just developed applies not only to short poems, but also to much larger documents. For example, when you uploaded this notebook, you should have also uploaded a text file with the name
USconstitution.txt
If so, then the next cell will show it as one of the files in the files list.
ls
The next cell uses a Python context ( the with structure) to read the USconstitution into Python in as a triple quote. A context allows Python to do "cleanup" (e.g., close the file) if anything strange happens it is not anticipating.
with open('USconstitution.txt','r') as TheFile:
USconstitution = TheFile.read()
USconstitution[:100]
Let's first do some pre-processing. As you can see in the first 100 characters listed above, the constitution doesn't follow modern conventions for capitalization. Also, there are characters (e.g., \n is the newline commmand) that have been inserted that were not in the original. We thus reduce all to lower case and remove strange characters.
usconstitution = USconstitution.lower()
# replace anything not a letter or number by a space
constitution = re.sub('[^a-z0-9]',' ',usconstitution)
constitution[:100]
Now we transform the constitution into a list of words.
WordList = constitution.split()
WordList[:10] #The first 10 words in the list
And now we can search for patterns. There are probably not any "trees" in the constitution, but words like "president" and "congress" are important!
PatternCounts = { 'president':0, 'congress':0 }
for word in WordList:
for pattern in PatternCounts.keys():
if( re.match(pattern,word) ):
PatternCounts[pattern] += 1
PatternCounts
The president gets more mention than congress, but there is an issue. The word "senate" also refers to congress.
PatternCounts = { 'president':0, 'congress':0, 'senate':0 }
for word in WordList:
for pattern in PatternCounts.keys():
if( re.match(pattern,word) ):
PatternCounts[pattern] += 1
PatternCounts
And what about the Judiciary? That is your job!
Exercise 3¶
Copy the previous cell's code into the next cell and add more search terms to PatternCounts. The terms you add should be justification for how balanced the constitution is in its discussion of the executive, legislative, and judicial branches -- at least as evidenced by word counts from the constitution.
PatternCounts = { 'president':0, 'congress':0, 'senate':0, 'executive':0, 'legislative':0, 'judicial':0 }
for word in WordList:
for pattern in PatternCounts.keys():
if( re.match(pattern,word) ):
PatternCounts[pattern] += 1
PatternCounts
First Letter Distributions¶
Amazingly, armed with nothing more than regular expressions and our simple dictionary-based algorithm, we can do all sorts of "text mining" types of applications, including visualizations!
First, however, let's set up the notebook for inline plots and load the pyplot module with the name plt.
%matplotlib inline
from matplotlib import pyplot as plt #Matlab like plotting
That is, once we have converted a document into words -- i.e., into individual 'strings' -- then we can count occurences of patterns we may be interested in.
For example, we could count the number of words that begin with a vowel, since re.match returns a match if pattern occurs at the beginning of a word.
PatternCounts = { 'a':0, 'e':0, 'i':0, 'o':0, 'u':0 }
Document = """Attention! Attention! This is a test of the emergency
notification system. If this had been an actual emergency,
we would not have said this is a test of the emergency
notification system. We would have said
'Attention! Attention! This is an emergency'.
This concludes our test of the emergency notification system."""
WordList = Document.lower().split()
for word in WordList:
for pattern in PatternCounts.keys():
if( re.match(pattern,word) ):
PatternCounts[pattern] += 1
PatternCounts
More importantly, if given a list of numbers or a document containing numbers, we can count how often the first digit of a number in the list is a 1, 2, 3, 4, 5, 6, 7, 8, or 9, respectively.
For example, the list below contains distances in light years to the 300 Brightest Stars.
LightYearsTo300BrightestStars = [9,310,4,37,25,42,770,11,144,430,530,17,320,65,600,260,34,25,3000,350,78,
430,52,700,88,240,130,111,1300,101,820,81,590,124,840,1800,145,101,270,630,
82,420,105,180,80,500,180,430,126,66,96,220,61,97,200,720,127,170,47,93,360,
36,610,130,1400,690,75,570,78,1500,230,148,920,55,380,400,65,550,310,79,210,
670,460,77,84,84,200,3000,49,540,72,140,220,460,390,58,530,1300,400,165,89,
160,73,60,77,270,170,140,99,37,520,510,116,310,520,1100,310,460,170,88,440,
39,1300,59,82,89,150,360,360,110,570,160,35,24,430,77,330,63,40,980,600,240,
370,102,39,71,170,230,180,200,160,440,170,460,540,730,390,610,260,1600,215,
185,88,760,220,420,250,96,132,1800,83,124,147,480,820,205,530,340,2600,300,
2000,310,85,340,900,250,100,149,150,139,410,101,86,48,570,78,540,220,240,370,
420,44,340,230,220,53,26,225,86,127,151,45,510,108,62,600,99,285,185,101,290,
630,215,101,560,160,185,370,102,175,500,530,350,380,1100,190,325,81,72,155,
120,165,180,1300,440,900,57,185,50,500,73,135,1100,740,200,730,150,235,370,
475,116,47,210,64,495,27,140,420,260,125,135,880,19,118,390,117,82,2100,112,
12,1200,420,220,250,130,790,59,115,270,117,29,135,1900,370,97,155,290,130,155,105
]
We want to count based on the first digit of each number. That is, we want to use the numbers as strings of digits , and consequently, we use list comprehension to transform each number into a string.
NumbersAsStrings = [ str(number) for number in LightYearsTo300BrightestStars ]
NumbersAsStrings[:10] #the first 10 numbers
Now we simply count using the same approach as above (we don't care about the 0 digit, so it is omitted).
PatternCounts = { '1':0, '2':0, '3':0, '4':0, '5':0,
'6':0, '7':0, '8':0, '9':0 }
for word in NumbersAsStrings:
for pattern in PatternCounts.keys():
if( re.match(pattern,word) ):
PatternCounts[pattern] += 1
PatternCounts
Notice that the digit 1 occurs more frequently than the other digits. We can in fact visualize this using a column graph in pyplot.
numdigits = [1,2,3,4,5,6,7,8,9]
strdigits = ['1', '2', '3', '4', '5', '6', '7', '8', '9']
Counts = [ PatternCounts[digit] for digit in strdigits ]
plt.bar( numdigits, Counts, align='center')
plt.xticks( numdigits ) # make sure all the digits are on the x-axis
plt.xlabel("Digits")
plt.ylabel("Counts");
That '1' occurs more often than '2', and '2' more often than '3', etcetera, is a well-known phenomenon known as Benford's law, which we will explore next.
Benford's Law¶
Benford's law says that in "naturally occurring" random numbers, the first digit of a number is 1 more often than it is 2, 3, 4, 5, 6, 7, 8, or 9. In particular, in base 10, the digit $d$ occurs with probability
$$\Pr(d) = \log_{10}(d+1) - \log_{10}(d) $$.
In the Python numerical package numpy, the base 10 logarithm is log10(x)
import numpy
numpy.log10(1000)
Since we will be using the logarithm rather often, let's import it itself so we don't always have to type "numpy" in front of it.
from numpy import log10
We can now use list comprehension to create the Benford's law distribution. In particular,
$$Pr(d) = \log_{10}\left( \frac{d+1}d \right) = \log_{10}\left(1+\frac 1 d \right) $$Let's try it out (remember, that from the \_\_future\_\_ we need to have imported division).
[ log10(1+1/d) for d in range(1,10) ]
Now notice how similar this is to the PatternCounts.values() divided by the number of stars (300).
[ count/300 for count in Counts ]
Not exactly, but the same type of falling off of frequency after the digit 1 which occurs about 30% of the time.
Benford's law is a phenomenon closely related to how numbers are represented in base 10 (or more generally, in base b). It is often "explained" in terms of how numerical values depend on units. In particular, "naturally occurring" random numbers tend to be associated with units like feet, inches, meters, kilograms, etcetera. A change of units implies a change in the numerical quantity. For example, 1 inch = 2.54 centimeters.
However, scaling a quantity implies a carry law. For example, multiplying by 2 is the same as adding a number to itself, and because of the possible final "carry" of 1, the digit 1 is more likely than any other digit to be the first digit of the result -- scaling 5.2 by 2 yields 10.4, which has a first digit of 1 because of the carry.
Let's look at another example, which are populations of the countries in the world in 2010 from largest to smallest.
CountriesPopulations = [1360720000,1240380000,317541000,249866000,201032714,185649000,173615000,152518015,143657134,
127220000,118395054,99133000,89708900,86613986,85964400,80586000,77217000,76667864,67514000,65926261,
65820916,63705000,59917907,53259000,52981991,50219669,47464000,46704314,45439822,44928923,44354000,
40117096,38700000,38502396,37964000,35357000,35295770,34035000,33174000,30475144,30183400,30004000,
29994272,28946101,26494504,25500100,25235000,24895000,24658823,23700715,23382836,23373517,23202000,
21898000,21263403,20609294,20386799,20277597,20121641,17322796,17165000,17129076,16838100,16634603,
16363000,15679000,15438384,15302000,15135000,14580290,13567338,12973808,12825000,11296000,11184873,
11167325,10886500,10824200,10815197,10537222,10513800,10496000,10487289,10413211,10323000,10163000,
10027254,9906000,9639741,9477100,9468100,9445281,8555072,8501502,8264070,8134100,8112200,8104000,
7282041,7184000,7181505,7059653,6783374,6580800,6547800,6340000,6333000,6202000,6191155,6190280,
6071045,5663133,5627235,5450614,5415459,5399200,5240000,5096300,4822000,4667096,4616000,4593100,
4483800,4512900,4448000,4420549,4294000,4290612,3957000,3791622,3615086,3559500,3461041,3405813,
3286314,3065850,3017000,2944459,2821977,2931300,2711476,2113077,2074000,2062294,2061644,2045239,
2024904,2005200,1849000,1815606,1704000,1672000,1622000,1328019,1311870,1257900,1250000,1234571,
1066409,873000,865878,858038,840974,784894,743798,744160,620029,598200,581344,567000,537000,
534189,491875,416055,405739,392291,393162,351461,349728,325620,317280,294906,285000,268270,
264652,258958,237549,212645,187820,187356,180000,159358,150563,109000,106461,106405,103328,
103036,101484,101351,99000,90945,86295,84497,76098,71293,64237,63085,56483,56086,55519,55456,
54000,53883,48244,37429,36979,36942,36136,33540,31458,30001,29537,28502,23296,20901,14974,
13452,13135,11323,9945,8938,6081,4922,4000,2655,2563,2302,2072,1613,1411,839,596,56]
Now we convert the numbers into strings...
NumbersAsStrings = [ str(number) for number in CountriesPopulations ]
NumbersAsStrings[:5] #first 5
... and apply our algorithm to NumbersAsStrings.
PatternCounts = { '1':0, '2':0, '3':0, '4':0, '5':0,
'6':0, '7':0, '8':0, '9':0 }
for word in NumbersAsStrings:
for pattern in PatternCounts.keys():
if( re.match(pattern,word) ):
PatternCounts[pattern] += 1
PatternCounts
We divide each value in PatternCounts by the total, thus producing a Frequency distribution.
This is shown in the column graph below:
numdigits = [1,2,3,4,5,6,7,8,9]
strdigits = ['1', '2', '3', '4', '5', '6', '7', '8', '9']
TotalCount = len( NumbersAsStrings )
Frequencies = [ PatternCounts[digit]/TotalCount for digit in strdigits ]
plt.bar( numdigits, Frequencies, align='center')
plt.xticks( numdigits ) # make sure all the digits are on the x-axis
plt.xlabel("Digits")
plt.ylabel("Frequencies")
plt.title("Does resemble Benford's law distribution");
Benford is not Always true!¶
Benford's law does not apply to all "naturally occurring" random numbers, however. For example, heights of human males in units of decimal feet mostly have either a 5 or a 6 as their first digit, with a few 7's and a few 4's. Execute the following sequence of cells to see an example where Benford's law does not hold.
MaleHeightsInAClassIOnceTaught = [5.5 ,
5.25, 5.62, 5.49, 5.52, 5.62, 5.41, 4.78, 6.26,
5.76, 5.99, 5.82, 6.02, 5.21, 5.55, 5.47, 5.57, 5.46,
5.4 , 5.36, 5.56, 5.81, 5.73, 5.19, 6.43, 6.24, 5.53,
5.47, 5.82, 5.5 , 5.72, 6.09, 5.74, 5.57, 5.77, 5.76,
5.85, 5.27, 5.23, 5.49, 5.71, 5.07, 6.04, 5.68, 5.67,
5.69, 5.86, 5.56, 6.1 , 5.84, 5.43, 5.67, 5.04, 5.91,
5.36, 5.99, 5.29, 6.94, 5.45, 5.24, 5.76, 5.85, 5.61,
5.64, 5.45, 5.07, 5.75, 5.9 , 5.85, 5.29, 5.69, 6.13,
5.54, 5.83, 5.29, 5.82, 5.42, 6.21, 5.74, 5.44, 5.95,
5.57, 5.51, 5.39, 5.8 , 5.55, 5.65, 5.71, 5.82, 5.56,
5.52, 5.59, 5.54, 5.56, 5.74, 5.7 , 5.73, 5.57, 5.79,
5.68, 5.66, 5.97, 5.95, 5.3 , 5.35, 5.57, 6.15, 5.23,
5.64, 6.01, 6.12, 5.46, 5.34, 5.46, 5.76, 5.37, 5.32,
5.47, 5.89, 5.52, 5.47, 5.51, 5.33, 5.51, 5.84, 5.67,
5.9 , 5.61 ]
NumbersAsStrings = [ str(number) for number in MaleHeightsInAClassIOnceTaught ]
NumbersAsStrings[:5] #first 5
PatternCounts = { '1':0, '2':0, '3':0, '4':0, '5':0,
'6':0, '7':0, '8':0, '9':0 }
for word in NumbersAsStrings:
for pattern in PatternCounts.keys():
if( re.match(pattern,word) ):
PatternCounts[pattern] += 1
PatternCounts
numdigits = [1,2,3,4,5,6,7,8,9]
strdigits = ['1', '2', '3', '4', '5', '6', '7', '8', '9']
TotalCount = len( NumbersAsStrings )
Frequencies = [ PatternCounts[digit]/TotalCount for digit in strdigits ]
plt.bar( numdigits, Frequencies, align='center')
plt.xticks( numdigits ) # make sure all the digits are on the x-axis
plt.xlabel("Digits")
plt.ylabel("Frequencies")
plt.title("Male heights do not seem to satisfy Benford's law.");
Gaussian distributed numbers ( those whose distributions resemble a "bell curve" ) in general do not tend to satisfy Benford's law, nor does uniformly randomly generated numbers such as those produced by a random number generator.
UniformlyRandom = numpy.random.randint(0,100,1000)
UniformlyRandom[:20]
NumbersAsStrings = [ str(number) for number in UniformlyRandom]
PatternCounts = { '1':0, '2':0, '3':0, '4':0, '5':0,
'6':0, '7':0, '8':0, '9':0 }
for word in NumbersAsStrings:
for pattern in PatternCounts.keys():
if( re.match(pattern,word) ):
PatternCounts[pattern] += 1
numdigits = [1,2,3,4,5,6,7,8,9]
strdigits = ['1', '2', '3', '4', '5', '6', '7', '8', '9']
TotalCount = len( NumbersAsStrings )
Frequencies = [ PatternCounts[digit]/TotalCount for digit in strdigits ]
plt.bar( numdigits, Frequencies, align='center')
plt.xticks( numdigits ) # make sure all the digits are on the x-axis
plt.xlabel("Digits")
plt.ylabel("Frequencies")
plt.title("Uniformly Random numbers do not satisfy Benford's law");
GaussianRandom = numpy.random.normal(50,25,1000)
GaussianRandom[:20]
NumbersAsStrings = [ str(number) for number in GaussianRandom]
PatternCounts = { '1':0, '2':0, '3':0, '4':0, '5':0,
'6':0, '7':0, '8':0, '9':0 }
for word in NumbersAsStrings:
for pattern in PatternCounts.keys():
if( re.match(pattern,word) ):
PatternCounts[pattern] += 1
numdigits = [1,2,3,4,5,6,7,8,9]
strdigits = ['1', '2', '3', '4', '5', '6', '7', '8', '9']
TotalCount = len( NumbersAsStrings )
Frequencies = [ PatternCounts[digit]/TotalCount for digit in strdigits ]
plt.bar( numdigits, Frequencies, align='center')
plt.xticks( numdigits ) # make sure all the digits are on the x-axis
plt.xlabel("Digits")
plt.ylabel("Frequencies")
plt.title("Gaussian (bell curve) Distributed numbers do not satisfy Benford's law");
Fraud Detection¶
In spite of its failures in some instances, Benford's law is a valuable tool -- and one that is necessarily implemented computationally. Population numbers, financial numbers, and many other sets of numbers can be expected to have a leading digit of 1 about 30% of the time.
To illustrate further, enter your username and then execute the cell to load additional data sets.
username = "hammese"
%run -i Assignment3Supplemental.py
If the cell executed correctly, there are now 7 additional data sets in memory. For example, CountyPopulations1990 is the population of each county in the United States in the 1990 census; and it satisfies Benford's law.
NumbersAsStrings = [ str(number) for number in CountyPopulations1990]
PatternCounts = { '1':0, '2':0, '3':0, '4':0, '5':0,
'6':0, '7':0, '8':0, '9':0 }
for word in NumbersAsStrings:
for pattern in PatternCounts.keys():
if( re.match(pattern,word) ):
PatternCounts[pattern] += 1
numdigits = [1,2,3,4,5,6,7,8,9]
strdigits = ['1', '2', '3', '4', '5', '6', '7', '8', '9']
TotalCount = len( NumbersAsStrings )
Frequencies = [ PatternCounts[digit]/TotalCount for digit in strdigits ]
plt.bar( numdigits, Frequencies, align='center')
plt.plot( range(1,10),BenProbs[1:], 'k-' )
plt.plot( range(1,10),BenProbs[1:], 'k.', markersize = 10 )
plt.xticks( numdigits ) # make sure all the digits are on the x-axis
plt.xlabel("Digits")
plt.ylabel("Frequencies")
plt.title("County Population 1990 distribution (blue) close to Benford's law (black). ");
If this was a statistics course, we would at this point introduce the Kolmogorov-Smirnov test for determining if two sets of data come from the sample distribution. However, we will simply notice visually that the distribution represented by the blue bars does resemble the theoretical (black lines and dots) Benford law distribution.
For example, the New York Stock Exchange closing prices on February 10, 2017 also have a Benford law distribution.
NumbersAsStrings = [ str(number) for number in NYSE10Feb2017]
PatternCounts = { '1':0, '2':0, '3':0, '4':0, '5':0,
'6':0, '7':0, '8':0, '9':0 }
for word in NumbersAsStrings:
for pattern in PatternCounts.keys():
if( re.match(pattern,word) ):
PatternCounts[pattern] += 1
numdigits = [1,2,3,4,5,6,7,8,9]
strdigits = ['1', '2', '3', '4', '5', '6', '7', '8', '9']
TotalCount = len( NumbersAsStrings )
Frequencies = [ PatternCounts[digit]/TotalCount for digit in strdigits ]
plt.bar( numdigits, Frequencies, align='center')
plt.plot( range(1,10),BenProbs[1:], 'k-' )
plt.plot( range(1,10),BenProbs[1:], 'k.', markersize = 10 )
plt.xticks( numdigits ) # make sure all the digits are on the x-axis
plt.xlabel("Digits")
plt.ylabel("Frequencies")
plt.title("NYSE closing prices distribution (blue) is close to Benford's law (black). ");
Financial data in general follow a Benford's law distribution -- unless it has been forged. That is,
_ __People tend to fabricate numbers that are not 'naturally random' and thus do not satisfy Benford's law__ _
It is, in fact, a standard practice for law enforcement and tax officials to use Benford's law to test if a suspect "cooked the books" -- i.e., made up numbers -- or if the numbers are naturally random as would be expected.
For example, here is financial data (sales figures) for a large corporation for its 2014 fiscal year.
NumbersAsStrings = [ str(number) for number in CorporateFinancial]
PatternCounts = { '1':0, '2':0, '3':0, '4':0, '5':0,
'6':0, '7':0, '8':0, '9':0 }
for word in NumbersAsStrings:
for pattern in PatternCounts.keys():
if( re.match(pattern,word) ):
PatternCounts[pattern] += 1
numdigits = [1,2,3,4,5,6,7,8,9]
strdigits = ['1', '2', '3', '4', '5', '6', '7', '8', '9']
TotalCount = len( NumbersAsStrings )
Frequencies = [ PatternCounts[digit]/TotalCount for digit in strdigits ]
plt.bar( numdigits, Frequencies, align='center')
plt.plot( range(1,10),BenProbs[1:], 'k-' )
plt.plot( range(1,10),BenProbs[1:], 'k.', markersize = 10 )
plt.xticks( numdigits ) # make sure all the digits are on the x-axis
plt.xlabel("Digits")
plt.ylabel("Frequencies")
plt.title("Typical Financial data lead-digit distribution (blue) resembles Benford's law (black). ");
But now consider fraudulent financial data. In particular, although many of the numbers are from actual transactions, a large fraction of these numbers were "made up" by a human in the act of committing financial fraud.
NumbersAsStrings = [ str(number) for number in FraudFinancial ]
PatternCounts = { '1':0, '2':0, '3':0, '4':0, '5':0,
'6':0, '7':0, '8':0, '9':0 }
for word in NumbersAsStrings:
for pattern in PatternCounts.keys():
if( re.match(pattern,word) ):
PatternCounts[pattern] += 1
numdigits = [1,2,3,4,5,6,7,8,9]
strdigits = ['1', '2', '3', '4', '5', '6', '7', '8', '9']
TotalCount = len( NumbersAsStrings )
Frequencies = [ PatternCounts[digit]/TotalCount for digit in strdigits ]
plt.bar( numdigits, Frequencies, align='center')
plt.plot( range(1,10),BenProbs[1:], 'k-' )
plt.plot( range(1,10),BenProbs[1:], 'k.', markersize = 10 )
plt.xticks( numdigits ) # make sure all the digits are on the x-axis
plt.xlabel("Digits")
plt.ylabel("Frequencies")
plt.title("Fraudulent data lead-digit distribution (blue) versus Benford's law (black). ");
This is, in fact, your task in the remainder of this notebook. Copy and paste code from above, and modify as needed, to answer the following:
Exercise 4¶
Do the first 100 numbers in the Fibonnacci sequence satisfy Benford's law? Explain why or why not.
First100FibonacciNumbers = [1,1,2,3,5,8,13,21,34,55,89,144,233,377,610,987,1597,2584,4181,6765,
10946,17711,28657,46368,75025,121393,196418,317811,514229,832040,
1346269,2178309,3524578,5702887,9227465,14930352,24157817,39088169,63245986,
102334155,165580141,267914296,433494437,701408733,1134903170,1836311903,
2971215073,4807526976,7778742049,12586269025,20365011074,32951280099,
53316291173,86267571272,139583862445,225851433717,365435296162,
591286729879,956722026041,1548008755920,2504730781961,4052739537881,
6557470319842,10610209857723,17167680177565,27777890035288,44945570212853,
72723460248141,117669030460994,190392490709135,308061521170129,
498454011879264,806515533049393,1304969544928657,2111485077978050,3416454622906707,
5527939700884757,8944394323791464,14472334024676221,23416728348467685,
37889062373143906,61305790721611591,99194853094755497,160500643816367088,
259695496911122585,420196140727489673,679891637638612258,1100087778366101931,
1779979416004714189,2880067194370816120,4660046610375530309,7540113804746346429,
12200160415121876738,19740274219868223167,31940434634990099905,
51680708854858323072,83621143489848422977,135301852344706746049,
218922995834555169026,354224848179261915075]
NumbersAsStrings = [ str(number) for number in First100FibonacciNumbers ]
PatternCounts = { '1':0, '2':0, '3':0, '4':0, '5':0,
'6':0, '7':0, '8':0, '9':0 }
for word in NumbersAsStrings:
for pattern in PatternCounts.keys():
if( re.match(pattern,word) ):
PatternCounts[pattern] += 1
numdigits = [1,2,3,4,5,6,7,8,9]
strdigits = ['1', '2', '3', '4', '5', '6', '7', '8', '9']
TotalCount = len( NumbersAsStrings )
Frequencies = [ PatternCounts[digit]/TotalCount for digit in strdigits ]
plt.bar( numdigits, Frequencies, align='center')
plt.plot( range(1,10),BenProbs[1:], 'k-' )
plt.plot( range(1,10),BenProbs[1:], 'k.', markersize = 10 )
plt.xticks( numdigits ) # make sure all the digits are on the x-axis
plt.xlabel("Digits")
plt.ylabel("Frequencies")
plt.title("First 100 Fibonacci Numbers (blue) versus Benford's law (black). ");
Based off of the graph shows that the First 100 Fibonacci Numbers do in face satisfy Benford's law. There are a few locations on the graph, such as digit 8, where there is a little discrepancy between the frequency of the Fibonacci and Benford's, but not enough for me to consider it to not satisfy.
Exercise 5¶
There are 3 "suspects" in a fraud investigation who submitted the following numbers in their financial documents -- KimmyTheKite, MannyTheMooch, and SammyTheSnake. Would you suspect any -- or all -- of committing fraud and why? (Show your work)
print( 'KimmyTheKite = {0}\n'.format(KimmyTheKite) )
print( 'MannyTheMooch = {0}\n'.format(MannyTheMooch) )
print( 'SammyTheSnake = {0}\n'.format(SammyTheSnake) )
NumbersAsStrings = [ str(number) for number in KimmyTheKite ]
PatternCounts = { '1':0, '2':0, '3':0, '4':0, '5':0,
'6':0, '7':0, '8':0, '9':0 }
for word in NumbersAsStrings:
for pattern in PatternCounts.keys():
if( re.match(pattern,word) ):
PatternCounts[pattern] += 1
numdigits = [1,2,3,4,5,6,7,8,9]
strdigits = ['1', '2', '3', '4', '5', '6', '7', '8', '9']
TotalCount = len( NumbersAsStrings )
Frequencies = [ PatternCounts[digit]/TotalCount for digit in strdigits ]
plt.bar( numdigits, Frequencies, align='center')
plt.plot( range(1,10),BenProbs[1:], 'k-' )
plt.plot( range(1,10),BenProbs[1:], 'k.', markersize = 10 )
plt.xticks( numdigits ) # make sure all the digits are on the x-axis
plt.xlabel("Digits")
plt.ylabel("Frequencies")
plt.title("KimmyTheKite (blue) versus Benford's law (black). ");
NumbersAsStrings = [ str(number) for number in MannyTheMooch ]
PatternCounts = { '1':0, '2':0, '3':0, '4':0, '5':0,
'6':0, '7':0, '8':0, '9':0 }
for word in NumbersAsStrings:
for pattern in PatternCounts.keys():
if( re.match(pattern,word) ):
PatternCounts[pattern] += 1
numdigits = [1,2,3,4,5,6,7,8,9]
strdigits = ['1', '2', '3', '4', '5', '6', '7', '8', '9']
TotalCount = len( NumbersAsStrings )
Frequencies = [ PatternCounts[digit]/TotalCount for digit in strdigits ]
plt.bar( numdigits, Frequencies, align='center')
plt.plot( range(1,10),BenProbs[1:], 'k-' )
plt.plot( range(1,10),BenProbs[1:], 'k.', markersize = 10 )
plt.xticks( numdigits ) # make sure all the digits are on the x-axis
plt.xlabel("Digits")
plt.ylabel("Frequencies")
plt.title("MannyTheMooch (blue) versus Benford's law (black). ");
NumbersAsStrings = [ str(number) for number in SammyTheSnake ]
PatternCounts = { '1':0, '2':0, '3':0, '4':0, '5':0,
'6':0, '7':0, '8':0, '9':0 }
for word in NumbersAsStrings:
for pattern in PatternCounts.keys():
if( re.match(pattern,word) ):
PatternCounts[pattern] += 1
numdigits = [1,2,3,4,5,6,7,8,9]
strdigits = ['1', '2', '3', '4', '5', '6', '7', '8', '9']
TotalCount = len( NumbersAsStrings )
Frequencies = [ PatternCounts[digit]/TotalCount for digit in strdigits ]
plt.bar( numdigits, Frequencies, align='center')
plt.plot( range(1,10),BenProbs[1:], 'k-' )
plt.plot( range(1,10),BenProbs[1:], 'k.', markersize = 10 )
plt.xticks( numdigits ) # make sure all the digits are on the x-axis
plt.xlabel("Digits")
plt.ylabel("Frequencies")
plt.title("SammyTheSnake (blue) versus Benford's law (black). ");
SammyTheSnake could easily be suspected for committing fraud based off of the results of the graph. KimmyTheKite could be a little iffy on whether or not she committed fraud... her graph is a little off when it comes to aligning with Benford's law, but not nearly as off as SammyTheSnake. MannyTheMooch seems to be safe. His pretty closely align with Benford's law.