Snippet: List Comprehension vs Set AND operation

timeit.timeit('a = range(50);b=range(25, 75);set(a)&set(b)')
# takes about 5 seconds
timeit.timeit('a = range(50);b=range(25, 75);[x for x in a if x in b]')
# takes about 15 seconds

Interesting Plot Snippet: Bitwise Or and Addition Operation

If the sum and the result of bitwise or operation of two numbers are the same, plotting these numbers in a scatter diagram produces a beautiful shape.

import matplotlib.pyplot as plt

l1 = []
l2 = []
for i in xrange(1000):
 for j in xrange(1000):
  if (i + j) == (i | j):
   l1.append(i)
   l2.append(j)

plt.scatter(l1,l2)
plt.show()

This is the figure that comes up:

Value Iteration in Gridworld, Reinforcement Learning

import operator
import copy
import math


class State:
    def __init__(self, rowid, colid):
        self.row = rowid
        self.col = colid

    def __str__(self):
        return str(self.row) + " " + str(self.col)


def get_value_from_state(state):
    return state_values[state.row][state.col]


def get_reward_from_state(state):
    return reward_values[state.row][state.col]


def move_to_new_state(state, action):
    if action == 'left':  # state will only change columnwise
        if state.col == 0 or (state.col == 2 and state.row == 1):  # then we can't move left
            return copy.deepcopy(state)
        else:
            return State(state.row, (state.col - 1))
    if action == 'right':  # state will only change columnwise
        if state.col == 3 or (
                        state.col == 0 and state.row == 1):  # can't move right when on the right boundary or at the first column
            return copy.deepcopy(state)
        else:
            return State(state.row, (state.col + 1))
    if action == 'up':
        if state.row == 0 or (state.row == 2 and state.col == 1):
            return copy.deepcopy(state)
        else:
            return State(state.row - 1, state.col)
    if action == 'down':
        if state.row == 2 or (state.row == 0 and state.col == 1):
            return copy.deepcopy(state)
        else:
            return State(state.row + 1, state.col)


def moveWithProbabilites(state, action):
    if action == 'left' or action == 'right':
        thisstate = move_to_new_state(state, action)
        upstate = move_to_new_state(state, 'up')
        downstate = move_to_new_state(state, 'down')
        return {thisstate: .8, upstate: .1, downstate: .1}

    elif action == 'up' or action == 'down':
        thisstate = move_to_new_state(state, action)
        leftstate = move_to_new_state(state, 'left')
        rightstate = move_to_new_state(state, 'right')
        return {thisstate: .8, leftstate: .1, rightstate: .1}


state_values = [
    [0, 0, 0, 1],
    [0, -100, 0, -1],  # -100 means we can't go this way
    [0, 0, 0, 0]
]
reward_values = [  # in each state, the agent gets a reward of -.03
    [-.03, -.03, -.03, -.03],
    [-.03, -.03, -.03, -.03],
    [-.03, -.03, -.03, -.03]
]
gamma = 1

states = [
    State(0, 0),
    State(0, 1),
    State(0, 2),
    State(1, 0),
    State(1, 2),
    State(2, 0),
    State(2, 1),
    State(2, 2),
    State(2, 3),
]
actions = ['left', 'right', 'up', 'down']

prev = state_values[0][0]

if __name__ == "__main__":
    while True:
        for s in states:
            mx = -9999999
            for a in actions:
                dic = moveWithProbabilites(s, a)
                tempSum = get_reward_from_state(s)  # to be added to the values later

                tempSum += gamma * dic.items()[0][1] * get_value_from_state(dic.items()[0][0])  # .9 * .1 * v(s)
                tempSum += gamma * dic.items()[1][1] * get_value_from_state(dic.items()[1][0])
                tempSum += gamma * dic.items()[2][1] * get_value_from_state(dic.items()[2][0])

                if tempSum > mx:
                    mx = tempSum

            state_values[s.row][s.col] = mx
            
        if abs(state_values[0][0] - prev) < .00003:
            break
        prev = state_values[0][0]

print state_values

Reinforcement Learning with Pygame Tutorial Part - 2

