Monday, June 10, 2013

Fitness function for multi objective optimization / Mapping vectors/tuples to numbers

When you want to let the computer discover how to best solve a problem, you can use a genetic algorithm to evolve a solution for a problem. For example, if you're trying to find which schedule of lessons will cause the least clashes (hopefully none), you create a population of possible schedules and let them evolve into better and better schedules. A more interesting example would be evolving a program which does something in particular, such as control a robot.

In genetic algorithms you have to supply a function called a fitness function which returns a number which quantifies how good a particular solution is at solving the problem. For example, in the case of the time table scheduler, a fitness function would be the number of clashes between lessons that are caused by a given schedule. The smaller this number, the better the schedule. If you're trying to evolve a program which generates prime numbers, then a fitness function would be the percentage of numbers generated which are actually prime numbers. The higher this number, the better the program.

The problems start when you want to maximize or minimize more than one thing, that is, when you have a number of separate fitnesses that you want to optimize together. An example of this is evolving your schedule such that both the number of lesson clashes (when students are expected to be in two different lessons at the same time) and the number of room clashes (when the same room is assigned to be used for more than one lesson at the same time) will be minimized. In the case of the program, you usually not only want to evolve a correct program, but also one which gives the correct output quickly. This is called multi objective optimization.

What is usually done is a weighted summation of the individual fitnesses. So you say that the correctness of a program is, say, twice as important as the speed of the program, so the fitnesses are combined as

fitness = 2 * percentage_correct_outputs + time_taken

Of course this is not correct since the correctness must go up whilst the time should go down. So a better combination might be

fitness = 2 * percentage_correct_outputs + 1/(time_taken + 1)

You might not think that this is the best fitness function but choosing a correct fitness function is part of the difficulty of using genetic algorithms so have fun finding a better one.

The problem with weighted summations is that the fitness will increase by simply increasing one of the sub fitnesses. So in the previous example, you can get a good fitness by simply creating a program that does nothing and takes 0 seconds to do so. This program is much easier to find than a program which is correct so the genetic algorithm tends to improve the execution time only and gets stuck as a program which increases the correctness will be much much slower and thus will reduce the fitness rather than increase it. A better way to combine the sub fitnesses is needed.

I found two ways to do this for two different needs. The first is for when the sub fitnesses are of equal importance, for example the time table scheduler minimizing the number of lesson and room clashes. The second is for when the sub fitnesses are prioritized such that the most important sub fitness should always be improved regardless of how the second is affected, for example the prime number generator being correct has the highest priority but between two equally correct programs, the faster one is preferred.

Equal priority fitness

If both sub fitnesses must be optimized together such that improving one without the other is not helpful, then we might add the two fitnesses together and subtract their absolute difference like this:

fitness = (sub1 + sub2) - |sub1 - sub2|

This way as one starts improving without the other, their difference will increase and the whole fitness will start becoming smaller. We can simplify this equation by calling the largest sub fitness "max" and the smallest sub fitness "min":

fitness = (max + min) - (max - min)
fitness = max + min - max + min
fitness = 2*min

So this fitness is equivalent to finding twice the minimum of the sub fitnesses, which makes sense since if you only care about increasing the minimum fitness, then both fitnesses will increase together. Of course multiplying by 2 is redundant since comparing a minimum with another minimum and comparing twice a minimum with twice another minimum will give the same result, so our final fitness function is:

fitness = min(sub1, sub2)

and you can add more sub fitnesses by just finding the minimum of all of them:

fitness = min(sub1, sub2, sub3, ..., subn)

So in our time table example, the fitness function would be:

fitness = min(-lesson_clashes, -room_clashes)

The negations are used so that in order to increase the fitness, the number of clashes have to be decreased.

Prioritized fitness

When one sub fitness is more important than the other, such as program correctness versus execution time, a different method is needed. What is needed is a function "f" which combines the two sub fitnesses ("sub_1" and "sub_2" where sub_1 has a greater priority than sub_2) such that the following conditions are met:

f(subA_1, subA_2) > f(subB_1, subB_2)
if subA_1 > subB_1 or
   subA_1 = subB_1 and subA_2 > subB_2

f(subA_1, subA_2) = f(subA_1, subB_2)
if subA_1 = subB_1 and subA_2 = subB_2

So for example f(2, 1) > f(1, 100), f(5, 20) > f(5, 1) and f(3, 3) = f(3, 3).

This is called a lexicographical ordering of the sub fitnesses. With strings this is natural. When sorting names, "John" comes before "Zach" because "J" comes before "Z", regardless of what the second letter is. But "James" comes before "John" because since the first letter is the same then we look at the second one and "a" comes before "o".

The problem we're facing is with finding a function which when given a pair of numbers will return a number which can be used to sort a list of pairs of numbers in lexicographical order. This is similar to Cantor's mapping of pairs of numbers to the natural numbers, except that this time we want a particular mapping which is lexicographic.

A natural lexicographic ordering in numbers comes when we look at numbers with a decimal point such as 2.1 and 3.5. The whole number part will always determine the ordering, unless they are equal, in which case the fractional part is then used. So if we can map our pairs of sub fitnesses to real numbers in such a way that the high priority sub fitness becomes the whole number part and the low priority sub fitness becomes the fractional part, then we would have found our "f".

What I found was that you can do this by squashing the low priority number into a proper fraction (a number between 0 and 1) and then adding it to the second number, assuming that it is a whole number. Unfortunately this method will only work if the high priority number is a whole number. The low priority number can be a real number however. In order to squash the low priority number, you can use a sigmoid function which given any number will return a number between 0 and 1 in such a way that ordering is preserved (sigmoid(a) > sigmoid(b) if and only if a > b). Another and perhaps simpler way would be to use the hyperbolic tangent or the arc tangent and then modify them so that their output is between 0 and 1 (y = (tanh(x)+1)/2 and y = (atan(x)+pi/2)/pi).

