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.