(Part 1: http://anixtech.blogspot.com/2016/12/reinforcement-learning-with-pygame-part.html)

I have created the second part of the tutorial as a video. It could be seen here:
https://www.youtube.com/watch?v=o5GiQkClbAY

The full code is also available in github:
https://github.com/hasanIqbalAnik/q-learning-python-example

Reinforcement Learning With Pygame Part - 1

Following is an example of a simple game which could be used to train agents. 

import pygame as pg
from pygame.locals import *
import sys
import random


def new_rect_after_action(newrct, act):
    if act == 'right':
        if newrct.right + newrct.width > windowWidth:
            return newrct
        else:
            return pg.Rect(newrct.left + newrct.width, newrct.top, newrct.width,
                           newrct.height)  # Rect(left, top, width, height)
    else:  # action is left
        if newrct.left - newrct.width < 0:
            return newrct
        else:
            return pg.Rect(newrct.left - newrct.width, newrct.top, newrct.width,
                           newrct.height)  # Rect(left, top, width, height)


def circle_falling(crclradius):
    newx = random.randint(0 + crclradius, 200 - crclradius)
    multiplier = random.randint(1, 4)
    newx *= multiplier
    return newx


windowWidth = 800
windowHeight = 400
FPS = 3  # frames per second setting
fpsClock = pg.time.Clock()

pg.init()

window = pg.display.set_mode((windowWidth, windowHeight))  # width, height
pg.display.set_caption('Catch the ball!')

# setup colors
RED = (255, 0, 0)
GREEN = (0, 255, 0)
WHITE = (255, 255, 255)
BLACK = (0, 0, 0)

# specify circle properties
crclCentreX = 50
crclCentreY = 50
crclRadius = 20

# specify rectangle properties
rctLeft = 400
rctTop = 350
rctWidth = 200
rctHeight = 50
rct = pg.Rect(rctLeft, rctTop, rctWidth, rctHeight)  # Rect(left, top, width, height)

action = 'left'
score = 0
font = pg.font.Font(None, 30)

while True:

    for event in pg.event.get():
        if event.type == QUIT:
            pg.quit()
            sys.exit()
        elif event.type == pg.KEYDOWN:
            if event.key == pg.K_LEFT:
                action = 'left'
                rct = new_rect_after_action(rct, action)
            elif event.key == pg.K_RIGHT:
                action = 'right'
                rct = new_rect_after_action(rct, action)

    window.fill(BLACK)
    # at this position, the rectangle should be here. else loses
    if crclCentreY >= windowHeight - rctHeight - crclRadius:
        if rct.left <= crclCentreX <= rct.right:
            score += 1
        else:
            score -= 1
        crclCentreX = circle_falling(crclRadius)
        crclCentreY = 50
    else:
        crclCentreY += 20

    pg.draw.circle(window, RED, (crclCentreX, crclCentreY),
                   crclRadius)  # circle(Surface, color, pos(x, y), radius, width=0)
    pg.draw.rect(window, GREEN, rct)  # rect(Surface, color, Rect, width=0)

    text = font.render('score: ' + str(score), True, (238, 58, 140))
    window.blit(text, (windowWidth - 100, 10))
    pg.display.update()
    fpsClock.tick(FPS)

This will create a game window that looks like the following screenshot. The green rectangle can be moved with left and right arrow of the keyboard. If the player can place it under the falling red circle, he gets a point, else loses one. That's all there is in this game. In the second part of this tutorial, we will add intelligence and train an agent to play this game.

Update: Second part is now available: http://anixtech.blogspot.com/2016/12/reinforcement-learning-with-pygame.html

Simplified Neural Network Backpropagation Example Python

This is an even simplified version of this fantastic tutorial: https://mattmazur.com/2015/03/17/a-step-by-step-backpropagation-example/
Let's consider a very simple neural network shown in the following image:

All we want is that given the input .1, the network would output the value .9.
We will be using $\frac{1}{1+e^{-x}}$ as our activation function. Initially both weight w1 and w2 are initialized to .5, though this can be any arbitrary value. However, different initialization shall converge to different final values.

First let's consider the net input in the hidden layer(neth). It's simply the input multiplied by w1. That is:(Here, $outh$ means the output of the hidden layer, $neto$ means the net input in the output layer, $outo$ is the output of the output layer. $E$ is the error. )

$neth = i1 * w1 = .1 * .5 = .05$
$outh = \frac{1}{1+e^{-neth}} = \frac{1}{1+e^{-.05}} = .51249$

$neto = outh * w2 = .51249 * .5 = .25624$
$outo = \frac{1}{1+e^{-neto}} = \frac{1}{1+e^{-.25624}} = .563711$

$E = \frac{1}{2}*(target - outo)^2 =  \frac{1}{2}*(.9 - .563711)^2 = .056545$

This concludes our forward pass of the network. We've selected that particular error function for our convenience, as it's convex and differentiable. Now, we can see that our target value was $.9$ but we've got $.5637$ . We are going to update our weights $w1, w2$ so that the network produces that desired output. How would we do that is the reason for using the famous Backpropagation algorithm. The basic idea is, we will calculate how much responsibility should each of the weights take for producing that error. Then the weight would be modified accordingly and the next iteration would begin.

First let's see what happens in case of $w2$. The partial differentiation expression $\frac{\partial E}{\partial w2}$ says the rate of change in the error with respect to $w2$. As it can't be calculated directly, we would use the chain rule to find it's constituting parts.

$\frac{\partial E}{\partial w2} = \frac{\partial E}{\partial outo} * \frac{\partial outo}{\partial neto} * \frac{\partial neto}{\partial w2}$


Let's find out each of those parts separately:

$\frac{\partial E}{\partial outo} = \frac{\partial}{\partial outo}(\frac{1}{2}*(target - outo)^2)  = \frac{1}{2} * 2 * ( target - outo) * (-1) $
$= (.9 - .563711) * (-1) = -.336289$

$\frac{\partial outo}{\partial neto} = \frac{\partial }{\partial neto} (\frac{1}{1+e^{-neto}}) = \frac{-1}{(1+e^{-neto})^2} * e^{-neto} * (-1) = \frac{e^{-neto}}{(1+e^{-neto})^2}$
$  = outo * (1-outo) = .563711 * ( 1- .563711) = .24594$

$\frac{\partial neto}{\partial w2} = \frac{\partial (w2*outh)}{\partial w2} = outh = .51249$

So finally we have found the responsibility of $w2$ for producing the error $E$:

$\frac{\partial E}{\partial w2} = -.33628 * .24594 * .51249 = -.042386$

Now let's turn our attention to $w1$, that is we want to calculate $\frac{\partial E}{\partial w1}$. Again, chain rule comes to our rescue.

$\frac{\partial E}{\partial w1} = \frac{\partial E}{\partial outo} * \frac{\partial outo}{\partial neto} *\frac{\partial neto}{\partial outh} * \frac{\partial outh}{\partial neth} * \frac{\partial neth}{\partial w1}$

We already have calculated $\frac{\partial E}{\partial outo}$ and $\frac{\partial outo}{\partial neto}$. Let's calculate the rest.

$\frac{\partial neto}{\partial outh} = \frac{\partial (w2*outh)}{\partial outh} = w2 = .5$
$\frac{\partial outh}{\partial neth} = outh * ( 1- outh) = .51249 * ( 1- .51249) = .24984$
$\frac{\partial neth}{\partial w1} = \frac{\partial (w1 * input)}{\partial w1} = input = .1 $

So,

$\frac{\partial E}{\partial w1} = -.336289 * .24594 * .5 * .24984 * .1 = -.001033$

Now it's time to update our initial guesses of $w1$ and $w2$. Let's say our learning rate $\eta = .9$
$w1 = w1 - (\eta * \frac{\partial E}{\partial w1}) =  .5 + .9 * 001033 = 0.5009298654$
Similarly,
$w2 = w2 - ( \eta * \frac{\partial E}{\partial w2}) = .5 + .9 * 042386 = 0.538148$

We can continue this process until our network's output is very close to our target. The following is a python implementation of that process:

from __future__ import division
import math

ctr = 0
eta = .9
w1 = .5
w2 = .5

i1 = .1
target = .9

while True:
    ctr += 1
    neth = i1 * w1
    outh = 1 / (1 + math.exp(-neth))

    neto = w2 * outh
    outo = 1 / (1 + math.exp(-neto))

    e = .5 * math.pow((target - outo), 2)
    print 'error now is: ', e
    if e < .00005:
        print '*' * 100
        print 'target is: ', target
        print 'output is now: ', outo
        print 'w1 and w2 is: ', w1, ',', w2
        print 'number of iterations: ', ctr
        break
    # update w2
    de_by_douto = (-1) * (target - outo)
    douto_by_dneto = outo * (1 - outo)
    dneto_by_dw2 = outh

    de_by_dw2 = de_by_douto * douto_by_dneto * dneto_by_dw2

    # update w1
    # de_by_dw1 = de_by_douto * douto_by_dneto * dneto_by_douth * douth_by_dneth * dneth_by_dw1
    dneto_by_douth = w2
    de_by_outh = de_by_douto * douto_by_dneto * dneto_by_douth
    douth_by_dneth = outh * (1 - outh)
    dneth_by_dw1 = i1

    de_by_dw1 = de_by_outh * douth_by_dneth * dneth_by_dw1

    # updated w1 and w2
    w1 = w1 - (eta * de_by_dw1)
    w2 = w2 - (eta * de_by_dw2)

Sample output:

Who Liked Your Last N Number of Facebook Posts Most

Previously it was possible to find out this information for any friend's profile too. But Facebook has recently revoked that access and apps now need to have this access directly from that friend. However, you can find out this info about yourself in the following way:

1. Go to this link: https://developers.facebook.com/tools/explorer/.
2. Click on 'Get Token' > 'Get User Access Token' > check user_posts > 'Get Access token'.
3. This would ask you for the permission to read your timeline. Click on 'Continue as your-name'.
4. You would get an access token like the following image. Copy the whole string.

5. Replace the 'your-token' part in the following code with that string.

from __future__ import division
import operator
import requests
from facepy import GraphAPI

large_dic = {}
token = 'your-token'
graph = GraphAPI(token)
n = 50  # number of posts in the history that we want to consider
ctr = 0
pcount = 0
res = graph.get('me/feed')

def update_progress(progress):
    print '\r{0}% [{1}]'.format(int(progress), '#'*int(progress/10)),

def merge_two_dicts(dict1, dict2):
    """ merges and sums up dict1 with the values from dict2. if key2 from dict2 is present in dict1, then dict1[key2] = dict1[key2] + dict2[key2]
    else it appends key2 with the value as it is
    example: if
    d1 = {'a': 1, 'b': 1, 'c': 1} and
    d2 = {'a': 4, 'd': 1}
    then it returns {'a': 5, 'c': 1, 'b': 1, 'd': 1}
    :param dict1, dict2
    :return merged dict1 and dict2
    """
    for key2 in dict2:
        if key2 in dict1:
            dict1[key2] = dict1[key2] + dict2[key2]
        else:
            dict1[key2] = dict2[key2]
    return dict1


def get_likes_list(post_id):
    """given a post id, this method queries the graphapi to get details of the post, gets the likes list
    and adds them into a dictionary by using the person's name as the key
    :param post_id
    :return dictionary containing likers
    """
    global graph
    global pcount
    d = {}  # to hold the liker's name with a count of 1

    likes_obj = graph.get(post_id + '/likes')
    while True:
        try:
            for liker in likes_obj['data']:
                d[liker['name']] = 1
            likes_obj = requests.get(likes_obj['paging']['next']).json()
        except KeyError:
            pcount += 1
            update_progress((pcount/n)*100)
            break
    return d


for item in res['data']:
    merge_two_dicts(large_dic, get_likes_list(item['id']))

while True:
    if ctr == (n - 25) / 25:  # as 25 is the default limit of /me/feed, to get the last n posts, set the value of ctr to (n-25)/25
        break
    try:
        res = requests.get(res['paging']['next']).json()
        for item in res['data']:
            merge_two_dicts(large_dic, get_likes_list(item['id']))
        ctr += 1
    except KeyError:
        print 'end of calculating posts'


sorted_x = sorted(large_dic.items(), key=operator.itemgetter(1))

print '\n'
print '*'*100
print 'total posts considered', pcount
print 'total likes generated', sum(large_dic.values())
print 'average likes per post', sum(large_dic.values()) / pcount
print '*'*100

for item in reversed(sorted_x):
    print item[0], 'liked your posts', item[1], 'times'

6. Change the value of n in the program to your preference, this is the number of posts that we want to consider while counting the likes.

7. Run that program, it might take a while depending on the response time of facebook server. The output would contain the list of likes in a sorted fashion. With the most liker at the top.

Sample output:

Comparing Radix, Merge, Selection and Bubble Sort

Implementation of Radix Sort:

bucket = {}

# creating bucket for each digits
def reset_bucket():
    global bucket
    for i in xrange(10): bucket[i] = []
# to handle variable length items`
def radix_group_numbers(numbers):
    joined = []
    max_length = max(len(str(x)) for x in numbers)
    bkt = {}
    for i in xrange(max_length + 1):
        bkt[i] = []
    for item in numbers:
        bkt[len(str(item))].append(item)
    for k, v in bkt.items():
        if v: joined.extend(radix_sort(v))
    return joined

def put_in_buckets(numbers, index):  # assuming that the numbers are of equal lenght
    global bucket
    try:
        for item in numbers:
            c = str(item)[index]
            if c is '0':
                bucket[0].append(item)
            elif c is '1':
                bucket[1].append(item)
            elif c is '2':
                bucket[2].append(item)
            elif c is '3':
                bucket[3].append(item)
            elif c is '4':
                bucket[4].append(item)
            elif c is '5':
                bucket[5].append(item)
            elif c is '6':
                bucket[6].append(item)
            elif c is '7':
                bucket[7].append(item)
            elif c is '8':
                bucket[8].append(item)
            elif c is '9':
                bucket[9].append(item)
    except IndexError:
        print 'error occured'
        print numbers, index

def radix_sort(lst):
    global bucket
    max_length = max(len(str(x)) for x in lst)
    for i in xrange(max_length):
        reset_bucket()
        put_in_buckets(lst, max_length - i - 1)
        lst = []
        for key, value in bucket.iteritems():
            if value is not None:
                for item in value:
                    lst.append(item)
        bucket.clear()
    return lst

Implementation of Merge Sort:

# merge two lists
class stack:
    def __init__(self, lst=None):
        if lst is None:
            self.items = []
        else:
            self.items = lst

    def pop(self):
        return self.items.pop()
    def peek(self):
        return self.items[len(self.items) - 1]
    def __str__(self):
        return str(self.items)
    def empty(self):
        return len(self.items) is 0
    def reverse(self):
        l = []
        while self.empty() is False:
            l.append(self.items.pop())
        return l

def merge(l1, l2):  # takes two sorted lists and merges them while keeping the order
    ll = []  # long list consists of l1 + l2
    l1.reverse()
    l2.reverse()
    st1 = stack(l1)
    st2 = stack(l2)
    while st1.empty() is False and st2.empty() is False:
        i1 = st1.peek()
        i2 = st2.peek()
        if i1 <= i2:
            ll.append(i1)
            l1.pop()
        else:
            ll.append(i2)
            l2.pop()
    if st1.empty() is True and st2.empty() is False:
        ll = ll + st2.reverse()
    elif st2.empty() is True and st1.empty() is False:
        ll = ll + st1.reverse()

    return ll


def merge_sort(lst):
    if len(lst) is 1:
        return lst
    else:
        mid = len(lst) // 2
        left = merge_sort(lst[:mid])
        right = merge_sort(lst[mid:])
        return merge(left, right)

Implementation of Selection Sort:

from utils import get_numbers_list

def selection_sort(unsorted_lst):
    for k in xrange(len(unsorted_lst) - 1):  # at max this much iteration required to fully sort

        current_minimum = unsorted_lst[k]  # first elem as the current minimum
        current_minimum_idx = k

        for i in xrange(k, len(unsorted_lst), 1):
            if current_minimum > unsorted_lst[i]:
                current_minimum = unsorted_lst[i]
                current_minimum_idx = i
        # swap the current minimum with the kth element
        tmp = unsorted_lst[k]
        unsorted_lst[k] = current_minimum
        unsorted_lst[current_minimum_idx] = tmp
    return unsorted_lst

Bubble Sort:

from utils import get_numbers_list

def bubble_sort(lst):
    max_iter = len(lst)
    counter = 0
    for m in xrange(max_iter):
        for i in xrange(max_iter - 1 - counter):  # -counter to execute faster
            if lst[i] > lst[i + 1]:  # should swap
                tmp = lst[i + 1]
                lst[i + 1] = lst[i]
                lst[i] = tmp
        counter += 1
    return lst

Now let's generate some sample data:

import sys
import random

# n is the number of samples generated
def generate_samples(n):
    lst = []
    while(n > 0):
        lst.append(random.randint(0, 1000000)) #random integer between 0 and a million
        n -= 1
    f = open("data.txt", "w+")
    for item in lst:
        f.write(str(item) + " ")
    f.close()
    print 'samples written at data.txt...'

def get_numbers_list(fname):
    num_arr = []
    f = open(fname, "r")
    for line in f:
        num_arr.append(list(int(x) for x in line.strip().split(" ")))
    return num_arr[0]

Then we measure the time required for their executions:

import time

from utils import generate_samples, get_numbers_list
from bubble_sort import bubble_sort
from merge_sort import merge_sort
from radix_sort import radix_sort, radix_group_numbers
from selection_sort import selection_sort

#first we need to generate a hugh amount of data

generate_samples(10000)

lst = get_numbers_list('data.txt')

start_time = time.time()
radix_group_numbers(lst)
end_time = time.time()

print 'radix sort took', end_time - start_time, 'seconds'

start_time = time.time()
merge_sort(lst)
end_time = time.time()

print 'merge sort took', end_time - start_time, 'seconds'

start_time = time.time()
bubble_sort(lst)
end_time = time.time()

print 'bubble sort took', end_time - start_time, 'seconds'

start_time = time.time()
selection_sort(lst)
end_time = time.time()

print 'selection sort took', end_time - start_time, 'seconds'

Sample outputs: for 1000 data points:

radix sort took 0.00797295570374 seconds
merge sort took 0.0257380008698 seconds
bubble sort took 0.0774970054626 seconds
selection sort took 0.0248939990997 seconds

for 10000 data points:

radix sort took 0.0430598258972 seconds
merge sort took 0.187104940414 seconds
bubble sort took 8.21462607384 seconds
selection sort took 2.59621596336 seconds

Fining all possible subsets of a given set

Three possible way to create all subsets of a given set is shown below:

1. Using Binary:

lst = ['a', 'b', 'c', 'd', 'e']
def interpret_combination(s):
 s = s[::-1] #reversing the string to match the index positions with the given list
 if not '1' in list(s): # all zeros
  return []
 else:
  indices_of_ones = [i for i, x in enumerate(list(s)) if x == '1'] # find the indices of 1s only
  l = list(s) #convert to list as you can't operate on strings
  for item in indices_of_ones:
   l[item] = lst[item] # replace the one's with the characters on the list
  l = [x for x in l if x != '0'] # replace the zeros
  return ''.join(map(str, l)) # join them to make a string again
all_subsets = []

for i in xrange(2 ** len(lst)): # until 2^(lenght of list)
 l1 = len(str(bin(i))[2:]) # binary representation without the preceding 0b
 l2 = len(lst)
 s = (l2 - l1)*'0' + str(bin(i))[2:] # prepending the zeros in front to make it more readable
 all_subsets.append(interpret_combination(s)) 

print all_subsets

2. Recursive:

def adder(elem, elem_set): # copy elem with each element of elem_set
 exp_lst = []
 for item in elem_set:
  exp_lst.append(elem + item)
 return exp_lst

def all_subsets(s): 
 if len(s) is 0:
  return ''
 elif len(s) is 1:
  return ['', s[0]]
 else:
  return all_subsets(s[:-1]) + adder(s[len(s)-1], all_subsets(s[:-1]))  

print all_subsets(['a','b', 'c', 'd'])  

3. Iterative:

lst = ['a', 'b', 'c', 'd', 'e']
sub_arr = []
sub_arr_tmp = []

for item in lst:
 for elem in sub_arr:
  sub_arr_tmp.append([item + elem])
 sub_arr.append(item)
 if len(sub_arr_tmp) is 0:
  sub_arr.append('')
 else:
   for e in sub_arr_tmp:
    if(isinstance(e, list)):
     sub_arr.append(e[0])
 del sub_arr_tmp[:]

print set(sub_arr)
print len(set(sub_arr))

Finding Gaussian Distribution Parameters Using Expectation Maximization Matlab

Expecation Maximization is a powerful approach to estimate parameters of statistical distributions. It consists of two steps, i) Expectation step, where the likelihood of the dataset is estimated using the current set of parameters. ii) Maximization step, where the parameters are recalculated based on the likelihood found in the expectation step.

