Saturday, September 8, 2012

Functional Programming


Today I was hacking on my TDT reading code when I came across the following ugly code:
if map(iscomplex, map(np.squeeze, map(np.array, args))):
    # ...
I thought to myself, "Hey why not throw some functional programming at this?", to get
if composemap(iscomplex, np.squeeze, np.array)(args):
    # ...
Definitely more readable if a little too terse.

What's great is that I can pass an arbitrary number of functions to compose map and it will run each function on all the arguments in args. Sweet.

How is it actually implemented?
def compose2(f, g):
    if not all(map(callable, (f, g))):
        raise TypeError('f and g must both be callable')
    return lambda *args, **kwargs: f(g(*args, **kwargs))
This little gem just returns a function that calls f and g with any number of arguments. Totally jacked from the Python functools module docs.
from functools import partial, reduce
from itertools import repeat

def composemap(*args):
    partial_map = map(partial, repeat(map, len(args)), args)
    return reduce(compose2, partial_map)
That's it! Now, to be sure there's a ton of stuff going on here and it took me a while to grok it, so let's break it down line by line.

The partial function takes a callable and optional positional and keyword arguments and returns another callable that has the given arguments already filled in. It's similar to currying, but a little different.

The first line of composemap results in
[partial(map, arg0), partial(map, arg1), ... partial(map, argN)].

This makes sense because what I want to do is actually run map, but with the functions I've provided as arguments first and whatever sequence i decided to map each of the functions I've provided over. However, I don't want to evaluate this right away so I create a list (actually a tuple) that contains all of the partial applications that I want to perform when I provide a sequence to map over.

OK whew! That was a lot to take in.

For the last line recall what partial does: it takes a callable with part of the sequence arguments that you want to provide and returns a callable with those arguments already set to go when you call the function with the final arguments (or more arguments).

Look at the reduce documentation for more details on reduce. If you're familiar with Haskell you can think of reduce as the Python version of foldl.

What I've done here is give reduce the list of partial functions that I want to eventually evaluate.

What the end result looks like is something like this

Given a function $h = f\circ g$ and an $n$-sequence of partial functions $\mathbf{p} = \left(p_1, p_2,\ldots,p_n\right)$ then the "reduce" function, in this case is given by
\begin{equation}
\textrm{reduce}\left(h, \mathbf{p}\right) = h\left(\ldots h\left(p_1, p_2\right), \ldots p_n\right)
\end{equation}

Recall that $h\left(f,g\right)$ returns a function. So running $h$ with $p_1,p_2$ as arguments gives $p_0\left(p_1\left(\cdot\right)\right)$. Pass that into $h$ with the third element of $\mathbf{p}$, $p_3$. Keep doing that until all of the arguments are eaten up.


















Sunday, February 26, 2012

Definition of a Confidence Interval

"95% of intervals of the form $$\left({\overline X} - \frac{2\sigma}{\sqrt{n}},{\overline X} + \frac{2\sigma}{\sqrt{n}}\right),$$ which are based on samples of size $n$ each, will contain the mean of the normal density in their interior."

Saturday, February 11, 2012

Convexity

Some things to prove:
  1. The distance to a set $S$ from some vector $\mathbf{x}$, given by $d\left(\mathbf{x}, S\right) = \displaystyle\min_{\mathbf{y}\in S}\,\,\|\mathbf{x} - \mathbf{y}\|$ is convex.
  2. The set of of all symmetric positive definite matrices is convex.
  3. Let $f\left(x\right)=\exp\left(Q\left(x\right)\right)$. Prove that $f$ is log concave if $Q\left(x\right)$ is concave.

To do this, I'm going to employ the use of some tools from Chapter 2 of the book Convex Optimization, by Steven Boyd and Lieven Vandenberghe. One of the great things about this book is that it is very accessible.

First, recall the definition of lines and line segments: If two points $x,y$ in $\mathbb{R}^{n}$ are not equal, then points of the form the form $y = \theta x_1 + \left(1 - \theta\right)x_2$ where $\theta \in \mathbb{R}$ form the line passing through $x_1$ and $x_2$.

This is nice and all, but it's not exactly the most intuitive way to define a line. The authors graciously give another interpretation which I find to be very intuitive: $$y=x_2+\theta\left(x_1-x_2\right).$$ They explain it in an equally intuitive way: "... $y$ is the sum of the base point $x_2$ (corresponding to $\theta = 0$) and the direction $x_1 - x_2$ (which points from $x_2$ to $x_1$) scaled by the parameter $\theta$." (p. 21)