So the combined fitness would be:

fitness = sub1 + sigmoid(sub2)

The nice this about this is that you can add more sub fitnesses in the following way, assuming that only the least important sub fitness is a real number whilst the rest are integers:

fitness = sub1 + sigmoid(sub2 + sigmoid(sub3 + ... + sigmoid(subn)...))

So in our program example, the fitness function would be:

fitness = percentage_correct_outputs + sigmoid(time_taken)

It should be noted that this is a way to map vectors/tuples of integers to real numbers whilst preserving lexicographic ordering.

Combined

Now we can even use both of the combinations in order to combine sub fitnesses in complex ways where some are of equal priority whilst others are of different priority. For example, say we are evolving a time table schedule which minimizes lesson clashes, minimizes room clashes and minimizes density such that you avoid cramming all the lessons in one day. The lesson and room clashes are of equal priority but the density has a lower priority than the other two. So the fitness function might be:

fitness = min(-lesson_clashes, -room_clashes) - sigmoid(density)

This will give us a real number which has a whole number part representing the number of undesirable clashes and a fractional part representing the density, both of which must be minimized. Since fitness should be increased, minimization is done by using the negations and subtraction.

Sunday, May 5, 2013

An equation for enumerating fractions

In a previous post I described Cantor's way of proving that some infinite lists are bigger than others. I said that his is done by showing that a number of infinite lists, which are called countably infinite, can be paired up with every natural number from 0 onward. I shall now describe an equation that gives the integer that is paired up with a particular fraction.

The fractions are:

1/1 1/2 1/3 1/4 ...
2/1 2/2 2/2 2/3 ...
3/1 3/2 3/2 4/3 ...
4/1 4/2 4/2 3/3 ...
...

and can be paired up with the natural numbers by doing the following:

Organize them in a grid like above, where the first row consists of all the fractions with a numerator of 1, the second row with a numerator of 2, etc.
Start from the corner (1/1) and pair it with 0
Go down each diagonal in turn, that is, [1/2, 2/1] would be the first diagonal, [1/3, 2/2, 3/1] would be the second, etc. and pair them with the next consecutive natural number

This will give us the following mapping:
[0:1/1, 1:2/1, 2:1/2, 3:3/1, 4:2/2, ...]

In order to find an equation which does this mapping, we first need to know in which diagonal a given natural number would be paired up in. So 0 is paired up with the corner (the diagonal 0), 1 is paired up with 1/2 which is in the second diagonal (diagonal 1), 2 is with 2/1 which is in the second (diagonal 1), 3 with 3/1 in the third (diagonal 2), etc. This will let us know what the nth fraction in the diagonal is.

Let's start by determining which natural number would be paired up with the first fraction in each diagonal. First fraction in diagonal 0 is paired up with 0, the first fraction in diagonal 1 is paired up with 1, in diagonal 2 with 3, in diagonal 3 with 6, etc. Notice that, since the size of each diagonal increases by 1 each time, the sequence of numbers just mentioned (0, 1, 3, 6, ...) is an arithmetic series. This means that the first fraction in diagonal 0 is paired up with 0, diagonal 1 with 0+1, diagonal 2 with 0+1+2, diagonal 3 with 0+1+2+3, etc. So in order to find which natural number n will be paired up with the first fraction in diagonal d, we use:

n = d(d+1)/2

d = 0 => n = 0
d = 1 => n = 1
d = 2 => n = 3
d = 3 => n = 6
d = 4 => n = 10

We can use this equation to find the diagonal d in which natural number n is paired up. This is done by making d subject of the formula using the completing the square method: d = -0.5 ± sqrt(2n + 1/4)

This new equation gives the following:

d = 0 => n = 0, -1
d = 1 => n = 1, -2
d = 2 => n = 1.56, -2.56
d = 3 => n = 2, -3
d = 4 => n = 2.37, -3.37
d = 5 => n = 2.7, -3.7
d = 6 => n = 3, -4
d = 7 => n = 3.27, -4.27
d = 8 => n = 3.53, -4.53
d = 9 => n = 3.77, -4.77
d = 10 => n = 4, -5

As you can see, the answer that comes out when using the minus sign is non-nonsensical, whilst the whole number part of the answer that comes out when using the plus sign is the diagonal position in which natural number y belongs. The plus part is what we want and we can extract the whole number part by using the floor function:

d = floor(-0.5 + sqrt(2n + 1/4))
or
d = floor((sqrt(8n + 1) - 1)/2)

This gives us:

n = 0 => d = 0
n = 1 => d = 1
n = 2 => d = 1
n = 3 => d = 2
n = 4 => d = 2
n = 5 => d = 2
n = 6 => d = 3
n = 7 => d = 3
n = 8 => d = 3
n = 9 => d = 3
n = 10 => d = 4

So natural number 0 will be paired up with a fraction in diagonal 0, 1 will be paired up with a fraction in diagonal 1, 2 in diagonal 1, 3 in diagonal 2, 4 in diagonal 2, etc.

Since the largest numerator and denominator in a given diagonal is 1 more than the position of the diagonal (diagonal 2 has a largest numerator, as well as a largest denominator, of 3, etc.), we now also know that the largest numerator and denominator in the diagonal in which the natural number n belongs is floor((sqrt(8n + 1) - 1)/2)+1. So the first fraction in that diagonal will be 1/(floor((sqrt(8n + 1) - 1)/2)+1) whilst the last fraction will be (floor((sqrt(8n + 1) - 1)/2)+1)/1.