Let's see an example: Suppose we have a dataset which has only one feature, namely, age. Suppose the data was taken from a country where most of the people are young and from another country where most of them are old. We are interested in finding out the mean and variance of these two distributions, assuming that they were normally distributed. 

Say our dataset looks like this: 

$[25, 27, 27, 28, 28, 30, 100, 105, 105, 107, 107, 108]$ . We shall implement the EM algorithm to find out the means and variances.

clc;clear all;
% sample dataset
d = [25 27 27 28 28 30 100 105 105 107 107 108];

%initial guesses of the two means and variances
mu1 = 10;
mu2 = 50;
sigma1 = 10;
sigma2 = 10;

% initial responsibilities
p1 = .5;
p2 = .5;

%until converges

for k = 1:100
    
% start expectation step    

v1 = p1 * normpdf(d, mu1, sqrt(sigma1)); 
v2 = p2 * normpdf(d, mu2, sqrt(sigma2));

r1 = v1 ./ (v1 + v2); %divide ith element of v1 by the sum of ith element of both v1 and v2
r2 = v2 ./ (v1 + v2) ;%divide ith element of v2 by the sum of ith element of both v1 and v2

% end expectation step    

% start maximization step

mu1 = sum(r1 .* d) / sum(r1); %recalculate mu1, weighted average of the datapoints
mu2 = sum(r2 .* d) / sum(r2); %recalculate mu2, weighted average of the datapoints

