June 28, 2013

Enabling formulas and code highlighting on Blogspot

This blog uses MathJax to embed LaTeX-style formulas and Google Code Prettify for syntax highlighting.

To install these tools select the Template tab from the blog admin panel and click Edit HTML. Insert the following code just after the opening head tag.
<script src='https://google-code-prettify.googlecode.com/svn/loader/run_prettify.js'/>    
  <script type='text/x-mathjax-config'>//<![CDATA[
    MathJax.Hub.Config({
      extensions: ["tex2jax.js"],
      jax: ["input/TeX", "output/HTML-CSS"],
      displayAlign: "left",
      styles: {
        ".MathJax_Display": {
          //"background-color": "rgb(230,230,230)",
          "padding-left": "4em",
          "float": "left",
          "display": "inline"
        }
      },
      tex2jax: {
        inlineMath: [ ['$','$'] ],
        displayMath: [ ['$$', '$$'] ],
        processEscapes: true
      },
      "HTML-CSS": { availableFonts: ["TeX"] }
    });
  //]]></script>
  <script src='http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=default' type='text/javascript'>
</script>   
The in posts you'll be able to use the following syntax to highlight code
<pre class="prettyprint">
int f(int x)
{
  return 2 * x;
}
</pre>
and \$ and \$\$ delimiters to insert LaTeX math inline or in display style.

A couple of off topic tips. To make the text area adapt to the browser window, in the same editor search for the text .content-outer and edit the code around it as follows:
body {
min-width: 100px /*$(content.width);*/
}
.content-outer, .content-fauxcolumn-outer, .region-inner {
min-width: 100px /*$(content.width);*/
max-width: 2000px /*$(content.width);*/
_width: 100% /*$(content.width);*/
}
Also, in math jax, it's possible to define macros as follows:
TeX: {
  Macros: {
    cov: '{\\operatorname{cov}}',
    reals: '{\\mathbb{R}}',
    T: '{\\mathrm{T}}',
    tr: '{\\operatorname{tr}}',
    paren: ['\\left( #1 \\right)', 1],
    bracket: ['\\left[ #1 \\right]', 1],
    brace: ['\\left\\{ #1 \\right\\}', 1]
  }
}

June 25, 2013

Pyson, something between JSON and Python

JSON, by virtue of its minimalist nature, lacks a few features which make it difficult to efficiently represent some instances of structured data.

In particular, JSON only allows the representation of trees, i.e. each node has exactly one parent node with exception of the root element. What this means is that data cannot be reused in multiple places within a document.

In addition, JSON's type system is very limited - there is no way to define custom types. This means that most structures must be represented as dictionaries, making it necessary to repeat the schema (i.e. the dictionary keys) along with every instance of the data (i.e. the dictionary values), or requiring ugly mechanisms like special keys (e.g "__type__") to give each dictionary a label indicating how to interpret the data.

So the idea is to extend JSON with a couple of additional light-weight constructs to solve these issues. In the same original spirit which gave rise to JSON, I decided to borrow some syntax from Python, of which JSON is a subset (kind of). So what we have is
$$
\text{JSON} \subset \text{Pyson} \subset \text{Python}
$$

Syntax

The extra syntactic features in Pyson are:
  • The root element is no longer any value but a document, i.e. a list of key-value pairs of the form key = value:
    x = 3
    y = {'a': 3, 'b': 2}
    z = [true, false]
    
  • It is possibile to reference values which have been defined previously
    x = [1,2,3]
    y = {'z': x}
    u = [y, x, y]
    
  • It is possibile to define new types and instantiate them
    DateTime = Type("DateTime", ["year", "month", "day"])
    today = DateTime(2013, 6, 24)
    
Parsing and semantic analysis

Modifying a JSON parser to accommodate the new syntax is very straightforward. It's more complicated to get the reference semantics right.