Since we also know which natural number is paired up with the first fraction in a given diagonal, we can know which natural number is paired up with the first fraction in the diagonal in which n is paired.

The natural number n which is paired up with the first fraction in diagonal d is: n = d(d+1)/2
The diagonal d in which natural number n is paired up is: d = floor((sqrt(8n + 1) - 1)/2)
So the natural number n which is paired up with the first fraction in the diagonal which natural number m is paired up is: n = floor((sqrt(8m + 1) - 1)/2)(floor((sqrt(8m + 1) - 1)/2) + 1)/2

This gives us:

m = 0 => n = 0
m = 1 => n = 1
m = 2 => n = 1
m = 3 => n = 3
m = 4 => n = 3
m = 5 => n = 3
m = 6 => n = 6
m = 7 => n = 6
m = 8 => n = 6
m = 9 => n = 6
m = 10 => n = 10

Observe that in a given diagonal such as [1/4, 2/3, 3/2, 4/1], you can find what any of the other fractions will be by knowing which is the first fraction and how far the fraction you want to find is from it (the numerator increases whilst the denominator decreases). For example, the 0th fraction in the diagonal is (1+0)/(4-0), the 1th fraction is (1+1)/(4-1), the 2th fraction is (1+2)/(4-2), etc.

We know how far m is from the first natural number n in the diagonal in which m resides:
m - floor((sqrt(8m + 1) - 1)/2)(floor((sqrt(8m + 1) - 1)/2) + 1)/2

And we know the first fraction in the diagonal in which m resides:
1/(floor((sqrt(8m + 1) - 1)/2)+1)

So the fraction which will be paired up with natural number n is:

numerator: 1 + (n - floor((sqrt(8n + 1) - 1)/2)(floor((sqrt(8n + 1) - 1)/2) + 1)/2)
denominator: floor((sqrt(8m + 1) - 1)/2)+1 - (n - floor((sqrt(8n + 1) - 1)/2)(floor((sqrt(8n + 1) - 1)/2) + 1)/2)

which I'm too lazy to simplify.

The equation in Python 3 is as follows:

for n in range(20):
    print(n, ":", int(1 + (n - math.floor((math.sqrt(8*n + 1) - 1)/2)*(math.floor((math.sqrt(8*n + 1) - 1)/2) + 1)/2)), "/", int(math.floor((math.sqrt(8*n + 1) - 1)/2)+1 - (n - math.floor((math.sqrt(8*n + 1) - 1)/2)*(math.floor((math.sqrt(8*n + 1) - 1)/2) + 1)/2)))

And here is the output:

So using this enumeration method, the fraction 100/173 will be paired with:

def findNat(num, den):
    while int(1 + (n - math.floor((math.sqrt(8*n + 1) - 1)/2)*(math.floor((math.sqrt(8*n + 1) - 1)/2) + 1)/2)) != num and int(math.floor((math.sqrt(8*n + 1) - 1)/2)+1 - (n - math.floor((math.sqrt(8*n + 1) - 1)/2)*(math.floor((math.sqrt(8*n + 1) - 1)/2) + 1)/2)) != den:
        n += 1
    return n
print(findNat(100, 173))

Friday, April 12, 2013

The Monty Hall problem

This is an interesting mathematical problem in probability which has an intuitive answer. The only problem is that the intuitive answer is wrong, which is what I love about it. It's called the Monty Hall problem and it goes like this:

Imagine you're in a game show and presented with three doors.

Behind one of the doors is a prize. You are to guess behind which door it is. At this point we should all agree that the chance of guessing the right door would be 1/3. You decide to pick door number 1.

You think that you have sealed you fate, but in an surprising turn of events, the game show host says that in order to help you, since at least one of the two doors left must be wrong, the host is going to tell you that door number 3 is surely wrong.

So now you know that the prize is either behind the door you chose or the other remaining door, door number 2. The host now asks you if you'd rather stick with your first choice or choose the other door.

Did the host actually help you? Is there some advantage to knowing which one of the three doors is surely wrong? Intuitively you'd say that your chances of winning by sticking with your first choice or choosing door number 2 are both 1/2. It's either behind one or the other right? So your chances of winning are equal right? Right guys? Guys?

Most people will not be convinced that this intuition is false. So here's a program I've written in Python 3 which simulates this situation 100 times and it will count the number of times you would have won if you stayed with the first choice and the number of time you would have won if you swapped your choice.

import random

stay_win = 0
swap_win = 0
doors = {"Door 1", "Door 2", "Door 3"}
for _ in range(100): #For 100 times do the following...
    #Select the correct door
    correct = random.choice(list(doors))
    #The contestant makes their choice
    first_choice = random.choice(list(doors))

    #The door which can be eliminated cannot be the correct door or the contestant's choice
    can_eliminate = doors - { correct, first_choice }
    #Select the eliminated door
    eliminated = random.choice(list(can_eliminate))

    #The contestant has only one door left to choose if they swap their door
    choice_left = doors - { first_choice, eliminated }
    #Make this unchosen, non-eliminated door the second choice for swapping the door
    second_choice = choice_left.pop()

    if first_choice == correct:
        stay_win += 1
    if second_choice == correct:
        swap_win += 1

    print("Correct:", correct, " | ", "First choice:", first_choice, " | ", "Eliminated:", eliminated, " | ", "Second choice:", second_choice, " | ", "Won if stay:", first_choice == correct)
print()
print("stay:", stay_win, "swap:", swap_win)

In order to let you scrutinize the result, all the individual situations will be outputted as well:

stay: 29 swap: 71

