Random Projects

Plotting Google Maps Data

Living a few months in Cambridge: (red is day, blue is night)

I love the blue mass at the astronomy centre and the movie theatre.

Inferring Gaussian Process Autocorrelation

It seems reasonable and intuitive to think that the sample autocovariance function of a stationary gaussian process to be a sufficient statistic for its covariance function, and I read that this is indeed true for certain stationary GPs with a rational spectrum. This condition is quite similar (if not the same) for GPs to possess a state space representation.

It’s also interesting that not all GPs are ergodic. GPs are mixing (and hence ergodic, I believe) if the covariance dies off to zero after a point, or if its spectrum is absolutely continuous. Loosely, this means that the distribution of the process can be inferred from just one long sample. GPs with an exponentiated sine squared (ESS) covariance function, for example, wouldn’t be ergodic.

As a consequence, the covariance function of a zero mean GP with an ESS kernel, and similar signals, cannot be inferred from a single sample. This is reasonable, as no matter how long the signal is, there’s no new information in it after a certain point. Intuitively, a sample from a zero mean GP with an ESS kernel might look like \((3, 2.5, 3, 2.5, ...)\). The ESS is a kernel which is periodic, and the correlation of points spaced half a period apart is closest to zero (compared to any other pair of points), but still strictly positive. Another sample from that GP may look like \((-2, -1.7, -2, -1.7, ...)\).

Points one and two are closer together within each sample than across samples due to the correlation, but given just one observation of the signal (and with no knowledge of the mean of the process), it would appear that points one and two are negatively correlated.

The image below shows this; the black line is the true autocovariance function of the zero mean GP with an ESS kernel, and the boxplots show the sampling distribution of the unbiased sample autocovariance function based on single samples.

Griffin Lim Algorithm and a Minimal Working Implementation

This minimal implementation below is based on the Librosa source on GitHub.

import numpy as np
from scipy.signal import stft, istft
from scipy.io import wavfile as wav

sr, audio = wav.read('audio.wav')

stft_of_audio = stft(audio)[2]
st_spectrum = np.abs(stft_of_audio)

angles = np.empty(st_spectrum.shape, dtype=np.complex64)
angles[:] = np.exp(2j * np.pi * np.random.rand(*st_spectrum.shape)) # angles[:] = 1.0

stft_recon = 0.; momentum = 0.99
for _ in range(2000):
	stft_recon_prv = stft_recon
	signal_recon = istft(st_spectrum * angles)[1]
	stft_recon = stft(signal_recon)[2]
	angles[:] = stft_recon - (momentum / (1 + momentum)) * stft_recon_prv
	angles[:] /= np.abs(angles) + 1e-16

audio_recon = np.array(signal_recon, dtype = np.int16)
audio_recon = np.repeat(audio_recon, 2).reshape(-1, 2)
wav.write("audio.wav", sr, audio_recon)

(my surface-level understanding follows!)

The GLA, in a few words, is an algorithm to find a matrix \(X\) such that the frobenius norm:

\[|X - S \circ S^{-1} (X)|_F\]

(with \(S\) corresponding to the STFT operation and \(S^{-1}\) to the inverse) is minimised (the absolute value of \(X\) - the spectrum - must be known and fixed). This norm is supposed to be zero for “proper” STFTs, but it isn’t when initiated with random phases like in the code above (the STFT has redundant information, so one point may be used in multiple windows - so when initiated with random phases, the matrix might not correspond to anything sensible and may be inconsistent due to the windows overlapping). It is quite interesting to me that, in a probabilistic sense, STFTs of signals would probably only ever admit a density where this norm is zero (i.e. if the STFT was consistent in the overlapping windows). And yet, by minimising this norm, we approach an STFT with phases more or less consistent with the spectrum we’ve got, so this “loss” doesn’t correspond to a likelihood at all.

Reinforcement Learning

I toyed around with the idea of fitting some reinforcement learning models for game balance two years ago but never wrote it up. I saw a Google AI blog post on this recently, which had me digging up some of this old work.

a. Tic-Tac-Toe Q Learning