In JSON, there is no difference between the parsing and semantic analysis because JSON values are a one-to-one representation of the abstract syntax tree (AST) generated by a parser.
$$
\text{JSON string} \quad\to_{\text{parser}}\quad \text{AST} = \text{data}
$$
This also means that converting a value stored in memory to a JSON string is immediate.
$$
\text{AST} = \text{data}\quad \to_{\text{to string}}\quad \text{JSON string}
$$

References, on the other hand, need to be resolved in the AST before it is a useful in-memory representation of the data.
$$
\text{Pyson string} \quad\to_{\text{parser}}\quad\text{AST}\quad\to_\text{linker}\quad \text{data}
$$
Serialization also becomes more complicated, because it's necessary to analyze the data, factor out common nodes and substitute in the references. It's also necessary that the variables are output in the correct order as required by the references.
$$
\text{data} \quad\to_\text{unlinker}\quad \text{AST} \quad\to_{\text{to string}}\quad \text{Pyson string}
$$

Note that when a Pyson string is parsed, linked, unlinked and converted back into a string, the output will be the same as the input (except for some irrelevant differences in variable order). When, however, the data being serialized has not been come from a Pyson string, it is possible that there are common nodes which don't have an explicit variable name - hence they must be generated automatically.

Examples: reading Pyson with C#

We'll read in this Pyson document
DateTime = Type("DateTime", ["year", "month", "day"])
today = DateTime(2013, 6, 24)
expiry = DateTime(2014, 6, 1)
period = [today, expiry]
The following C# code steps through parsing, linking, unlinking and conversion to string
List<KeyValuePair<string, object>> inputAST = PysonParser.Parse(inputString);

Dictionary<string, object> data = PysonParser.Link(inputAST);

List<KeyValuePair<string, object>> outputAST = PysonParser.Unlink(data);

string outputString = PysonParser.ToString(outputAST);
The output string comes out exactly the same as the input string. We'll rarely be interested in the intermediate ASTs so we can alternatively just write
Dictionary<string, object> data = PysonParser.Link(inputString);

string outputString = PysonParser.ToString(data);
From C# we can easily access the data using the dynamic keyword.
dynamic data = PysonParser.Link(inputString);
double startYear = data["period"][0]["year"];
double endYear = data["period"][1]["year"];
bool test1 = data["period"][0] is PysonInstance;
bool test2 = data["DateTime"] is PysonType;
bool test3 = data["period"][0].Type == data["DateTime"];

Examples: generating Pyson with C#

Here's an example of how to generate a Pyson string from an in-memory object.
var dateTime = new PysonType("DateTime", new object[] { "year", "month", "day" });
var today = dateTime.New(2013, 6, 24);
var expiry = dateTime.New(2014, 6, 1);

dynamic data = new Dictionary<string, object>();
data["period"] = new object[] { today, expiry };
 
string outputString = PysonParser.ToString(data);
The output isn't as readable as we might expect:
Type1 = Type("DateTime", ["year", "month", "day"])
period = [Type1(2013, 6, 24), Type1(2014, 6, 1)]
This is because the Dictionary object doesn't directly contain the DateTime type, so the rather ugly name 'Type1' is generated automatically. This can be fixed by adding it to the Dictionary:
var dateTime = new PysonType("DateTime", new object[] { "year", "month", "day" });
var today = dateTime.New(2013, 6, 24);
var expiry = dateTime.New(2014, 6, 1);

dynamic data = new Dictionary<string, object>();
data["period"] = new object[] { today, expiry };
data["DateTime"] = dateTime;
 
string outputString = PysonParser.ToString(data);
which generates:
DateTime = Type("DateTime", ["year", "month", "day"])
period = [DateTime(2013, 6, 24), DateTime(2014, 6, 1)]
Note that objects which appear more than once in the data are automatically factored out and placed before the first reference.
var dateTime = new PysonType("DateTime", new object[] { "year", "month", "day" });
var today = dateTime.New(2013, 6, 24);
var expiry = dateTime.New(2014, 6, 1);