What? That can't be right! You win by swapping almost twice as often as when you stay with your first choice. Let's try it again!

stay: 39 swap: 61

Um...

stay: 37 swap: 63

Hmm... Why would you win twice as often when you swap than when you stay? The only way to know this is by seeing what needs to happen in order to win by sticking with your first choice and to win by swapping.

Here is what the different possibilities are when the correct door is door number 1:

And here is what the different possibilities are when the correct door is door number 2:

And finally here is what the different possibilities are when the correct door is door number 3:

Did you get it? No? Well, look carefully at what needs to happen in order to win using the first choice. You need to guess the correct door on first try. What are the chances of you doing that? It's 1/3 of course. Could you ever win by choosing a wrong door on first try and stick to it? No.

Well then what do you need to do in order to win by swapping doors? Clearly, you will only lose if you choose the correct door on first try, but if it's the wrong door you've chosen on first try then you'll win. What are the chances of you choosing the wrong door on first try? 2/3, twice as many as choosing the correct one.

Clearly you're twice as likely to win by swapping than by staying. Obviously this doesn't mean that you'll always win, it just means that you're more likely to win. This illustrates a very important point which is that intuition isn't perfect. Sometimes it's just not enough to use intuition in order to predict what will happen, you will need to get down and dirty and actually do boring work in order to make a correct judgement. Never overestimate yourself.

Friday, March 15, 2013

Fractions that approximate pi

We all know that 22/7 is an approximation of pi, although some think that it is exactly pi. It is not, it is just close. But why do we use 22/7 when there are other fractions that approximate pi better? Let's see what some of these fractions are.

To find these fractions I've written a Python 3 program which calculates the errors of various fractions when compared to pi. It works in the following way:

First, it considers a range of denominators to use for fractions. For each denominator, the numerator which makes the fraction closest to pi is found (for example, with a denominator of 7, the closest fraction would be 22/7).

Each fraction will get a number of digits at the beginning the same as pi. We would like to find which is the closest fraction to pi which gets, say, the first 5 digits exactly like pi's. So the best fraction for every number of exact digits is found.

Approximations are only useful if they are easy to use and remember and so we will be considering different lengths of numbers when performing the above process. First we shall find the best fractions which consist of a denominator which is less than 10, then which is less than 100, then 1000, etc. in order to control the sizes of the fractions.

Here is the code:

import math
import collections
 
def closestFraction(denominator, target):
    #The closest fraction which uses a particular denominator d to a target number t requires that we find the best numerator n such that n/d = t, so n = td, the closest whole number of which is round(td)
    numerator = round(denominator*target)
    fraction = numerator/denominator
    return (numerator, denominator, fraction)
 
def closestFractions(maxDenominatorRange, target):
    fractions = []
    fractionsFound = set()
    for denominator in range(1, maxDenominatorRange):
        (numerator, denominator, fraction) = closestFraction(denominator, target)
        if fraction not in fractionsFound: #avoid different forms of the same fraction due to non-simplification
            fractionsFound.add(fraction)
            fractions.append((numerator, denominator, fraction))
    return fractions
 
print("pi =", math.pi)
print()
for maxDenominatorRange in [10, 100, 1000, 10000, 100000]:
    #Group fractions by the number of first digits which are exactly like pi's
    numFirstDigitsDict = collections.defaultdict(list)
    for (numerator, denominator, fraction) in closestFractions(maxDenominatorRange, math.pi):
        numFirstDigits = 0
        while round(fraction, numFirstDigits) == round(math.pi, numFirstDigits):
            numFirstDigits += 1
        numFirstDigits -= 1
        numFirstDigitsDict[numFirstDigits].append((numerator, denominator, fraction))
 
    #Keep only the best fraction in each first digits group
    bestNumFirstDigitsDict = dict()
    for numFirstDigits in numFirstDigitsDict:
        bestNumFirstDigitsDict[numFirstDigits] = min(numFirstDigitsDict[numFirstDigits], key=lambda x:abs(x[2] - math.pi))
 
    print("==========")
    print("Fractions with denominator less than", maxDenominatorRange)
    for numFirstDigits in bestNumFirstDigitsDict:
        (numerator, denominator, fraction) = bestNumFirstDigitsDict[numFirstDigits]
        print("Best fraction with", numFirstDigits, "starting digits like pi's (rounded):", numerator, "/", denominator, "=", fraction)
    print()

And here are the results:

pi = 3.141592653589793

==========
Fractions with denominator less than 10
Best fraction with 0 starting digits like pi's (rounded): 19 / 6 = 3.1666666666666665
Best fraction with 1 starting digits like pi's (rounded): 25 / 8 = 3.125
Best fraction with 2 starting digits like pi's (rounded): 22 / 7 = 3.142857142857143

==========
Fractions with denominator less than 100
Best fraction with 0 starting digits like pi's (rounded): 167 / 53 = 3.150943396226415
Best fraction with 1 starting digits like pi's (rounded): 195 / 62 = 3.1451612903225805
Best fraction with 2 starting digits like pi's (rounded): 311 / 99 = 3.1414141414141414

==========
Fractions with denominator less than 1000
Best fraction with 0 starting digits like pi's (rounded): 167 / 53 = 3.150943396226415
Best fraction with 1 starting digits like pi's (rounded): 412 / 131 = 3.145038167938931
Best fraction with 2 starting digits like pi's (rounded): 2975 / 947 = 3.1414994720168954
Best fraction with 3 starting digits like pi's (rounded): 3085 / 982 = 3.1415478615071284
Best fraction with 4 starting digits like pi's (rounded): 2818 / 897 = 3.141583054626533
Best fraction with 6 starting digits like pi's (rounded): 355 / 113 = 3.1415929203539825

