618.sagews

Brief study related to the Weiss-Weinstein bound, by Ahmed Fasih

This brief demo studies the $\eta(s, h)$ function described by Weiss and Weinstein in their 1985 IEEE Trans. Info. Theory paper (vol. IT-31, no. 5, p. 680-682).

We want to estimate the frequency of an exponentially-damped sine-wave (the attenuation factor is proportional to the frequency) in noise from deliberately undersampled data. We use aliasing to give us large peaks in the likelihood far away from the true frequency.

Assume uniform priors on $\theta$. Then, give its bounds and a vector containing it. This vector should be fairly fine because we'll do trapezoidal integration over it.

import numpy

theta_min = 0.
theta_max =  40.
theta = numpy.linspace(theta_min,theta_max, 250)

Next, set up the signal's sample rates and times. We deliberately undersample.

fs=10.
Tmax = 1.
t=numpy.arange(0, Tmax, 1/fs)

These are our test points $\lbrace h \rbrace$. We let them be a grid for plotting purposes, and let their extents be different from those of $\theta$ because we can.

hs = numpy.arange(0, 25, 0.1)

This is our true signal, $s(t)$, an exponentially-damped sine wave. The full data model would thus be

$$x(t)=\sin(2\pi\theta t) \exp(-0.02 \theta t) + \text{AWGN}$$

[Nota bene: Numpy's sin() will return a matrix or an array depending on what type the time vector "t" is. The distinction between array and matrix types is more important here than in Matlab (where it is non-existant) because * does element-wise multiply for arrays, and does matrix product for matrixes.]

print "Make sure 't' has type array, so * does element-wise multiply:", type(t)
s = lambda x: numpy.sin(t*x*2.*numpy.pi) * numpy.exp(-t*x*.02)

Make sure 't' has type array, so * does element-wise multiply: <type 'numpy.ndarray'>

Here, we assume $s=\frac{1}{2}$ and unit noise variance (that's why there's no "sigma_n" variable), and evaluate $\eta(h)$ for the range of $h$ specified previously. Each evaluation thereof involves an integral:

$$\begin{align}\eta(h) &= \int_\Theta \sqrt{p(\theta+h) p(\theta)} \int_\Omega \sqrt{p(x|\theta) p(x|\theta+h)} dx d\theta \\

&\propto \int_\Theta \exp \left( -\frac{1}{4 \sigma_n^2} || s(\theta+h) - s(\theta) ||^2 \right) d\theta

\end{align}$$

from numpy.linalg import norm

def eta(h):
    return numpy.trapz( 
    [numpy.exp(-1/4.*norm(s(th+h)-s(th))**2.) for th in theta] , # NB: this is a list!
    theta)

all_eta = [eta(h) for h in hs] # NB: this is a list!

Having stored $\eta(h, s=\frac{1}{2})$ in the above list, we plot below. The reason we point out the lists is that they are not arrays or matrixes, and arithmetic operators on them will fail. [We could convert them all to arrays with "numpy.array(...)" but both numpy.trapz() above and pylab.plot() below work fine with lists. As a sidenote, these functions accept a variety of inputs like lists, matrixes, arrays, etc., because they are all members of the same class family and are iteratable.]

from matplotlib import rc
rc('text', usetex=True)

import pylab
pylab.figure()
pylab.plot(hs, all_eta, '.-')

pylab.xlabel('$h$', fontsize=18)
pylab.ylabel('$\propto \eta(h, s=0.5)$', fontsize=18)
pylab.title('$\propto \eta(h, s=0.5)$, with uninformative priors', fontsize=18)
pylab.savefig('etas', dpi=72)