sigma1 = sum((d-mu1) .* (d-mu1) .* r1) / sum(r1); %recalculate weighted variance first distribution
sigma2 = sum((d-mu2) .* (d-mu2) .* r2) / sum(r2); %recalculate weighted variance second distribution

p1 = sum(r1) / length(d); % recalculate responsibilities, total weights normalized by number of data points
p2 = sum(r2) / length(d);

% end maximization step

end

mu1 % 27.5000, mean of the first distribution
mu2 % 105.3333, mean of the second distribution
sigma1 % 2.2500, variance of the first distribution
sigma2 % 6.8889, variance of the second distribution

Basic 2D, 3D Function Plotting in Matlab

There are a number of ways you can plot a function in matlab. Let's see some examples:

1. Let's start with one variable. Plot of the function $f(x) = x^2 -3$ would be:

fplot(@(x) x^2 - 3, [-5 5]) //fplot requires limits

2. Plotting the sum of the outcomes of two dice throw:
We know that the outcome of a dice throw is uniformly distributed over the range [1,6]. So the sum of the outcomes should be in the range [2, 12]. If we continue this experiment many times, we will start to notice that the output is not uniformly distributed, it's actually normal. We can see the simulation by running the command :

hist(randi([1 6], 1, 10000)+randi([1 6], 1, 10000))%adding two uniformly distributed r.v of the range(1,6) 10000 times