==========
Fractions with denominator less than 10000
Best fraction with 0 starting digits like pi's (rounded): 167 / 53 = 3.150943396226415
Best fraction with 1 starting digits like pi's (rounded): 412 / 131 = 3.145038167938931
Best fraction with 2 starting digits like pi's (rounded): 15541 / 4947 = 3.1414998989286436
Best fraction with 3 starting digits like pi's (rounded): 28497 / 9071 = 3.1415499944879284
Best fraction with 4 starting digits like pi's (rounded): 26669 / 8489 = 3.1415950053009776
Best fraction with 5 starting digits like pi's (rounded): 31218 / 9937 = 3.1415920297876623
Best fraction with 6 starting digits like pi's (rounded): 355 / 113 = 3.1415929203539825

==========
Fractions with denominator less than 100000
Best fraction with 0 starting digits like pi's (rounded): 167 / 53 = 3.150943396226415
Best fraction with 1 starting digits like pi's (rounded): 412 / 131 = 3.145038167938931
Best fraction with 2 starting digits like pi's (rounded): 15541 / 4947 = 3.1414998989286436
Best fraction with 3 starting digits like pi's (rounded): 28497 / 9071 = 3.1415499944879284
Best fraction with 4 starting digits like pi's (rounded): 294069 / 93605 = 3.14159500026708
Best fraction with 5 starting digits like pi's (rounded): 198379 / 63146 = 3.1415924999208182
Best fraction with 6 starting digits like pi's (rounded): 308429 / 98176 = 3.141592649934811
Best fraction with 7 starting digits like pi's (rounded): 209761 / 66769 = 3.141592655274154
Best fraction with 8 starting digits like pi's (rounded): 208341 / 66317 = 3.1415926534674368
Best fraction with 9 starting digits like pi's (rounded): 104348 / 33215 = 3.141592653921421
Best fraction with 10 starting digits like pi's (rounded): 312689 / 99532 = 3.1415926536189365

In my opinion, the best fraction is 355/113 which is easy to remember and is accurate up to 6 digits:

pi      = 3.141592653589793
355/113 = 3.1415929203539825

Tuesday, March 12, 2013

Montecarlo method of finding the area of a circle and pi

You may know that the digits of pi look random but did you know that you can use randomness to find pi? You can actually approximate the area of any circle using random points in a method called the Montecarlo method. The principle of this method uses something called "sampling" which works as follows:

Say you want to find the number of people in a country who are female. You don't have time to check every single person so instead you take a random sample of people, say 10% of the population, and check only those. You find that 40% of the people in this sample are female so you generalize this percentage to the whole population and say that approximately 40% of the country is female. If the sample was random, this would be a reasonable approximation and this sort of thing is done all the time in surveys.

Now what we want to do is to find the area of a circle such as this:

The Montecarlo method works in the following way:

Place that circle inside a square such as this:

Keep in mind that the area of a square is easy to find. Now throw a point inside the square at random, which would result in something like this:

Now if we repeat this process with many points and count how many of the points lie inside the circle, we can approximate what fraction of the area of the square is taken by the circle by treating the points as a sample of the total infinite number of points in the square.

For example, in the above diagram there are 17 points in total but only 14 of the points lie inside the circle. So we can say that approximately 14/17 of the area inside the square is taken by the circle.

Now lets assume that we have a decent number of points and have a good approximation of the fraction of area that is taken by the circle. How can we use this fraction to find the area of the circle? Just multiply the fraction by the area in the square. Since we know the fraction of the area of the square that lies inside the circle, we can find the area of the circle by finding what this fraction of the area actually is.

So in the above example, the length of the side of the square is about 6cm, so the area of the square is 6x6 = 36cm^2, so the area of the circle is approximately 14/17 x 36 = 29.64cm^2.

In reality, since the circle has a radius which is half as big as the side of the square, its radius is 3cm, so its area is pi x 3^2 = 28.27cm^2. Not bad for a sample of 17 points which weren't exactly placed randomly.

How can this be used to find an approximation of pi? Just approximate the area of a circle of radius 1. The area would then be pi x 1^2 = pi.

Now let's use a computer to create lots of points. How can we check if a point is inside the circle using a computer? We can simply check if the distance of the point to the center of the circle is less than or equal to the radius of the circle. If the distance to the center is greater than the radius, than the point cannot be inside the circle.

We find the distance of a point to the center by using Pythagoras' theorem as follows:

Here is the Python 3 program we'll be using:

import math, random

def isPointInCircle(x, y, Cx, Cy, radius):
    return math.sqrt((x - Cx)**2 + (y - Cy)**2) <= radius

def approximateCircleArea(radius, numberOfPoints):
    squareSide = radius*2
    Cx = radius
    Cy = radius

    pointsInside = 0
    for i in range(numberOfPoints):
        x = random.random()*squareSide
        y = random.random()*squareSide

        if (isPointInCircle(x, y, Cx, Cy, radius)):
            pointsInside = pointsInside + 1

    return pointsInside / numberOfPoints * squareSide**2

and we'll be finding what the area of a circle of radius 1 is for different numbers of points. Since the radius is 1, the area should approximate pi:

Number of points	Approximate area (should be 3.14159...)
10	2.4
100	3.0
1000	3.152
10000	3.14
10000	3.14896
100000	3.14112

Thursday, February 7, 2013

Maximum Likelihood Estimate / Laplace's law

What is the maximum likelihood estimate in statistics? Let's use an example from natural language processing.