dynamic data = new Dictionary<string, object>();
data["period"] = new object[] { today, expiry };
data["DateTime"] = dateTime;
data["simulation"] = new Dictionary>string, object<();
data["simulation"]["start"] = today;
 
string outputString = PysonParser.ToString(data);
and the output
DateTime = Type("DateTime", ["year", "month", "day"])
obj3 = DateTime(2013, 6, 24)
period = [obj3, DateTime(2014, 6, 1)]
simulation = { "start": obj3 }

May 21, 2013

Using COM constants in Python

The Python library win32com allows you to instantiate and manipulate COM objects transparently from Python. One aspect which is not automatic is the treatment of constants, for example the following code will generate an error unless you've performed some necessary steps.
import win32com.client
xl = win32com.client.Dispatch(“Excel.Application”)
print win32com.client.constants.xlWorkbookNormal
To make this work you first need to run the script makepy.py passing the relevant COM type in the command line.
C:\Python27\Lib\site-packages\win32com\client>makepy.py "Excel.Application"
After this step the above Python code should print -4143.

April 26, 2013

Visual Studio-style docking windows with DockPanel Suite

Installation

These instructions work with Visual Studio 2010.

Download the binaries from the DockPanel Suite website. In Visual Studio, select Choose Items from the Toolbox window's context menu (it may be hidden, display it with View > Toolbox). In .NET Framework Components click Browse and select the WeifenLuo.WinFormsUI.Docking.dll assembly downloaded from the website. This will add the DockPanel component to the Toolbox.

Usage

The container form must be an MDI form containing a DockPanel component. Add a new form to your Visual Studio project, and set its IsMdiContainer property to True. Place a DockPanel component on it and make it fill up the form's client area by setting its Dock property to Fill.

To add docking windows to the container form you can instantiate them in code and associate them to the container form by calling the Show method specifying the DockPanel as the owner.

public partial class Form1 : Form
{
  Form2 formA = new Form2();
  Form2 formB = new Form2();
  Form2 formC = new Form2();
 
  public Form1()
  {
    InitializeComponent();

    formA.Show(this.dockPanel1);
    formB.Show(this.dockPanel1);
    formC.Show(this.dockPanel1);
  }
}

Here's the result.





April 20, 2013

The Metropolis algorithm applied to Bayesian analysis

As the dimensionality of a Bayesian model increases, numerically integrating the parameter space becomes computationally intensive and, in all but the most premeditated cases, analytic computation is unfeasible. An alternative is to sample the posterior distribution
$$p(\theta|D) = \frac{p(D|\theta) p(\theta)}{\int p(D|\theta)p(\theta) d\theta}$$
via Monte-Carlo. The Metropolis algorithm, first discussed in Metropolis, Rosenblush, Teller (1953), allows you to sample from any distribution $p(x)$ providing you can compute a function $Cp(x)$ proportional to $p(x)$. In this case, the algorithm allows us to sample the posterior without performing the integral in the denominator.

The algorithm

