%pylab inline
from scipy.misc import derivative
from scipy.optimize import curve_fit
import inspect

def signal(t, a, f, ph):
    return a * sin(2 * pi * f * t + ph)

p0 = [2, 8, 0]
noise = 0.1

T = arange(0, 1, 0.01)
plot(T, signal(T, *p0), '.-k')
xlabel('time (s)')

labels = inspect.getargspec(signal)[0][1:]
p0dict = dict(zip(labels, p0))
print('labels:', labels)
print('p0dict:', p0dict)

D = zeros((len(p0), len(T)))

# for every parameter
for i, argname in enumerate(labels):
    # for every sample
    for k, t in enumerate(T):
        # Define a function to be derivated.
        # It's our signal at a given time t with all parameters fixed, except the one named 'argname'.
        func = lambda x: signal(t, **dict(p0dict, **{argname: x}))
        
        # Calulate derivative of func at point given by p0
        # Be careful with dx, since it's 1 by default!
        D[i,k] = derivative(func, p0dict[argname], dx=0.0001)

plot(T, signal(T, *p0), '--k', lw=2, label='signal')

for Di, argname in zip(D, labels):
    plot(T, Di, '.-', label=argname)
    
legend(loc='best')
xlabel('time (s)')

I = 1/noise**2 * einsum('mk,nk', D, D) # How cool is that?
print(I)

iI = inv(I)

for argname, variance in zip(labels, iI.diagonal()):
    print('{}: {:.2g}'.format(argname, sqrt(variance)))

S = signal(T, *p0) + randn(T.size) * noise
plot(T, S, '.-k')

popt, pcov = curve_fit(signal, T, S, p0=p0)

for argname, variance in zip(labels, pcov.diagonal()):
    print('{}: {:.2g}'.format(argname, sqrt(variance)))

Tl = linspace(T[0], T[-1], 10000)
plot(Tl, signal(Tl, *popt))

def cramer_rao(model, p0, X, noise, show_plot=False):
    """Calulate inverse of the Fisher information matrix for model
    sampled on grid X with parameters p0. Assumes samples are not
    correlated and have equal variance noise^2.
    
    Parameters
    ----------
    model : callable
        The model function, f(x, ...).  It must take the independent
        variable as the first argument and the parameters as separate
        remaining arguments.
    X : array
        Grid where model is sampled.
    p0 : M-length sequence
        Point in parameter space where Fisher information matrix is
        evaluated.
    noise: scalar
        Squared variance of the noise in data.
    show_plot : boolean
        If True shows plot.
    
    Returns
    -------
    iI : 2d array
        Inverse of Fisher information matrix.
    """
    labels = inspect.getargspec(model)[0][1:]
    p0dict = dict(zip(inspect.getargspec(model)[0][1:], p0))
    
    D = zeros((len(p0), X.size))
    for i, argname in enumerate(labels):
        D[i,:] = [ derivative(lambda p: model(x, **dict(p0dict, **{argname: p})),
                              p0dict[argname],
                              dx=0.0001)
                  for x in X ]
        
    if show_plot:
        plot(X, model(X, *p0), '--k', lw=2, label='signal')
        for d, label in zip(D, labels):
            plot(X, d, '.-', label=label)
        legend(loc='best')
        title('Parameter dependence on particular data point')
    
    I = 1/noise**2 * einsum('mk,nk', D, D)
    iI = inv(I)
    
    return iI

N = []
T0 = arange(0, 11)
for t0 in T0:
    T = arange(-t0, 10 - t0, 0.01)
    N.append( cramer_rao(signal, p0, T, noise)[2,1] )
    
plot(T0, N, 'o-')
axhline(0, color='black')
xlabel('zero location - point with fixed phase')
ylabel('frequency - phase correlation')

T = arange(-0.5, 0.51, 0.01)
cramer_rao(signal, p0, T, noise, show_plot=True);
xlabel('time (s)')