Let's say that we want to estimate the probability/likelihood of every word in the English language vocabulary. What this means is that you want to estimate the probability that a randomly chosen word happens to be a particular word. There are many uses for such an estimate, called a unigram, such as attempting to predict what the next word to be written by a user will be (in order to type it for the user). Morse code gives the most frequent letters (letters with a highest probability of being used) the shortest codes in order to make sending messages faster. Something similar can be used for words in order to speed up the searching of words and to compress text files.

How can we estimate such a probability? What is usually done is that a large amount of text is collected, called a corpus, and the words in the text are counted. By dividing the number of times each word occurs by the total number of words in the corpus we will have an estimate of the probability for each word. This assumes there is a large number of words and that the words are used naturally (not a single word repeated over and over for example).

So P(w) = F(w)/F(*),
where P(w) is the probability of using a word called "w", F(w) is the frequency or number of times that the word "w" is found in the corpus and F(*) is the total number of words in the corpus.

What is the problem with this estimation? The problem is that any words which are not found in the corpus are assumed to have a probability of zero, that is, that they don't exist. This is not realistic, as even a very large corpus would not have every possible word ever, and new words are invented all the time. What we need to do is give a little less probability to the words we have in the corpus in order to share it with the words which are not in the corpus. This is called smoothing.

As it is, the way we are estimating the probability gives all the probability to the words in the corpus, which is the maximum probability you can give them. For this reason, this method is called the maximum likelihood estimate.

Now how can we smooth out our probability to give the unknown words a little of the probability too? A simple way we can do this is by assuming that each unknown word has an equal probability which is less than the probabilities of the words in the corpus. We can do this easily by assuming that each unknown word has a frequency of 1 and that each known word has a frequency of one more than it actually is. By increasing the number of times that each word occurs in the corpus by 1, including the words which aren't there, our probabilities will change such that the unknown words will all have an equal probability which is less than any known word's probability. In effect, no word will have a probability of zero, since every word would now look like it occurs at least once.

We might be tempted to say that the new probability is "(F(w) + 1)/F(*)", but this would be wrong because when you add together the probabilities of each word the total should come to 1, which it will not using this equation. As it was before, if you add together all the probabilities in the form of "F(w)/F(*)", the total would come to 1. For example if our corpus consisted of just one sentence "the boy played with the girl", then the probabilities of each word would be 2/6 for "the", 1/6 for "boy", 1/6 for "played", 1/6 for "with" and 1/6 for "girl". 2/6 + 1/6 + 1/6 + 1/6 + 1/6 = 1. If we use the wrong formula on the other hand, the probabilities would be 3/6 for "the", 2/6 for "boy", 2/6 for "played", 2/6 for "with" and 2/6 for "girl". 3/6 + 2/6 + 2/6 + 2/6 + 2/6 = 11/6 ≈ 1.8.

What we need to do is to add 1 to the denominator of the probability fraction for every different word we can find the probability of. What this means is that we need to know how many unknown words together with the known words there can be. This is called the vocabulary size and assuming that you know it you just need to add it to the denominator of the probability fraction like this:

P(w) = (F(w) + 1)/(F(*) + V),
where V is the vocabulary size or the total number of words including unknown ones.

This method is called Laplace's law, or rule of succession or add-one smoothing. It's not a very realistic estimation either however, since unknown words will not have equal probabilities and will not necessarily have a probability less than the rarest words in the corpus. Also, giving the unknown words a frequency of 1 might look like it's small, but when you consider that it is just one less than the smallest frequency in the corpus then it might be too large an estimate, especially considering that we might be including words in our vocabulary size which just don't exist (maybe no one uses them anymore).

There are other smoothing methods which can be used, such as using synonyms of unknown words. By using the probability of synonyms which are found in the corpus you can find an estimation of the unknown word's probability.

Let's see another example of how maximum likelihood estimate and Laplace's law can be used. Let's say that we want to do something more interesting than finding the probabilities of unigrams (single words) and instead want to find the probabilities of bigrams, that is, pairs of words. The nice thing about bigrams is that you can use a word in order to predict what the next word might be in a sentence. This would be useful for making writing text messages faster by having your cell phone predict and write for you the next word you want to write. Using more words than just the previous word in order to predict the next word would make for more accurate guesses, but let's just consider bigrams instead of the full "ngrams".

The probability we want to find is that of the user writing a given word provided that the previous word he/she wrote is known. So for example, we want to know what the probability of the user writing "cream" is after he/she has written the word "ice". We can then do this for every word and sort them by their probability in order to present the user with a list of suggested next words. Using maximum likelihood estimate, the probability is found by taking a corpus and counting the number of times the word "ice cream" is found and dividing that with the number of times the word "ice" is found. That way we'll estimate the probability of the word "ice" being followed by the word "cream".

P(w1 w2) = F(w1 w2)/F(w1),
where w1 is the first word and w2 is the next word.

Using Laplace's law we need to add 1 to the numerator and add the total number of possible word pairs you can have (known and unknown) to the denominator. The total number of possible word pairs might be calculated simply as the vocabulary size squared, that is, the vocabulary size multiplied by itself. This assumes that every possible word can follow every possible other word, include the same word. This isn't realistic at all but it might do.

P(w1 w2) = (F(w1 w2) + 1)/(F(w1) + V^2)

Saturday, January 12, 2013

Summary of research paper "Discovering Word Senses from Text" by Patrick Pantel and Dekang Lin

This is a summary of the research paper http://webdocs.cs.ualberta.ca/~lindek/papers/kdd02.pdf. This paper describes a way to group together similar words according to the different meanings of each word and all by just using the contexts of the words. You should also check the PhD thesis of the algorithm used in this paper by Panten at http://www.patrickpantel.com/download/papers/2003/cbc.pdf. You might also be interested in seeing this Java implementation of the algorithm at https://github.com/fozziethebeat/S-Space/blob/master/src/main/java/edu/ucla/sspace/clustering/ClusteringByCommittee.java.