We start by defining a transition probability
$$M(x,x') = q(x,x')\min\left(\frac{p(x')}{p(x)}, 1\right)$$
where $q$ is any arbitrary transition probability. In practice, to create a chain of samples based on $M$ we
  1. Sample $q(x, \cdot)$ to choose $x'$
  2. Move to $x'$ if $p(x') > p(x)$, else move to $x'$ with probability $p(x')/p(x)$, else remain in $x$.
As the number of iterations increases, the probability of being in state $x$ converges to $p(x)$. Intuitively, this happens because $p$ is an eigenfunction of $M$ with eigenvalue $1$. Further details can be found in Robert, Casella (2004).

Practical issues

Although the algorithm works in principle with any non-degenerate proposal distribution $q$, the optimal choice minimizes the burn-in period, i.e. the number of initial iterations we discard before considering the chain to have 'converged' to the target distribution.

April 17, 2013

The simplest possible Bayesian model

In chapter 5, Kruschke goes through the process of using Bayes' rule for updating belief in the 'lopsidedness' of a coin. Instead of using R I decided to try implementing the model in Excel.

The Model

Coin tosses are IID with a $\mathrm{Bernoulli}(\theta)$ distribution, thus the probability of observing $N$ of which $z$ are heads is
$$p(z|\theta, N) = \theta^z (1-\theta)^{(N-z)}$$
Using Bayes' rule we get
$$p(\theta|z,N) = \frac{\theta^z (1-\theta)^{N-z} p(\theta)}{p(z|N)}$$
Note that if $p(\theta)$ is already of the form
$$f(\theta;a,b) = \frac{\theta^{a-1} (1-\theta)^{b-1}}{B(a,b)}$$
where the denominator is a normalizing constant then $p(\theta|N) = f(\theta;a+z,b+N-z)$.

Excel Implementation

This is very straight-forward using the BETADIST function, which is cumulative, so we divide the unit interval into equally spaced subintervals, calculate the probability of $\theta$ being in the interval and update on the basis of the new observations.

Here's the result:


Discrete model

Using a discrete distribution for $\theta$ rather than a continuous parametric form we can allow for much more freedom in the choice of initial prior distribution. The update rule becomes
$$p(\theta|z_N,N)=\frac{p(z_N|z_{N-1},N-1,\theta)p(\theta|z_{N-1},N-1)}{\sum_{i=1}^n p(z_N|z_{N-1},N-1,\theta_i) p(\theta_i|z_{N-1},N-1)}$$
in other words
$$p_{N,j} = \frac{\tilde p_{N,j}}{\sum_{i=1}^n \tilde p_{N,i}}$$
where
$$\tilde p_{N,k} = (\theta_k y_N + (1-\theta_k)(1-y_N)) p_{N-1,j}$$
and $y_N = z_N - z_{N-1}$ is the result of the $N$th coin toss.

As expected the result is very similar to the continuous model.

In the discrete model we can assign any initial prior. Here is an example with an initial prior that assigns non-zero probability to only 3 values of $\theta$ (22.5%, 52.5% and 77.5%) and random tosses generated with a real theta of 20%. The distribution quickly converges to probability 1 that $\theta=22.5\%$.


April 12, 2013

Doing Bayesian Data Analysis

Assuming I can keep at it, I'll be making my way through Kruschke's Doing Bayesian Data Analysis. Here's a few concepts he goes through in Chapter 4.

The Bayes factor

This is a ratio which allows you to compare which out of two models best fits the data. By introducing a binary parameter which determines the choice of model, Bayes' rule
$$p(M_i|D)=\frac{p(D|M_i)p(M_i)}{p(D)}$$
gives us the Bayes factor
$$\frac{p(M_1|D)}{p(M_2|D)}=\frac{p(D|M_1)}{p(D|M_2)}\frac{p(M_1)}{p(M_2)}$$
or also
$$\frac{p(D|M_1)}{p(D|M_2)}=\frac{p(M_1|D)}{p(M_2|D)}\frac{p(M_2)}{p(M_1)}$$
Note that if we have no prior preference between the two models, these ratios are equal. The Bayes' factor is useful because on their own $p(D|M_i)$ scale on the basis of the model - the more complex the model, the lower the absolute value - but this is actually a good thing because by considering the ratio we automatically trade off complexity (when there's little data) against descriptive value (when the simpler model doesn't fit the data).

Data order invariance

I don't get this part - it's always true that $p(\theta|D', D)$ is invariant of the order in which Bayes' rule is applied, by definition. The factorized rule cited at the end is just the result of applying the definition of data independence given the parameters.

The problem with Bayesian mathematics

Computing the posterior given some new data usually means performing an integral, which, even using approximations, can be computationally intensive (especially since the posterior will be fed into an optimizer). Numerical integration on the other hand would impose a limit on the dimensionality of the model. Thus we winning method is Markov Chain Monte Carlo (MCMC).