Let's look at this in vector notation. $y$ should actually be $\mathbf{y}$ because we're not talking about scalar values here, we're talking about points/vectors on the $n$-dimensional real line, which have more than one component. The same goes for $x_1$ and $x_2$. So, it should be written $$\mathbf{y} = \mathbf{x}_2 + \theta\left(\mathbf{x}_1 - \mathbf{x}_2\right).$$ It's one of my math pet peeves when people don't use something to indicate whether they are talking about vectors or scalars, even if you define it that way from the start. We're so used to working with symbols like $x$ and $y$ as scalar values, that it's nice to be able to let the paper take care of the notation and not have to think "Am I looking at a scalar or vector here?". When vectors are bold you can immediately see that the data type is a vector, whereas with the other convention $x$ and $y$ could be a scalar or a vector depending on how it was defined at the beginning of the [proof, book, paper, etc.]. OK, rant over!

Assuming you have NumPy and matplotlib installed you can use the following Python code to demonstrate the definition of a line i just described.
from pylab import *

def randbetween(a, b, *size):
    assert a < b, 'a must be less than b'
    for s in size:
        assert isinstance(s, int), \ 
            'the elements of size must all be integers'
    bma = b - a
    return bma * rand(*size) + a

b = 10
a = 1
x1 = randbetween(a, b, 1, 2)
x2 = randbetween(a, b, 1, 2)

# 0 <= theta <= 1
theta = linspace(0, 1, 10000)[newaxis]

# compute the line over theta
y = dot(theta.T, x1) + dot(1 - theta.T, x2)

# plot the points and the line
figure()
hold('on')
plot(x1[0, 0], x1[0, 1], '.', ms=20)
plot(x2[0, 0], x2[0, 1], '.', ms=20)
plot(y[:, 0], y[:, 1])
hold('off')

# show the plots
show()
You can literally just copy this code, then in an IPython console type: paste. This will run all of the code I just listed and should result in a nice plot. Obviously for $\theta > 1$ and $\theta < 0$ the line will extend beyond $\mathbf{x}_1$ and beyond $\mathbf{x}_2$ respectively. What's more, you can rewrite the second equation in vector-matrix form. Given two $k$-dimensional row vectors $\mathbf{x}_1$ and $\mathbf{x}_2$ and given a column vector $\boldsymbol\theta = \left(\theta_1, \theta_2, \dotsc, \theta_n\right)$ where $\theta_1 = 0$ and $\theta_n = 1$ and monotonically increasing elements, then a line in $\mathbb{R}^{n}$ can be written as $$\mathbf{Y} = \boldsymbol\theta\cdot\mathbf{x}_1 + \left(\mathbf{1}_n - \boldsymbol\theta\right)\cdot\mathbf{x}_2,$$ where $\mathbf{1}_n$ is an $n\times{1}$ vector with each element equal to 1, and $\mathbf{Y}$ is an $n\times{k}$ matrix with each column being the points of the line in the $k$th dimension. The next important concept is the notion of an affine set. A set $C\subseteq{\mathbb{R}^{n}}$ is said to be affine if the line through any points in $C$ not equal to each other is in $C$. Said differently, "$C$ contains the linear combinations of any two points in $C$, provided the coefficients in the linear combination sum to one (Boyd & Vandenberghe)." (p. 22). In symbolic terms, $$x_1, x_2\in{C}\textrm{ and }\theta\in{\mathbb{R}}\Rightarrow{\theta{x_1}+\left(1-\theta\right)x_2\in{C}}$$ OK, now on to convex sets. The definition is exactly the same as an affine set except that $\theta$ is bounded by 0 and 1. A set $C$ is convex if, $$x_1,x_2\in{C}\textrm{ and }\theta\in\left[0,1\right]\Rightarrow{\theta{x_1}+\left(1-\theta\right)x_2\in{C}}$$ Let's now define a convex function. A function $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is said to be convex if the domain of $f$ is a convex set and if for all $x$, $y$ is an element of the domain of $f$, and $\theta$ with $0\leq\theta\leq{1}$ we have, $$f\left(\theta x + \left(1-\theta\right)y\right)\leq\theta f\left(x\right) + \left(1-\theta\right)f\left(y\right).$$ Strict convexity is when $x\neq{y}$ and $0\lt\theta\lt{1}$ and concavity is when $-f$ is convex, with strict concavity if $-f$ is strictly convex.