Which should produce an image resembling a normal distribution:

2. Drawing a circle from the equation ( $x^2 + y^2 = 5^2 $):

f = @(x) sqrt(5^2 - x.^2)
fplot(@(x) sqrt(5^2 - x.^2), [-5 5])
hold on
fplot(@(x) (-1)*sqrt(5^2 - x.^2), [-5 5])
hold off
axis equal %to preserves axis ratio, else would end up looking like an ellipse 

3. Let's consider a function $z = x^2 + y^2$, for a range of values of x and y:

[x, y] = meshgrid(linspace(-10, 10, 100), linspace(-10, 10, 100))
z = x.^2 + y.^2
mesh(x,y,z)

4. Another example:

[x, y] = meshgrid(linspace(-10, 10, 100), linspace(-10, 10, 100))
xcosx = x.*cos(x)
ysiny = y.*sin(y)
z = xcosx' * ysiny
surfc(x, y , z)

Interesting Plots

1. General shape of an even element Magic Square:

plot(magic(16))

2. General shape of an odd element Magic Square:

plot(magic(15))

3.
x = -100:100;
y = -100:100;
[X, Y] = meshgrid(x, y);
figure
surf(x, y, X./Y)

4.  //2 b contin

Solving a simple Linear Regression problem using CVX and Matlab