I wrote up an implementation of an inefficient table-based Q-learning algo that has a dataframe contain all possible states and future states in tic-tac-toe. The (pretty ugly) code is below. Here, a reward is +1 if x wins, -1 if o wins, 0 otherwise.

Since it’s a solved game, the Qs only take the value -1, 0 or 1 after fitting.

Tic-Tac-Toe Q Learning
import numpy as np
import pandas as pd
from itertools import product

def validate(string):
    blanks = string.count(' ')
    xs = string.count('x')
    os = string.count('o')

    return (blanks + xs + os == 9) and \
           ((xs == os) or (xs - os == 1))

boards = list(map(lambda x: ''.join(x), product(*[' xo'] * 9)))
boards = [board for board in boards if validate(board)]

df = pd.DataFrame(dict(boards=boards))

def event(string):
    x_wins = 0; o_wins = 0

    select = lambda x: ''.join([string[i] for i in x])

    winning_lines = (
        string[:3], string[3:6], string[6:9], # rows
        select([0, 3, 6]), select([1, 4, 7]), select([2, 5, 8]), #cols
        select([0, 4, 8]), select([2, 4, 6]) # diag

    for line in winning_lines:
        unique_chars = ''.join(set(line))
        if len(unique_chars) == 1 and unique_chars != ' ':
            if unique_chars == 'x':
                x_wins += 1
            if unique_chars == 'o':
                o_wins += 1

    if x_wins + o_wins == 1:
        if x_wins: return 'x'
        if o_wins: return 'o'
    elif x_wins + o_wins > 2:
        return 'i' # impossible boards

    if string.count(' '):
        return ''
        return 'd'

df['event'] = df.apply(lambda x: event(x['boards']), axis=1)
df = df.loc[df.event != 'i'].reset_index(drop=True)

df['Q'] = np.random.uniform(-0.1, 0.1, size=df.shape[0])
df.loc[df.event == 'd', 'Q'] = 0.0
df.loc[df.event == 'x', 'Q'] = 1.0
df.loc[df.event == 'o', 'Q'] = -1.0
df['Q_next'] = df['Q'].copy()

def next_state(string, return_index=True):
    next_char = 'x' if string.count('x') == string.count('o') else 'o'
    chars = tuple(string)
    next_states = [''.join(chars[:i] + (next_char,) + chars[(i + 1):]) \
                   for i in range(9) if chars[i] == ' ']

    if return_index:
        next_states = list(df.loc[df.boards.isin(next_states)].index)
    return next_states

df['next_state'] = None
df.loc[df.event == '', 'next_state'] = \
    df.loc[df.event == ''].apply(lambda x: next_state(x['boards']), axis=1)

def next_q(string, index):
    qs = np.array(df.Q)[df.loc[index, 'next_state']]
    return qs.max() if string.count('x') == string.count('o') else qs.min()

for _ in range(100):
    df.loc[df.event == '', 'Q_next'] =\
        df.loc[df.event == ''].apply(lambda x: next_q(x['boards'], x.name), axis=1)
    df.loc[:, 'Q'] = df.Q*0.9 + 0.1*df.Q_next

df.loc[df.boards == 'xxo o x  ', 'Q'] # expected reward is 0 because it's o's turn
df.loc[df.boards == 'xxo o x o', 'Q'] # expected reward is 1 because the next turn, x wins

b. DnD / Pokemon Battle Simulations

I coded up some classes in python but haven’t gotten around to coding up the RL side of it.

One neat trick: as moves can be special and do weird stuff, I coded them up as dicts of metadata, with a field named “special_effects” that contains a lambda function that takes in the opponent class instance and applies the weird stuff to it (iirc).

Another point, sometimes one can simplify decisions as other choices may be strictly worse. Then, simple simulation can lead to the probability that a party wins, as a function of input params.

Simulation out of Boredomville
sim = function(eagle_hp=26, dragon_hp=22, eagle_ac=13,
               dragon_ac=16, eagle_str=3, eagle_dex=3,
               dragon_dex=1, dragon_str=2, can_use_bw=T) {

    d = function(n=1, r=20) sample(1:r, n, T)

    adv = function(roll) ifelse(roll == 20, 2, 1)

    eagle_turn = function(ac) {
        dmg = 0; hit = d()
        if(hit + 2 + eagle_str >= ac)
            dmg = dmg + adv(hit)*d(1, 6) + eagle_str
        hit = d()
        if(hit + 2 + eagle_str >= ac)
            dmg = dmg + adv(hit)*sum(d(2, 6)) + eagle_str

    dragon_turn = function(ac) {
        if(can_use_bw) {
            if(sample(0:2, 1)) can_use_bw = F
            save_mod = ifelse(d() + eagle_dex >= 11, 0.5, 1)
            return(floor(save_mod * sum(d(4, 8))))
            hit = d()
            if(hit + 2 + dragon_str >= ac)
                return(adv(hit)*d(1, 10) + dragon_str)
            else return(0)

    if(d() + eagle_dex >= d() + dragon_dex)
        dragon_hp = dragon_hp - eagle_turn(dragon_ac)
    while(T) {
        if(dragon_hp <= 0) return(1)
        eagle_hp = eagle_hp - dragon_turn(eagle_ac)
        if(eagle_hp <= 0) return(0)
        dragon_hp = dragon_hp - eagle_turn(dragon_ac)

mean(replicate(10000, sim()))

Voice Conversion

I got two speakers’ data from the VoxCeleb dataset (about 500 5-second samples, which perhaps is not enough) and tried to fit some voice conversion models. First attempt was a VAE type model, with one encoder that encodes mel-specs from both speakers, and a decoder for each speaker that decodes the mel-specs from some latent variables. I also added a “cycle-consistency loss”. It didn’t work - the converted mel-specs (mel-spec of speaker A -> Encoder -> Decoder_B) weren’t understandable at all, although there was vague feature-transfer going on. I also tried fitting the CycleGAN to this data, but no luck there either. Again, some vague transfer going on but the converted samples just sounded like the source speaker.

New York Conditional Taxi Dropoff Probabilities

I fit a twenty component mixture of multivariate normals, using scikit-learn, to the four dimensional new york taxi pickup/dropoffs dataset.

The dimensions look like (pickup_lat, pickup_lon, dropoff_lat, dropoff_lon). The aim is to predict the distribution of (dropoff_lat, dropoff_lon) by conditioning on (pickup_lat, pickup_lon).

Fancy ways to do this might include fitting a neural net or some kind of a gp to the conditional density, but here, I literally just fit a 4d mvn mixture to the whole dataset. To condition, we just plug in the pickup position and renormalize (Bayes rule).

Envelope Modelling

Google’s Quick Draw dataset contains multiple observations of quickly drawn envelopes. I fit a 256-component restricted boltzmann machine (heavily overparameterised; not a great model - I know) to the data, which represents a big nasty distribution over the random field that represents an envelope image. Now, starting off with a completely random image, using Gibbs sampling, we can make our way to the typical set of the distribution, which hopefully looks like an envelope. Here’s what the burn in looks like:

Inferring the Extent of Differentiability

Let’s say that we have an observation of a noiseless function but we don’t know how smooth it is. You could probably fit a Matern GP with different smoothness parameters to see which parameter maximises the log marginal likelihood (the matern parameter corresponds to the number of times one can differentiate a sample from the gp).

Below, I’ve simulated a Matern GP with a particular parameter, and fit it using parameters ranging from \(\\{0.5, ..., 5\\}\). The color & label correspond to the parameter while sampling.

Modelling Audio using GPs

I used the S-PAD and the GP-PAD models from Richard Turner’s thesis to make these plots using some random audio data from the internet.

Sample stan code for this.
// S-PAD
data {
	int n;     // len
	real x[n]; // audio
parameters {
	real<lower = 0, upper = 1> l;
	real<lower = 0, upper = 10> s;
	vector<lower = -10, upper = 2>[n] sigma;
model {
	sigma[1] ~ normal(0, s);
	sigma[2:n] ~ normal(l * sigma[1:(n - 1)], s*(1 - l^2)^0.5);
	x ~ normal(0, exp(sigma));

data {
	int n;
	int n_s; // nrow of S22
	real seg_a[n]; // audio sample
	matrix[n, n_s] factor; // mvn conditional distribution shift factor S12 * S22^(-1)
	cholesky_factor_cov[n_s] K22c; // cholesky factor of S22
parameters {
	real<lower = -7.5, upper = 10> mu;
	vector<lower = -7.5, upper = 10>[n_s] sigma;
transformed parameters {
	vector[n] sigma_vec;
	sigma_vec = rep_vector(mu, n) +
		factor*(sigma - rep_vector(mu, n_s));
model {
	sigma ~ multi_normal_cholesky(rep_vector(mu, n_s), K22c);
	seg_a ~ normal(0, log(1 + exp(sigma_vec)));

Modelling my d20 dice

I fit a spline on a sphere representing my d20. The color represents the model probabilities of landing on a face. The distance from the centre represents the proportion of times my dice fell on a particular number during six hundred trials (empirical data in other words). I was mainly testing out a plotting idea here (and if it’d be possible to fit splines on sphere - the latter is pretty easy, with a GP, just have the distance be something like a great-circle distance).

Changes in Park-Going

… w.r.t. baseline, as a result of the pandemic (as of 23rd Apr 2020). Based on the Google mobility dataset.


Efficient Gaussian Process Computation

I’ll try to give examples of efficient gaussian process computation here, like the vec trick (Kronecker product trick), efficient toeliptz and circulant matrix computations, RTS smoothing and Kalman filtering using state space representations, and so on.

4 min read

Gaussian Processes in MGCV

I lay out the canonical GP interpretation of MGCV’s GAM parameters here. Prof. Wood updated the package with stationary GP smooths after a request. I’ve run through the predict.gam source code in a debugger, and mainly, the computation of predictions follows:

~1 min read


I wanted to see how easy it was to do photogrammetry (create 3d models using photos) using PyTorch3D by Facebook AI Research.

1 min read

Dead Code & Syntax Trees

This post was motivated by some R code that I came across (over a thousand lines of it) with a bunch of if-statements that were never called. I wanted an automatic way to get a minimal reproducing example of a test from this file. While reading about how to do this, I came across Dead Code Elimination, which kills unused and unreachable code and variables as an example.

1 min read
Back to Top ↑



I used to do a fair bit of astrophotography in university - it’s harder to find good skies now living in the city. Here are some of my old pictures. I’ve kept making rookie mistakes (too much ISO, not much exposure time, using a slow lens, bad stacking, …), for that I apologize!

1 min read

Probabilistic PCA

I’ve been reading about PPCA, and this post summarizes my understanding of it. I took a lot of this from Pattern Recognition and Machine Learning by Bishop.

1 min read

Modelling with Spotify Data

The main objective of this post was just to write about my typical workflow and views rather than come up with a great model. The structure of this data is also outside my immediate domain so I thought it’d be fun to write up a small diary on making a model with it.

6 min read

Morphing with GPs

The main aim here was to morph space inside a square but such that the transformation preserves some kind of ordering of the points. I wanted to use it to generate some random graphs on a flat surface and introduce spatial deformation to make the graphs more interesting.

2 min read

SEIR Models

I had a go at a few SEIR models, this is a rough diary of the process.

4 min read

Speech Synthesis

The initial aim here was to model speech samples as realizations of a Gaussian process with some appropriate covariance function, by conditioning on the spectrogram. I fit a spectral mixture kernel to segments of audio data and concatenated the segments to obtain the full waveform. Partway into writing efficient sampling code (generating waveforms using the Gaussian process state space representation), I realized that it’s actually quite easy to obtain waveforms if you’ve already got a spectrogram.

5 min read
Back to Top ↑


Gaussian Process Middle C

First of my experiments on audio modeling using Gaussian processes. Here, I construct a GP that, when sampled, plays middle c the way a grand piano would.

1 min read
Back to Top ↑


Back to Top ↑