Let's use the definition of convexity to show that $d$ is convex.

The triangle inequality and homogeneity show us that $$f\left(\theta\mathbf{x} + \left(1-\theta\right)\mathbf{y}\right)\leq f\left(\theta\mathbf{x}\right) + f\left(\left(1-\theta\right)\mathbf{y}\right)=\theta f\left(\mathbf{x}\right) + \left(1-\theta\right)f\left(\mathbf{y}\right)$$ is true. Then substituting $d$ for $f$ allows us to see that $d$ is convex.

In my next few posts I'll show the 2 and 3 parts of the list at the beginning of this post.

Uninstall All Ruby Gems

Today I ran into an issue when trying to do a nice Rails tutorial. Having two versions of Ruby installed on my machine was (I think) causing some conflicts. At best, it was confusing me. So, here's how to uninstall all Ruby gems:

gem list | cut -d" " -f1 | xargs gem uninstall -aIx

Tuesday, February 7, 2012

Gamma

I was working on some mathematical statistics problems and I came across a slick estimation problem that I thought I would share.

Problem:
Given $$\left(f\mid k\right) = \frac{x^{k-1}\exp\left(-x/\theta\right)}{\Gamma\left(k\right)\cdot\theta^{k}},\,\,\,\textrm{ where } x \gt 0,\,\,\, k \gt 0,\,\,\, \theta = 1,$$Determine whether it is possible to find a value of $c$ such that $d\left(X\right) = cX^{2}$ will be an unbiased estimator of $k$.

My attempt at a solution:

$$E\left[c\cdot X^{2}\right] = k$$ for the left hand side to be an unbiased estimator of $k$. $$c\cdot E\left[X^{2}\right] = k$$ is true via the properties of $E\left[\cdot\right]$. Rearranging terms we can see that $$c = \frac{k}{E\left[X^{2}\right]}$$ $$c = k\cdot\left(\int_{0}^{\infty}x^{2}\frac{x^{k-1}e^{-x}}{\Gamma\left(k\right)}\mathrm{d}\,\,x\right)^{-1} = k\cdot\left(\int_{0}^{\infty}\frac{x^{k+1}e^{-x}}{\Gamma\left(k\right)}\mathrm{d}\,\,x\right)^{-1}$$ $$k\cdot\Gamma\left(k\right)\cdot\left(\int_{0}^{\infty}x^{k+1}e^{-x}\mathrm{d}\,\,x\right)^{-1} = \frac{\Gamma\left(k+1\right)}{k\cdot\Gamma\left(k+1\right)}=\frac{1}{k}.$$ So the value of $c$ such that $E\left[cX^{2}\right]$ is NOT an unbiased estimator of $k$ is $1/k$ because by definition, an unbiased estimator of some parameter $\theta$ cannot depend on the value of $\theta$.

Saturday, February 4, 2012

Setting up Git for my website

Just wanted to log the process of setting up the Git version control system for editing the code of my website.

  • Local Machine
    • Make a directory for your code
    • Run git init inside of it
    • Add some content and commit it
  • Server
    • Password stuff
      • Make an SSH public encryption key if you haven't already
        • ssh -t rsa -C 'someone@somewebsite.com'
        • Follow the instructions it spits out
        • Give a password if you want
      • ssh -p 2222 me@myserver.com 'mkdir -p .ssh'
        • You'll have to enter a password here
      • cat ~/.ssh/id_rsa.pub | ssh -p 2222 me@myserver.com 'cat >> .ssh/authorized_keys'
    • Git stuff
      • ssh -p 2222 me@myserver.com
      • cd public_html && mkdir website.git && cd website.git && git init --bare
      • cd hooks && cp post-receive.sample post-receive
      • echo 'GIT_WORK_TREE=/dir/where/you/want/to/store/code git checkout -f' >> post-receive
  • Local Machine
    • git remote add website ssh://ip.add.re.ss:portnumber/~/public_html/website.git
    • git push -u website +master:refs/heads/master && git push
Now, if I want to edit the code of my website on my local machine I can do so and have changes pushed to the website whenever I want to. I just have to be careful not to push breaking changes. Time to get intimate with Git branches!

Funny

Saw this picture while trolling some links from Hacker News.

#4 reminds me of some arguments I've had...