Caution: Contains over-(simplification, generalization, and some other *tions)

Homer is a hawker who sells batman toys. Kids love those very much and his business is on the rise. He noticed that for the past five days, these are his costs and profits:

Items Sold	Profit
100	960
80	760
70	660
150	1460
90	860

Now, a competitor, Ned, has stolen these data and wants to know a little more about his daily operating expense and per item price. Ned is a very good statistician and feels that these data should form a nice regression problem expressed in the following equation:

y = b0 + b1*x

Where 'y' is the value of profit, 'b0' is the fixed cost that Homer incurs everyday and 'b1' is the per price item and 'x' is the number of items Homer sells everyday.

Ned recently bought a license of Matlab and downloaded cvx, a nice modelling tool and thought to give it a try to find out Homer's secrets.

He writes the following matlab code:

b = [960 760 660 1460 860]; % profit vector
x = [100 80 70 150 90]; % item sold each day

cvx_begin
    variables b1 b0 % b1 = per item price, b0 = fixed cost per day
    minimize(norm(((b1*x) + b0) - b )) % minimizing the estimated and actual values
cvx_end

And yes, he finds that per day cost is 40 and per item price is 10. (For Algebraist: ya, I know)

The Last Dumpling: Another Game Of Probability

BBT? Sheldon? Leonard? Rings a bell?