In order to judge the similarity between two words, the grammatical relations between those words and other words as used in a corpus were used. For example, US state names can be recognized as similar because they all share the following contexts:
___ appellate court
campaign in ___
___ capital
governor of ___
___ driver's license
illegal in ___
primary in ___
___'s sales tax
senator for ___

So, for example, by counting the number of times the verb "campaign" is grammatically related to each word in a corpus, that count can be used to find other words which are similar because they have similar counts. The grammatically related word is called a feature of the word. By finding all the different features of each word in a corpus, we can give each word a list of its features and the number of times that feature was encountered as a context of the word. For example, we can find the number of times the word "campaign" appears as a feature of "California", as well as the word "capital", "governor", etc. All these features when grouped together are called a feature vector of the word.

However, in this paper, it is not the counts, or frequencies, that are associated with the features. Using raw counts to judge the importance of a feature is misleading as the feature might just be a common feature in general and not useful to find similar words by checking if they have features in common. So instead the pointwise mutual information is used. This compares the probability of the two word co-occurring together with the probability that they do so by coincidence (independently). In this paper it is calculated as:

mi(w,c) = p(w,c) / (p(w) × p(c))

where "w" is the word, "c" is a context of "w", "p(w,c)" is the probability of finding "c" as a context of "w", "p(w)" is the probability of finding the word "w" in some context and "p(c)" is the probability of finding "c" as a context of some word. These probabilities are calculated as follows:

p(w,c) = f(w,c) / N
p(w) = sum(f(w,c) for c in contexts) / N
p(c) = sum(f(w,c) for w in words) / N
N = sum(sum(f(w,c) for w in words) for c in contexts)

where "f(w,c)" is the number of times "c" is a context of "w" and "N" is the total number of word-context combinations.

However the pointwise mutual information gives bigger values to rare words and contexts so it is made fair by being multiplied by a discounting factor:

df(w,c) = f(w,c) / (f(w,c)+1) × min(f(w), f(c)) / min(f(w), f(c))

where

f(w) = sum(f(w,c) for c in contexts)
f(c) = sum(f(w,c) for w in words)

In order to measure the similarity between two words, the cosine coefficient is used:

sim(a,b) = sum(mi(a,c) × mi(b,c) for c in contexts) / sqrt(sum(mi(a,c)^2 for c in contexts) × sum(mi(b,c)^2 for c in contexts))

Words with different meanings, or senses, are called polysemous words. In order to cluster the words together according to their senses, using a different cluster for each sense, three phases are used.

In the first phase, the top 10 similar words for each word were found. In order to speed up the search, the feature vectors were indexed by their features, such that you can find all the feature vectors that have a given feature. The index is then used to find feature vectors that share features in common with the word we are finding similar words to. In this way you skip all the words which share nothing in common with the word of interest. In order to further speed things up, only some of the word's features are considered. The features which have a small value (pointwise mutual information) are ignored when looking for similar words. So only words which share these high valued features are retrieved from the index. So if we are looking for similar words to a word with the following feature vector:

{
  ("object of", "carry") : 5.118,
  ("object of", "bring") : 3.364,
  ("object of", "look at") : 2.753,
  ("object of", "have") : 1.353
}

and we only consider features with a value higher than 3.0 then we will use the index to find all the words which have a feature vector with the features

{
  ("object of", "carry") : 5.118,
  ("object of", "bring") : 3.364
}

This narrows down the number of words to check the similarity of (using the sim function) and so speeds up the search greatly. The top 10 most similar words are then retrieved. The fact that words which don't share high valued features are ignored did not seem to affect the result.

In the second phase, a clustering algorithm called "clustering by committee" is applied to the words using the similarity lists of the first phase. The clustering algorithm tries to assign the words into small clusters called "committees". Words which don't fit into any committee are called "residue". The point of these small clusters is to obtain centroids from which to cluster all the words, just like the randomly chosen centroids for k-means clustering, but instead the centroids are chosen intelligently. These centroids would originate from committees of words all of which mean a particular word sense. In this way the committee would decide whether to allow a word in the cluster or not. The way each committee is chosen is such that it is composed of words which do not have multiple meanings as used in the corpus. This makes the centroid of its feature vectors tend to only have features which belong to a particular word sense rather than to multiple senses, which leads to purer clusters of word senses. For example, when clustering USA state names, state names which are also names of cities (such as New York) are not included in the committee and only state names which mean nothing else are included which would avoid clustering city names in with the state name cluster when full clustering is performed. This is the clustering algorithm:

def clustering_by_committee(words, simlists, committeeSimilarityThreshold, residueSimilarityThreshold):
    highest_scoring_clusters = []
    #STEP 1.1
    for word in words:
        #STEP 1.2
        clusters = averagelink_cluster(simlists[word]) #use average link hierarchical agglomerative clustering to cluster the words into similarity clusters such that clusters can only be merged if their similarity is below a minimum similarity threshold
        #STEP 1.3
        for cluster in clusters:
            scores[cluster] = size(cluster) * average_pairwise_similarity(cluster)
        #STEP 1.4
        highest_scoring_clusters.append(arg_max(scores)) #append the highest scoring cluster

    #STEP 2.1
    sort_descending(highest_scoring_clusters) #sort highest_scoring_clusters in descending order of score

    #STEP 3.1
    committees = []
    #STEP 3.2
    for cluster in highest_scoring_clusters:
        #STEP 3.3
        centroid = compute_centroid(cluster) #the centroid is computed by averaging the frequencies (not the pointwise mutual information) of the features of each word in the cluster. these average frequencies are used to compute the pointwise mutual information which will be used to create the centroid vector (just as the feature vectors were constructed).
        #STEP 3.4
        if similarity(centroid, any committee in committees) < committeeSimilarityThreshold: #a committee cannot be very similar to other committees
            committees.append(cluster)

    #STEP 4.1
    if committees == []:
        return committees

    residues = []
    #STEP 5.1
    for word in words:
        #STEP 5.2
        if similarity(word, any committee in committees) < residueSimilarityThreshold: #if a word is not sufficiently close to any committee then it is a residue and more committees are needed for this word to belong to. otherwise the word is "covered" by a committee.
            residues.append(word)

    #STEP 6.1
    if residues == []:
        return committees
    #STEP 6.2
    else:
        return committes + clustering_by_committee(residues, simlists, threshold1, threshold2)

In the third phase, the committees discovered in the second phase are used to actually cluster words by their different senses. The words are not added to committees themselves, rather to clusters each of which is represented by its own committee. The similarity of a word to a cluster is the similarity between the word's feature vector and the centroid of the cluster's committee (which remains the same regardless of which words are added to the cluster). If hard clustering is used then words are only added to the most similar cluster, but if soft clustering is used, as is the case here, then the same word can be added to more than one cluster, one for each sense of the word. Here is the algorithm which assigns words to clusters:

def assign_word_to_clusters(word, committeeClusters, minWordSimilarity, maxClusterSimilarity):
    clusters = []
    similar_clusters = top_200_most_similar(committeeClusters, word) #find the top 200 clusters whose committee's centroid is most similar to the word
    while similar_clusters != []:
        best_cluster = get_most_similar_cluster(similar_clusters, word)
        if similarity(word, best_cluster) < minWordSimilarity:
            break #there is no cluster left that is similar enough to assign the word to
        else if similarity(best_cluster, any cluster in similar_clusters) < maxClusterSimilarity: #a word cannot be assigned to a cluster whose committee's centroid is similar to another cluster's committee's centroid
            best_cluster.add(word)
            remove_common_features(word, best_cluster) #remove the features from the word's feature vector which are also found in best_cluster's committee's centroid
        similar_clusters.remove(best_cluster)

The fact that the features in the word which are also found in the cluster's committee's centroid are removed before it is assigned to another cluster is important to allow it to be assigned to other less common sense clusters. For example, the word "plant" would first be assigned to the "factory" sense cluster in a newspaper corpus because that is the meaning of the word most of the time in such a corpus. After it is assigned to the factory sense cluster, the features which have to do with this sense are removed (or their values are set to 0) from the feature vector of "plant". Without these features, the remaining features would be of a different sense such as of living plant and so now it can be easily assigned to the living plant sense as these features are not overshadowed anymore. So this system will output a list of clusters to which each sense of every word will belong. In order to evaluate it, the sense clusters are compared to WordNet synonym groups called synsets. This is how the comparison is made: First, the similarity measure proposed by Lin which measures the similarity between two WordNet synsets is used. This is described in another post here.

sim_synset(s1, s2) = 2 log P(most_specific_common_synset(s1, s2)) / (log P(s1) + log P(s2))

The similarity between a word and a synset is the maximum similarity between the synset and one of the synsets in which the word belongs:

sim_synset_word(s, w) = max(sim_synset(s, t) for t in synsets_of(w))

Let "top_words(c, k)" be the top "k" most similar words in a sense cluster "c" to the committee which the cluster was based on. The similarity between a synset and a sense cluster is the average similarity between the synset and "top_words(c, k)".

sim_synset_cluster(s, c, k) = sum(sim_synset_word(s, w) for w in top_words(c, k)) / k

We can check if a word belongs to a correct sense cluster by checking the similarity of the sense cluster to all the synsets of the word. If the maximum similarity found is greater than a threshold, then it is considered to be a correct sense cluster.

is_correct_sense_cluster(w, c, k, threshold) = max(sim_synset_cluster(s, c, k) for s in synsets_of(w)) >= threshold

"k" was set to 4 and the threshold was varied in the experiments. By checking for which synset is most similar to the cluster, we can find the WordNet sense that corresponds to the cluster. If more than one cluster corresponds to the same WordNet sense, then only one of them is considered correct. The precision of the system's assignment of a particular word is the percentage of clusters the word is assigned to that are a correct sense. The precision of the system in general is the average precision of all words. The recall of the system's assignment of a particular word is the percentage of correct clusters the word is assigned to out of all the actual correct senses of the word which are found in the corpus. Since the corpus was not sense tagged this is a tricky situation. So the next best thing was to use several clustering methods to extract as many correct sense clusters (as defined above) for the words as possible. The union of all these correct senses are used as the number of different senses each word has in the corpus. The overall recall of the system is the average recall of every word. The precision and recall scores are combined into an F measure:

F = 2 × precision × recall / (precision + recall)

In order to do the evaluation, 144 million words were parsed using Minipar in order to extract the grammatical relations between words. The intersection of words which are found in WordNet with those which in the corpus have a pointwise mutual information with all their features greater than 250 was used for evaluation. 13403 words met the criteria and the average number of features per word was 740.8. As was already mentioned, the threshold described in "is_correct_sense_cluster" was varied during the experiment. Here are the results of the F measure using different thresholds with different clustering algorithms. The clustering algorithms were modified to just give the centroids of the clusters and then associated all words which have a similarity greater than 0.18. This is in order to associate the same word with different clusters for each sense. The system described is "CBC" which stands for "Clustering By Committee" which outperformed all the other clustering methods for all thresholds.