Anyways, let's say, Sheldon and Leonard are two very good friends. While sharing a dinner, Sheldon notices that only one dumpling is remaining on the plate. Being his friend, he doesn't want to upset Leonard but also wants that delicious dumpling!

So, he proposes a game.  As an assertive character, he sets the rules. He and Leonard would roll two dices. If in both of the two dices, a 1, 2, 3 or 4 comes up, Leonard wins the last dumpling. If any of them contains a 5 or 6, Sheldon wins. Leonard greedily accepts the terms, thinking that he would have a much higher probability of winning it.

But we all know Sheldon wouldn't want to play if he doesn't have a more than half chance of winning. So what do you think? Who gets the last dumpling?

Ok, let's see what happens when two dices are rolled.

The probability that one of the dices would have the value between 1-4 is the summation of their individual probability of happening. Let's call the outcome of the roll X. So it would be:

P(Dice 1 being 1-4)  = P(1<=X<=4) = P(X=1) + P(X=2)  + P(X=3) + P(X=4) = 1/6 + 1/6 + 1/6 + 1/6 [as each of them has an equal chance of happening] = 4/6

The second dice would have the same probability of this happening. So both of them having this outcome is the product of the probabilities. It's simple counting principle. If one event can happen in 5 ways and another can happen in 6 ways, together they can happen in 5*6 = 30 ways.

P(Both Dices being 1-4) = P(Dice 1 being 1-4) * P( Dice 2 being 1-4) = 4/6 * 4/6 = 16/36 = 4/9.

So if the game is played 9 times, Leonard would win 4 of the times and rest of the times Sheldon would win. Meaning, probability of Leonard's winning is about 44% and Sheldon's is 56%.

Bad luck Leonard!

Chuck-a-Luck: A fine game involving probability

I stumbled onto this game while reading a fantastic book, 'Innumeracy' by John Allen Paulos. The rules of the game is simple: you play against a host. The  host has three dices with number 1-6 on them. Then you would be asked to choose a number in that range. Let's say you have chosen number 4. The host would then roll the three dices at the same time. Now there are three outcomes:

1. If all three dices have 4 on them, you get 3 dollars.
2. If two of the dices come up with 4, you get 2 dollars.
3. If only dice has 4 on it, you get 1 dollar.
4. Else(none of the dices has came up with 4) then you pay the host 1 dollar.

Now, what do you think about your chances are in this game? Obviously, as it is played in casinos, the casino owners won't have allowed this game if there were no profit for them. So even if this game seems quite winnable, you should re-calculate your chances.

Here is the expected outcome of this game:

expected outcome(EO) = sum(all the ways you can win * associated rewards) - sum(all the ways you can lose * associated cost)

EO = sum(P(all three dices has 4) * 3 + P(any two dices has 4) * 2 + P(any dices has 4) * 1) - P(no dice has 4) * 1

= (3C3 * (1/6)^3 * (5/6)^0)*3 + (3C2 * (1/6)^2 * (5/6)^1)*2 + (3C1 * (1/6)^1 * (5/6)^2)*1) - (3C3 * (5/6)^3 * (1/6)^0)*1

= 3/216 + 10/72 + 25/72 - 125/216 = -.08

As you can see, the outcome is quite counter intuitive. It seems that you would lose .08 dollars every time you play the game.

Following is a simulation of the game play. It shows how many times you can play before you lose all of your money in the cold-hearted probabilistic eventuality:

'''
code is also hosted at https://github.com/hasanIqbalAnik/probabilities/blob/master/monty_hall_simul.py
'''
import random
def dice_roll():
    return random.randint(1, 6)

balance = 10 #you start with 10 dollars
count_games = 0
choice = random.randint(1,6)

while(balance > 0):
    outcomes = [dice_roll(), dice_roll(), dice_roll()]
    if(outcomes.count(choice) == 3):
        balance = balance + 3
    elif(outcomes.count(choice) == 2):
        balance = balance + 2
    elif(outcomes.count(choice) == 1):
        balance = balance + 1
    else:
        balance = balance - 1
    count_games = count_games + 1

print balance # 0 :(
print count_games # sample output, 74 times

This also matches the probability calculation that you would be able to play around 80 times with 10 dollars, losing .08 dollars on average each time.

HITS Algorithm Calculation Example

The equations for HITS algorithm are:

a = L'h (L' denoting transpose of adjacency matrix L)
h = La

a, h denoting the authority score and hub score respectively.

The following python code shall compute the hub score and authority score for this particular graph:

import numpy as np

adjacency_mtx = np.matrix([

  [0, 0, 0, 1, 0],
  [0, 0, 0, 1, 1],
  [0, 0, 0, 1, 0],
  [0, 0, 0, 0, 0],
  [0, 0, 0, 0, 0]
 ])
hub_score = np.matrix([ #initial guess
  [2],
  [2],
  [2],
  [2],
  [2]
 ])
authority_score = np.matrix.transpose(adjacency_mtx)*hub_score
hub_score = adjacency_mtx*authority_score

print authority_score
print hub_score

#output:
[[0]
 [0]
 [0]
 [6]
 [2]]

[[6]
 [8]
 [6]
 [0]
 [0]]

This matches our intuition that, node 4 has a higher authority over node 5 and node 2 is a better hub than nodes 1 and 3.

Maximum Independent weight and traceback set in Graph

The following code shall find the maximum weight independent set and weight in a given graph:

def wis(arr):
 if(len(arr) == 0):
  return 0
 A = []
 B = {}
 for i in xrange(len(arr)+1):
  A.insert(i, 0)
 A[0] = 0
 A[1] = arr[0]
 B[arr[0]] = [arr[0]]
 for i in xrange(2, len(arr)+1): 
  A[i] = max(A[i-1], A[i-2]+arr[i-1])
  try:
   if(A[i-2]+arr[i-1] > A[i-1]):
    B[A[i]] = [B[A[i-2]], arr[i-1]]
  except KeyError:
   B[A[i-2]] = 0
   B[A[i]] = [B[A[i-2]], arr[i-1]]
 return A,B

res, trace = wis([5, 7, 100, -15, 3, 2])
print res[len(res)-1]
print trace[res[len(res)-1]]

Beauty of Dynamic Programming

Simply printing the 37th term of Fibonacci using recursion took me 10 seconds(8 GB, Core i5). Tried to find the 50th term, was waiting for like eternity.

def fib(n):
 if n<=2:
  f = 1
 else:
  f = fib(n-1) + fib(n-2)
 return f

print fib(37) #10sec, 24157817

But,

def fibo(n):
 fib = {}
 for k in range(1, n+1):
  if k<=2:
   f = 1
  else:
   f = fib[k-1] + fib[k-2]
  fib[k] = f
 return fib[k]
 

print fibo(100)#instant, 354224848179261915075

Even the 100th Term comes along instantaneously!

Artificial Neural Network in Octave: Backpropagation to predict test scores

This post shall be using the same code of Programming assignment 4(week 5) in the online course:
https://www.coursera.org/learn/machine-learning/
As Backpropagation is a mathematically complex algorithm, using a simpler dataset and reviewing each step along the way would give us a better intuition. That is the goal of this post.

Here, we shall be taking the following steps.

1. Generate some random data points. Let's say, we are going to appear for the GRE test. As a part of the preparation, we would appear for two preliminary tests of the same format, let's say powerprep1(p1) and powerprep2(p2). We would try to predict the final score based on these two. I have written a function to generate test scores as much as needed in a random fashion.
https://gist.github.com/hasanIqbalAnik/6aa2af7138595d2ba85d

Here, p1, p2 would consist our input matrix X, and final scores would be Y. These are the first two and the last columns of the data matrix returned by our data generating function. For example:

P1	P2	Final
317	318	319
305	306	307
302	303	303

2. Structure of our Neural Network: For simplicity, It would have 3 layers: 1 input layer(3 nodes(including bias)), 1 hidden layer(6 nodes(including bias) and 1 output layer.

The full code is available here:
https://gist.github.com/hasanIqbalAnik/bd51dbf3e91550c69620

Now, to fit this dataset in the code, we need to care about just the following things:

• Our num_labels would be 340.
• We do not have pre-initialized Theta1 and Theta2 as we did in the assignment, so we would need to initialize them randomly from the beginning. 
• Handle lambda carefully, a higher value of lambda would result in less overfitting but high bias and vice versa.

The rest can just be left as it is. The prediction performance would depend on a number of things like, how your data is distributed, number of hidden layers, your handling of bias and variance etc.

Natural Language Processing with Python: Chapter 6 Excercise Answers

#ex1
#Too descriptive

#ex 2
def gender_features(word):
    return {
        'suffix1': word[-1:],
        'suffix2': word[-2:],
        'startswith': word[0].lower(),
        'length':len(word),
        'first2char':word[0:2].lower(),
        'containsyn':'yn' in word
        }

names = ([(name, 'male') for name in names.words('male.txt')] +
[(name, 'female') for name in names.words('female.txt')])
random.shuffle(names)
featuresets = [(gender_features(n), g) for (n,g) in names]

train_set = nltk.apply_features(gender_features, names[500:])
devtest_set = nltk.apply_features(gender_features, names[500:1000])
test_set = nltk.apply_features(gender_features, names[1000:len(names)])

classifier = nltk.NaiveBayesClassifier.train(train_set)

p(nltk.classify.accuracy(classifier, test_set)) # 81% accuracy
p(classifier.show_most_informative_features(5))