In [1]:
# gdapy06_00 clears all variables and import various modules
# you must run this cell first, before running any of others

# note that these examples use inverse theory techniques that
# are only covered in later chapters of the book

# History:
# 2023/04/24 - Created by W. Menke from Matlab version
# 2024/05/23 -  Patches to accomodate new verions of numpy, scipy, matplotlib
#               patches dot product to return scalar value

# A reset deletes any previously-created variables from memory
%reset -f

# import various packages

import numpy as np                    # matrices & etc
from matplotlib import pyplot as plt  # general plotting 
from matplotlib import patches        # plotting shapes
from datetime import date             # date and time arithmetic
from math import exp, pi, sin, cos, tan, sqrt, floor, log10, nan, log   # math functions
import scipy.linalg as la             # linear algebra functions
import os                             # operating system
import matplotlib.cm as cm
from matplotlib.colors import ListedColormap
import scipy.signal as sg
import scipy.stats as st
import scipy.sparse.linalg as las
from scipy import sparse
import matplotlib

# function to make a numpy N-by-1 column vector
# c=eda_cvec(v1, v2, ...) from a list of several
# array-like entities v1, v2, including a number
# a list of numbers, a tuple of numbers, an N-by-0 np array
# and a N-by-1 np array. The function
# also insures that, if v is an np array, that
# c is a copy, as contrasted to a view, of it
# It promotes integers to floats, and integers
# and floats to complex, by context.
# This version concatenates many argments,
# whereas c=eda_cvec1(v1) takes just one argiment.
# I recommend always using eda_cvec(v1, v2, ...)
def gda_cvec(*argv):
    t = int;
    Nt = 0;
    for a in argv:
        v = gda_cvec1(a);
        N,M = np.shape(v);
        Nt = Nt + N;
        if( N==0 ):
            continue; # skip vector of zero length
        if (t==int) and isinstance(v[0,0],float):
            t=float;
        elif isinstance(v[0,0],complex):
            t=complex;
    w = np.zeros((Nt,1),dtype=t);
    Nt = 0;
    for a in argv:
        v = gda_cvec1(a);
        N,M = np.shape(v);
        w[Nt:Nt+N,0:1] = v;  # patch 20230418 was #w[Nt:Nt+N,0] = v[0:N,0];
        Nt = Nt + N;
    return w;

# function to make a numpy N-by-1 column vector
# c=gda_cvec1(v) from entity v that is array-like,
# including a number, a list of numbers, a tuple
# of numbers, an N-by-0 np array and a N-by1 np array.
# It promotes integers to floats, and integers
# and floats to complex, by context. The function
# also insures that, if v is an np array, that
# c is a copy, as contrasted to a view, of it.
# This version takes just one input argment.
# whereas c=gda_cvec(v1,v2,...) concatenates
# many argiments.
def gda_cvec1(v):
    if isinstance(v, int) or isinstance(v, np.int32):
        w = np.zeros((1,1),dtype=int);
        w[0,0] = v;
        return w;
    elif isinstance(v, float):
        w = np.zeros((1,1),dtype=float);
        w[0,0] = v;
        return w;
    elif isinstance(v, complex):
        w = np.zeros((1,1),dtype=complex);
        w[0,0] = v;
        return w;
    elif isinstance(v, np.ndarray):
        s = np.shape(v);
        if len(s) == 1:
            return np.copy(np.reshape(v,(s[0],1)));
        else:
            [r,c]=s;
            if( c==1 ):
                return(np.copy(v));
            elif(r==1):
                return(np.copy(v.T));
            else:
                raise TypeError("gda_cvec: %d by %d ndarray not allowed" % (r, c));
    elif isinstance(v, list):
        r = len(v);
        t = int;
        for vi in v:
            if isinstance(vi,int) or isinstance(vi, np.int32): #patch v->vi 20230418
                pass;
            elif isinstance(vi,float):
                t=float;
            elif isinstance(vi,complex):
                t=complex;
                break;
            else:
                raise TypeError("gda_cvec: list contains unsupported type %s" % type(vi));
        w = np.zeros((r,1),dtype=t);
        w[:,0] = np.array(v); # patch v -> np.array(v)
        return w;
    elif isinstance(v, tuple):
        r = len(v);
        t = int;
        for vi in v:
            if isinstance(vi,int) or isinstance(vi, np.int32): #patch v->vi 20230418
                pass;
            elif isinstance(vi,float):
                t=float;
            elif isinstance(vi,complex):
                t=complex;
                break;
            else:
                raise TypeError("gda_cvec: tuple contains unsupported type %s" % type(vi));
        w = np.zeros((r,1),dtype=t);
        w[:,0] = np.array(list(v)); # patch v -> np.array(list(v));
        return w;
    else:
        raise TypeError("gda_cvec: %s not supported" % type(v));

        
# gda_draw function makes a "pictorial matrix equation"
# arguments are vectors, matrices and strings
# which are plotted in the order that the appear
# except that strings starting with 'title ' are plotted
# under the subseqeunt matrix or vector
# always returns a status of 1
def gda_draw(*argv):
    DOCOLOR=True;
    if( DOCOLOR ):
        bwcmap = matplotlib.colormaps['jet'];
    else:
        bw = np.zeros((256,4));
        v = 0.9*(256 - np.linspace( 0, 255, 256 ))/255;
        bw[:,0] = v;
        bw[:,1] = v;
        bw[:,2] = v;
        bw[:,3] = np.ones(256);
        bwcmap = ListedColormap(bw);
    # size of plot
    W = 16;
    H = 4;
    fig1 = plt.figure(1);
    # figsize width and height in inches
    fig1.set_size_inches(W,H);
    ax1 = plt.subplot(1,1,1);
    plt.axis([0, W, -H/2, H/2]);
    plt.axis('off');
    LM = W/6;    # matrix width and heoght
    LV = W/40;   # vector width
    FS = 0.12;    # character width
    TO = 0.4;    # title vertical offset
    SP = 0.2;    # space between objects
    LS = 0.2;    # leading space
    p = LS; # starting x-position
    istitle=0; # flags presence of a title
    for a in argv:
        if isinstance(a,np.ndarray):
            sh = np.shape(a);
            if len(sh) == 1:  # conversion to nx1 array
                n = sh[0];
                m = 1;
                ap = a;
                a = np.zeros((n,1));
                a[:,0] = ap;
            else:
                n = sh[0];
                m = sh[1];
            if m==1:
                pold=p;
                left=p;
                right=p+LV;
                bottom=-LM/2;
                top=LM/2;
                plt.imshow( a, cmap=bwcmap, vmin=np.min(a), vmax=np.max(a), extent=(left,right,bottom,top) );
                p = p+LV;
                pm = (p+pold)/2;
                if istitle:
                    plt.text(pm,-(LM/2)-TO,titlestr,horizontalalignment='center');
                    istitle=0;
                p = p+SP;
            else:
                pold=p;
                left=p;
                right=p+LM;
                bottom=-LM/2;
                top=LM/2;
                plt.imshow( a, cmap=bwcmap, vmin=np.min(a), vmax=np.max(a), extent=(left,right,bottom,top) );
                p = p+LM;
                pm = (p+pold)/2;
                if istitle:
                    plt.text(pm,-(LM/2)-TO,titlestr,horizontalalignment='center');
                    istitle=0;
                p = p+SP;
        elif isinstance(a,str):
            ns = len(a);
            istitle=0;
            if( ns>=6 ):
                if 'title ' in a[0:6]:
                    istitle=1;
                    titlestr=a[6:];
            if( istitle != 1):
                plt.text(p,0,a);
                p = p + ns*FS + SP;
    plt.show();
    return 1;

# bandpass filter, used in seismological example, but hand
# in a variety of settings involving time series
def gda_chebyshevfilt(d, Dt, flow, fhigh):
    # (dout,u,v)=gda_chebyshevfilt(d, Dt, flow, fhigh);
    # chebyshev IIR bandpass filter
    # d - input array of data
    # Dt - sampling interval
    # flow - low pass frequency, Hz
    # fhigh - high pass frequency, Hz
    # dout - output array of data
    # u - the numerator filter
    # v - the denominator filter
    # these filters can be used again using dout=filter(u,v,din);

    # make sure input timeseries is a column vector
    s = np.shape(d);
    N = s[0];
    if(N==1):
        dd = np.zeros((N,1));
        dd[:,0] = d;
    else:
        dd=d;
        
    # sampling rate
    rate=1/Dt;

    # ripple parameter, set to ten percent
    ripple=0.1;  

    # normalise frequency
    fl=2.0*flow/rate;
    fh=2.0*fhigh/rate;

    # center frequency 
    cf = 4 * tan( (fl*pi/2) ) * tan( (fh*pi/2) );

    # bandwidth
    bw = 2 * ( tan( (fh*pi/2) ) - tan( (fl*pi/2) ) );

    # ripple parameter factor
    rpf = sqrt((sqrt((1.0+1.0/(ripple*ripple))) + 1.0/ripple));
    a = 0.5*(rpf-1.0/rpf);
    b = 0.5*(rpf+1.0/rpf);

    u=np.zeros((5,1));
    v=np.zeros((5,1));
    theta = 3*pi/4;
    sr = a * cos(theta);
    si = b * sin(theta);
    es = sqrt(sr*sr+si*si);
    tmp= 16.0 - 16.0*bw*sr + 8.0*cf + 4.0*es*es*bw*bw - 4.0*bw*cf*sr + cf*cf;
    v[0,0] = 1.0;
    v[1,0] = 4.0*(-16.0 + 8.0*bw*sr - 2.0*bw*cf*sr + cf*cf)/tmp;
    v[2,0] = (96.0 - 16.0*cf - 8.0*es*es*bw*bw + 6.0*cf*cf)/tmp;
    v[3,0] = (-64.0 - 32.0*bw*sr + 8.0*bw*cf*sr + 4.0*cf*cf)/tmp;
    v[4,0] = (16.0 + 16.0*bw*sr + 8.0*cf + 4.0*es*es*bw*bw + 4.0*bw*cf*sr + cf*cf)/tmp;
    tmp = 4.0*es*es*bw*bw/tmp;
    u[0,0] = tmp;
    u[1,0] = 0.0;
    u[2,0] = -2.0*tmp;
    u[3,0] = 0.0;
    u[4,0] = tmp;

    dout = sg.lfilter(u.ravel(),v.ravel(),dd.ravel());
    return (gda_cvec(dout),gda_cvec(u),gda_cvec(v));

# function to perform the multiplication FT (F v)
def gdaFTFmul(v):
    # this function is used by bicongugate gradient to
    # solve the generalized least squares problem Fm=F.
    # Note that "F" must be called "gdaFsparse".
    global gdaFsparse;
    s = np.shape(v);
    if(len(s)==1):
        vv = np.zeros((s[0],1));
        vv[:,0] = v;
    else:
        vv=v;
    # weirdly, "*" multiplies sparse matrices just fine
    temp = gdaFsparse*vv;
    return gdaFsparse.transpose()*temp;

# function to perform the multiplication FT (F v)
def gda_FTFmul(v):
    # this function is used by bicongugate gradient to
    # solve the generalized least squares problem Fm=F.
    # Note that "F" must be called "gdaFsparse".
    global gdaFsparse;
    s = np.shape(v);
    if(len(s)==1):
        vv = np.zeros((s[0],1));
        vv[:,0] = v;
    else:
        vv=v;
    # weirdly, "*" multiplies sparse matrices just fine
    temp = gdaFsparse*vv;
    return gdaFsparse.transpose()*temp;


# function to perform the multiplication FT f
def gda_FTFrhs(f):
    # this function is used by bicongugate gradient to
    # solve the generalized least squares problem Fm=F.
    # Note that "F" must be called "gdaFsparse".
    global gdaFsparse;
    return gdaFsparse.transpose()*f;

def gda_matrixmin(E):
    E1 = np.min(E,axis=0);
    k = np.argmin(E,axis=0);
    Emin = np.min(E1);
    ic = np.argmin(E1);
    ir = k[ic];
    return ir, ic, Emin;
In [2]:
# gdapy06_01
# data point inside of 3D box with 1:1:1 refernce line

# this version differs from the MATLAB(R) version in that it uses
# a wirefram grid to plot a sphere, as contrasted to plotting a
# shaded patch.

print("gdapy06_01");

# x-axis
Nx = 31;
xmin = 0.0;
xmax = 5.0;
Dx = (xmax-xmin)/(Nx-1);
x = gda_cvec(xmin + Dx*np.linspace(0,Nx-1,Nx));
 
# y-axis
Ny = 31;
ymin = 0.0;
ymax = 5.0;
Dy = (ymax-ymin)/(Ny-1);
y = gda_cvec(ymin + Dy*np.linspace(0,Ny-1,Ny));
 
# z-axis
Nz = 31;
zmin = 0;
zmax = 5;
Dz = (zmax-zmin)/(Nz-1);
z = gda_cvec(zmin + Dz*np.linspace(0,Nz-1,Nz));
         
# the point is at a randomly chosen point in d-space
# using a Normal distribution with specified mean and variance
rbar = np.random.normal(loc=2.5,scale=0.5,size=(3,1));

# setup for 3D figure
fig1 = plt.figure(1,figsize=(12,5));   # smallish figure
ax = fig1.add_subplot(111, projection='3d');
plt.axis('off');
plt.grid(b=None);
ax.axes.set_xlim3d(left=xmin, right=xmax);
ax.axes.set_ylim3d(bottom=ymin, top=ymax);
ax.axes.set_zlim3d(bottom=zmin, top=zmax);
 
# improvise 3D box
ax.plot3D( gda_cvec(xmin,xmin).ravel(), gda_cvec(ymin,ymin).ravel(), gda_cvec(zmin,zmax).ravel(), "k-", lw=2 );
ax.plot3D( gda_cvec(xmin,xmin).ravel(), gda_cvec(ymin,ymax).ravel(), gda_cvec(zmin,zmin).ravel(), "k-", lw=2 );
ax.plot3D( gda_cvec(xmin,xmax).ravel(), gda_cvec(ymin,ymin).ravel(), gda_cvec(zmin,zmin).ravel(), "k-", lw=2 );
ax.plot3D( gda_cvec(xmax,xmax).ravel(), gda_cvec(ymax,ymax).ravel(), gda_cvec(zmax,zmin).ravel(), "k-", lw=2 );
ax.plot3D( gda_cvec(xmax,xmax).ravel(), gda_cvec(ymax,ymin).ravel(), gda_cvec(zmax,zmax).ravel(), "k-", lw=2 );
ax.plot3D( gda_cvec(xmax,xmin).ravel(), gda_cvec(ymax,ymax).ravel(), gda_cvec(zmax,zmax).ravel(), "k-", lw=2 );
ax.plot3D( gda_cvec(xmax,xmin).ravel(), gda_cvec(ymin,ymin).ravel(), gda_cvec(zmax,zmax).ravel(), "k-", lw=2 );
ax.plot3D( gda_cvec(xmax,xmax).ravel(), gda_cvec(ymin,ymin).ravel(), gda_cvec(zmax,zmin).ravel(), "k-", lw=2 );
ax.plot3D( gda_cvec(xmin,xmin).ravel(), gda_cvec(ymax,ymin).ravel(), gda_cvec(zmax,zmax).ravel(), "k-", lw=2 );
ax.plot3D( gda_cvec(xmin,xmin).ravel(), gda_cvec(ymax,ymax).ravel(), gda_cvec(zmax,zmin).ravel(), "k-", lw=2 );
ax.plot3D( gda_cvec(xmax,xmax).ravel(), gda_cvec(ymax,ymin).ravel(), gda_cvec(zmin,zmin).ravel(), "k-", lw=2 );
ax.plot3D( gda_cvec(xmax,xmin).ravel(), gda_cvec(ymax,ymax).ravel(), gda_cvec(zmin,zmin).ravel(), "k-", lw=2 );
 
# plot 1:1:1 line
ax.plot3D( gda_cvec(0.0,5.0).ravel(), gda_cvec(0.0,5.0).ravel(), gda_cvec(0.0,5.0).ravel(), "b-", lw=3 );

# draw sphere using wireframe technique
Nu = 21;
u0 = 2.0*pi*np.linspace(0,Nu-1,Nu)/(Nu-1);
Nv = 11;
v0 = 2.0*pi*np.linspace(0,Nv-1,Nv)/(Nv-1);
u, v = np.meshgrid(u0,v0);
diameter=0.2;
x = diameter*np.cos(u)*np.sin(v)+rbar[0,0];
y = diameter*np.sin(u)*np.sin(v)+rbar[1,0];
z = diameter*np.cos(v)+rbar[2,0];
ax.plot_wireframe(x, y, z, color="r")
plt.show();

gdapy06_01
No description has been provided for this image
In [3]:
# gdapy06_02
# 3D Normal pdf plotted along along 1:1:1 line

# this version differs from the MATLAB(R) version in that it uses
# a dots to visualize the pdf, whereas MATLAB(R) uses multiple surface

print('gdapy06_02\n');

# x-axis
Nx = 31;
xmin = 0.0;
xmax = 5.0;
Dx = (xmax-xmin)/(Nx-1);
x = gda_cvec(xmin + Dx*np.linspace(0,Nx-1,Nx));
 
# y-axis
Ny = 31;
ymin = 0.0;
ymax = 5.0;
Dy = (ymax-ymin)/(Ny-1);
y = gda_cvec(ymin + Dy*np.linspace(0,Ny-1,Ny));
 
# z-axis
Nz = 31;
zmin = 0;
zmax = 5;
Dz = (zmax-zmin)/(Nz-1);
z = gda_cvec(zmin + Dz*np.linspace(0,Nz-1,Nz));

# setup for 3D figure
fig1 = plt.figure(1,figsize=(12,5));   # smallish figure
ax = fig1.add_subplot(111, projection='3d');
plt.axis('off');
plt.grid(b=None);
ax.axes.set_xlim3d(left=xmin, right=xmax);
ax.axes.set_ylim3d(bottom=ymin, top=ymax);
ax.axes.set_zlim3d(bottom=zmin, top=zmax);
 
# improvise 3D box
ax.plot3D( gda_cvec(xmin,xmin).ravel(), gda_cvec(ymin,ymin).ravel(), gda_cvec(zmin,zmax).ravel(), "k-", lw=2 );
ax.plot3D( gda_cvec(xmin,xmin).ravel(), gda_cvec(ymin,ymax).ravel(), gda_cvec(zmin,zmin).ravel(), "k-", lw=2 );
ax.plot3D( gda_cvec(xmin,xmax).ravel(), gda_cvec(ymin,ymin).ravel(), gda_cvec(zmin,zmin).ravel(), "k-", lw=2 );
ax.plot3D( gda_cvec(xmax,xmax).ravel(), gda_cvec(ymax,ymax).ravel(), gda_cvec(zmax,zmin).ravel(), "k-", lw=2 );
ax.plot3D( gda_cvec(xmax,xmax).ravel(), gda_cvec(ymax,ymin).ravel(), gda_cvec(zmax,zmax).ravel(), "k-", lw=2 );
ax.plot3D( gda_cvec(xmax,xmin).ravel(), gda_cvec(ymax,ymax).ravel(), gda_cvec(zmax,zmax).ravel(), "k-", lw=2 );
ax.plot3D( gda_cvec(xmax,xmin).ravel(), gda_cvec(ymin,ymin).ravel(), gda_cvec(zmax,zmax).ravel(), "k-", lw=2 );
ax.plot3D( gda_cvec(xmax,xmax).ravel(), gda_cvec(ymin,ymin).ravel(), gda_cvec(zmax,zmin).ravel(), "k-", lw=2 );
ax.plot3D( gda_cvec(xmin,xmin).ravel(), gda_cvec(ymax,ymin).ravel(), gda_cvec(zmax,zmax).ravel(), "k-", lw=2 );
ax.plot3D( gda_cvec(xmin,xmin).ravel(), gda_cvec(ymax,ymax).ravel(), gda_cvec(zmax,zmin).ravel(), "k-", lw=2 );
ax.plot3D( gda_cvec(xmax,xmax).ravel(), gda_cvec(ymax,ymin).ravel(), gda_cvec(zmin,zmin).ravel(), "k-", lw=2 );
ax.plot3D( gda_cvec(xmax,xmin).ravel(), gda_cvec(ymax,ymax).ravel(), gda_cvec(zmin,zmin).ravel(), "k-", lw=2 );
 
# plot 1:1:1 line
ax.plot3D( gda_cvec(0.0,5.0).ravel(), gda_cvec(0.0,5.0).ravel(), gda_cvec(0.0,5.0).ravel(), "b-", lw=3 );

# Normal pdf parameter
rbar = gda_cvec(2.5, 2.5, 2.5);  # center
sd=0.5; # sqrt(variance);

# render normal pdf as a bunch of dots in 3D space
Nr = 10;
r = gda_cvec(0.05, 0.1, 0.15, 0.2, 0.3, 0.3, 0.5, 0.8); # radius of dots
Nth = 10;
Dth = 2*pi/Nth;
Nph = 5;
Dph = pi/Nph;
for ith in range(Nth):  # loop over horizontal angle
    th = Dth*ith;
    sth = sin(th);
    cth = cos(th);
    for iph in range(Nph): # loop over vertical angle
        ph = Dph*iph;
        sph = sin(ph);
        cph = cos(ph);
        xr = rbar[0,0]+sin(th)*sin(ph)*r;
        yr = rbar[1,0]+cos(th)*sin(ph)*r;
        zr = rbar[2,0]+cos(ph)*r;
        ax.scatter(xr, yr, zr, marker=".", color="r");
plt.show();

gdapy06_02
No description has been provided for this image
In [4]:
# gdapy06_03
# Likelihood surface for mean/variance
 
print("gdapy06_03");

# random data
N = 100;
meantrue=2.5;
sigmatrue=1.5;
d = np.random.normal(loc=meantrue,scale=sigmatrue,size=(N,1));
 
# m1-axis, with m1=mean
Nm = 51;
mmin = 1;
mmax = 5.0;
Dm = (mmax-mmin)/(Nm-1);
m = gda_cvec(mmin + Dm*np.linspace(0,Nm-1,Nm));
 
#% s-axis, with s=sqrt(variance)
Ns = 51;
smin = 1;
smax = 5.0;
Ds = (smax-smin)/(Ns-1);
s = gda_cvec(smin + Ds*np.linspace(0,Ns-1,Ns));
 
# (m,s) grid
[X,Y]=np.meshgrid( m, s);
 
# likelihood surface
L=np.zeros((Nm,Ns));
 
# tabulate likelihood surface
# Normal P = (1/sqrt(2pi)) * (1/s) * exp( -0.5 (d-m)^2 / s^2 )
#        L = -0.5log(2*pi) - log(s) - ( -0.5 (d-m)^2 / s^2 )
for i in range(Nm):
    for j in range(Ns):
        mp=X[i,j];
        sp=Y[i,j];
        L[i,j] = -N*log(2*pi)/2 - N*log(sp) - 0.5*np.sum(np.power(d-mp,2))/(sp**2);

# find (i,j) of point of maximum
Li, Lj, Lmax = gda_matrixmin(-L);
Lmax=-Lmax;

# minimum for plotting purposes
i, j, Lmin = gda_matrixmin(L);

# maximum likelihood point
mbest = X[Li,Lj];
sbest = Y[Li,Lj];
print("True: mean %f sigma %f" % (meantrue,sigmatrue) );
print("Est:  mean %f sigma %f Likelihood %f" % (mbest,sbest,Lmax) );

fig1 = plt.figure(1,figsize=(12,5));   # smallish figure
ax1 = fig1.add_subplot(111, projection='3d');
mycmap=matplotlib.colormaps['jet'];
surf = ax1.plot_surface(X, Y, L, cmap=mycmap);
fig1.colorbar(surf, shrink=0.5);
Dm=0.4;
Ds=0.4;
DL=50.0;
ax1.plot3D( gda_cvec(mbest,mbest).ravel(), gda_cvec(sbest,sbest).ravel(), gda_cvec(Lmin,Lmax).ravel(), "k-", lw=3 );
ax1.plot3D( gda_cvec(mbest,mbest).ravel(), gda_cvec(sbest-Ds,sbest+Ds).ravel(), gda_cvec(Lmin,Lmin).ravel(), "k-", lw=3 );
ax1.plot3D( gda_cvec(mbest-Dm,mbest+Dm).ravel(), gda_cvec(sbest,sbest).ravel(), gda_cvec(Lmin,Lmin).ravel(), "k-", lw=3 );
ax1.plot3D( gda_cvec(meantrue,meantrue).ravel(), gda_cvec(sigmatrue-Ds,sigmatrue+Ds).ravel(), gda_cvec(Lmin,Lmin).ravel(), "r-", lw=1 );
ax1.plot3D( gda_cvec(meantrue-Dm,meantrue+Dm).ravel(), gda_cvec(sigmatrue,sigmatrue).ravel(), gda_cvec(Lmin,Lmin).ravel(), "r-", lw=1 );


plt.xlabel("mean");
plt.ylabel("sigma");
plt.show();

print("Caption: Depiction of likelihood surface L(mean,sigma) (colors)");
print("with estimated value (black) and true value (red) superimposed.");

gdapy06_03
True: mean 2.500000 sigma 1.500000
Est:  mean 2.280000 sigma 1.400000 Likelihood -175.567400
No description has been provided for this image 
Caption: Depiction of likelihood surface L(mean,sigma) (colors)
with estimated value (black) and true value (red) superimposed.
In [5]:
# gdapy06_04
# examples of probability density function p(d1,d2)
# one is peaked, the other had a ridge
 
print("gdapy06_04");

# d1-axis
Nd1 = 51;
d1min = 0.0;
d1max = 5.0;
Dd1 = (d1max-d1min)/(Nd1-1);
d1 = gda_cvec(d1min + Dd1*np.linspace(0,Nd1-1,Nd1));
 
# d2-axis
Nd2 = 51;
d2min = 0.0;
d2max = 5.0;
Dd2 = (d2max-d2min)/(Nd2-1);
d2 = gda_cvec(d2min + Dd2*np.linspace(0,Nd2-1,Nd2));
 
# (d1,d2) grid
X, Y = np.meshgrid( d1, d2);
 
# distribution 1, has a peak
P1=np.zeros((Nd1,Nd2));
dbar = gda_cvec(2.5, 2.5);
sd1 = 0.5;
sd2 = 0.5;
C = np.diag( gda_cvec(sd1, sd2).ravel() );
CI = la.inv(C);
DC = la.det(C);
norm = (1.0/(2.0*pi)) * (1.0/sqrt(DC));
for i in range(Nd1):
    for j in range(Nd2):
        d = gda_cvec( X[i,j], Y[i,j] ) - dbar;
        r2 = np.matmul(d.T,np.matmul(CI,d)); r2=r2[0,0];
        P1[i,j] = norm*exp( -0.5 * r2 );

# for test purposes
# A = Dd1*Dd2*sum(sum(P1)
 
# distribution 2, ridge
# Note distribution not normalizable
P2=np.zeros((Nd1,Nd2));
sda = 0.5;
d0a = gda_cvec(2.5, 2.5);
w = np.ones((2,1));
for i in range(Nd1):
    for j in range(Nd2):
        dr =gda_cvec( X[i,j]-d0a[0,0], Y[i,j]-d0a[1,0] );
        r2= np.matmul(dr.T,w)**2; r2=r2[0,0];
        P2[i,j] = 0.2*exp(-0.5*r2/sda);

# plot
fig1 = plt.figure(1,figsize=(12,5));   # smallish figure
ax1 = fig1.add_subplot(121, projection='3d');
mycmap=matplotlib.colormaps['jet'];
surf = ax1.plot_surface(X, Y, P1, cmap=mycmap);
plt.xlabel("x");
plt.ylabel("y");

ax2 = fig1.add_subplot(122, projection='3d');
mycmap=matplotlib.colormaps['jet'];
surf = ax2.plot_surface(X, Y, P2, cmap=mycmap);
plt.xlabel("x");
plt.ylabel("y");
plt.show();

print("Caption: (left) Pdf p(x1,x2) with a peak. (right) Pdf");
print("p(x1,x2) with a ridge.");

gdapy06_04
No description has been provided for this image 
Caption: (left) Pdf p(x1,x2) with a peak. (right) Pdf
p(x1,x2) with a ridge.
In [6]:
# gdapy06_05
# examples of probability density function p(m1,m2)
# uncorrelated, equal variance, one with small
# variance, the other with large variance
 
print("gdapy06_05");

# m1-axis
Nm1 = 51;
m1min = 0.0;
m1max = 5.0;
Dm1 = (m1max-m1min)/(Nm1-1);
m1 = gda_cvec(m1min + Dm1*np.linspace(0,Nm1-1,Nm1));
 
# m2-axis
Nm2 = 51;
m2min = 0;
m2max = 5.0;
Dm2 = (m2max-m2min)/(Nm2-1);
m2 = gda_cvec(m2min + Dm2*np.linspace(0,Nm2-1,Nm2));
 
# setup for distribution 1, small variance
P1=np.zeros((Nm1,Nm2));
mbar1 = gda_cvec(2.5, 2.5);
sd1 = 0.5;
C1 = np.diag( gda_cvec(sd1**2, sd1**2).ravel() );
CI1 = la.inv(C1);
DC1 = la.det(C1);
 
# setup for distributions 2, large variance
P2=np.zeros((Nm1,Nm2));
mbar2 = gda_cvec(2.5, 2.5);
sd2 = 1.0;
C2 = np.diag( gda_cvec( sd2**2, sd2**2).ravel() );
CI2 = la.inv(C2);
DC2 = la.det(C2);
 
# make distributions
norm1 = (1/(2*pi)) * (1/sqrt(DC1));
norm2 = (1/(2*pi)) * (1/sqrt(DC2));
for i in range(Nm1):
    for j in range(Nm2):
        x1 =gda_cvec(m1[i,0], m2[j,0]) - mbar1;
        x2 =gda_cvec(m1[i,0], m2[j,0]) - mbar2;
        r2 = np.matmul(x1.T,np.matmul(CI1,x1)); r2=r2[0,0];
        P1[i,j] = norm1*exp( -0.5 * r2 );
        r2 = np.matmul(x1.T,np.matmul(CI2,x1)); r2=r2[0,0];
        P2[i,j] = norm2*exp( -0.5 * r2 );

# for test purposes
# A1 = Dd1*Dd2*sum(sum(P1));
# A2 = Dd1*Dd2*sum(sum(P2));
 
# draw distributions
gda_draw(P1,"   ",P2 );
print("Caption: (Left) Pdf p(m1,m2) with small and equal sigmas. (right)");
print("(Right) Pdf p(m1,m2) with large and equal sigmas. (right)");

gdapy06_05
No description has been provided for this image 
Caption: (Left) Pdf p(m1,m2) with small and equal sigmas. (right)
(Right) Pdf p(m1,m2) with large and equal sigmas. (right)
In [7]:
# gdapy06_06

# examples of probability density function p(m1,m2)
# unequal variance, uncorrelated

print("gdapy06_06");

# m1-axis
Nm1 = 51;
m1min = 0.0;
m1max = 5.0;
Dm1 = (m1max-m1min)/(Nm1-1);
m1 = gda_cvec( m1min + Dm1*np.linspace(0.0,Nm1-1,Nm1));
 
# m2-axis
Nm2 = 51;
m2min = 0.0;
m2max = 5.0;
Dm2 = (m2max-m2min)/(Nm2-1);
m2 = gda_cvec( m2min + Dm2*np.linspace(0.0,Nm2-1,Nm2));
 
# setup for distribution 1, C11<C22
P1=np.zeros((Nm1,Nm2));
mbar1 = gda_cvec(2.5, 2.5);
sd1 = 0.5;
sd2 = 1.0;
C1 = np.diag( gda_cvec(sd1**2, sd2**2).ravel() );
CI1 = la.inv(C1);
DC1 = la.det(C1);
 
# normalization
norm1 = (1.0/(2.0*pi)) * (1.0/sqrt(DC1));
 
# build distribution
for i in range(Nm1):
    for j in range(Nm2):
        x1 = gda_cvec( m1[i,0], m2[j,0] ) - mbar1;
        r2 = np.matmul(x1.T,np.matmul(CI1,x1)); r2=r2[0,0];
        P1[i,j] = norm1*exp( -0.5 * r2 );
 
# for test purposes
# A1 = Dd1*Dd2*np.sum(P1);
 
gda_draw(P1);
print("Caption: Pdf p(m1,m2) with unequal sigmas and zero covarince.");

gdapy06_06
No description has been provided for this image 
Caption: Pdf p(m1,m2) with unequal sigmas and zero covarince.
In [8]:
# gdapy06_07

# examples of probability density function p(m1,m2)
# one is a ridge, the other is Normal but mimics a ridge

print("gdapy06_07");

# m1-axis
Nm1 = 51;
m1min = 0.0;
m1max = 5.0;
Dm1 = (m1max-m1min)/(Nm1-1);
m1 = gda_cvec(m1min + Dm1*np.linspace(0,Nm1-1,Nm1));
 
# m2-axis
Nm2 = 51;
m2min = 0.0;
m2max = 5.0;
Dm2 = (m2max-m2min)/(Nm2-1);
m2 = gda_cvec(m2min + Dm2*np.linspace(0,Nm2-1,Nm2));
 
# setup for distribution 1, Normal
P1=np.zeros((Nm1,Nm2));
mbar1 = gda_cvec(2.5, 2.5);
sd12 = 2.0;
cov = 2.0-0.1;
C1 = np.array( [ [sd12, cov], [cov, sd12] ] );
CI1 = la.inv(C1);
DC1 = la.det(C1);
norm1 = (1.0/(2.0*pi)) * (1.0/sqrt(DC1));
 
# calculate pdf
for i in range(Nm1):
    for j in range(Nm2):
        x1 = gda_cvec( m1[i,0], m2[j,0] ) - mbar1;
        r2 = np.matmul(x1.T,np.matmul(CI1,x1)); r2=r2[0,0];
        P1[i,j] = norm1*exp( -0.5 * r2 );

# setup for distribution 2, ridge
# Note distribution not normalizable
P2=np.zeros((Nm1,Nm2));
sda2 = 0.5;
d0a = gda_cvec(2.5, 2.5);
w = gda_cvec(1.0,-1.0);
for i in range(Nm1):
    for j in range(Nm2):
        dr = gda_cvec( m1[i,0], m2[j,0] ) - d0a;
        r2 = np.matmul(dr.T,w)**2; r2=r2[0,0];
        P2[i,j] = 0.2*exp(-0.5*r2/sda2);

# plot distributions
gda_draw(P2, "   ", P1 );

gdapy06_07
No description has been provided for this image
In [9]:
# gdapy06_08

# example of probability density function p(m1,m2)
# implements inequality constraint m1<m2
 
print("gdapy06_08");

# m1-axis
Nm1 = 51;
m1min = 0.0;
m1max = 5.0;
Dm1 = (m1max-m1min)/(Nm1-1);
m1 = gda_cvec(m1min + Dm1*np.linspace(0,Nm1-1,Nm1));
 
# m2-axis
Nm2 = 51;
m2min = 0.0;
m2max = 5.0;
Dm2 = (m2max-m2min)/(Nm2-1);
m2 = gda_cvec(m2min + Dm2*np.linspace(0,Nm2-1,Nm2));
 
# compute pdf
P1=np.zeros((Nm1,Nm2));
for i in range(Nm1):
    for j in range(Nm2):
        P1[i,j] = float(m1[i,0]<=m2[j,0]);

# plot distribution
gda_draw(P1 );
print("Caption: Pdf p(m1,m2) implementing the inequality constraint m1<m2.");

gdapy06_08
No description has been provided for this image 
Caption: Pdf p(m1,m2) implementing the inequality constraint m1<m2.
In [10]:
# gdapy06_09
 
# relative entropy of onenormsl pdf p(m)
# with respect to another q(m)
 
print("gdapy06_09");

# m-axis
N=101;
mmin = -25;
mmax = 25;
Dm = (mmax-mmin)/(N-1);
m = gda_cvec(mmin + Dm*np.linspace(0,N-1,N));
 
# p(m) = pA(m)
mbarp = 0.0;
sigmamp = 1;
p = (1.0/(sqrt(2.0*pi)*sigmamp))*np.exp(-(np.power(m-mbarp,2)/(2.0*sigmamp*sigmamp)));
normp = Dm*np.sum(p);
 
# q(m) = pN(m)
mbarq = 0.0;
sigmamq = 5.0;
q = (1.0/(sqrt(2.0*pi)*sigmamq))*np.exp(-np.power(m-mbarq,2)/(2.0*sigmamq*sigmamq));
normq = Dm*np.sum(q);
 
# plot
fig1 = plt.figure(1,figsize=(12,5));
ax1 = plt.subplot(1,1,1);
plt.axis( [mmin, mmax, 0, max(p)] );
plt.plot( m, p, "r-", lw=3);
plt.plot( m, q, "g-", lw=3);
plt.xlabel("m");
plt.ylabel("p and q");
plt.show();
print("Caption: Two Normal pdfs p(m) (red) and q(m) (green)");

# relative entropy S
r = (sigmamp**2)/(sigmamq**2);
Sanalytic = ((mbarp-mbarq)**2)/(2*sigmamq*sigmamq) + 0.5*(r-1.0-log(r));
Snumeric = Dm*np.sum(np.multiply(p,np.log(np.divide(p,q))));
print("S = %f (analytic) and %f (numeric)" % (Sanalytic, Snumeric) );

# relative entropy as a function of sigmamp
Nv=100;
sigmampv = gda_cvec(sigmamq*np.linspace(1.0,Nv,Nv)/Nv);
rv = np.power(sigmampv,2)/(sigmamq**2);
Sv = ((mbarp-mbarq)**2)/(2*sigmamq*sigmamq) + 0.5*(rv-1-np.log(rv));

fig1 = plt.figure(1,figsize=(12,5));
ax1 = plt.subplot(1,1,1);
plt.axis( [0, sigmamq, 0, np.max(Sv)] );
plt.plot( sigmampv, Sv, "k-", lw=3);
plt.plot( sigmamp, Sanalytic, "ro", lw=3);
plt.xlabel('sigmamp');
plt.ylabel('S')
plt.show();
print("Caption: Relative entropy of p(m) and q(m) as the");
print("sigma of p id varied. One point on curve (red circle)");
print("determined by an analytic method");

gdapy06_09
No description has been provided for this image 
Caption: Two Normal pdfs p(m) (red) and q(m) (green)
S = 1.129438 (analytic) and 1.129438 (numeric)
No description has been provided for this image 
Caption: Relative entropy of p(m) and q(m) as the
sigma of p id varied. One point on curve (red circle)
determined by an analytic method
In [11]:
# gdapy06_10

# examples of probability density function p(d1,m1)
print("gdapy06_10");

# d1-axis
Nd1 = 51;
d1min = 0.0;
d1max = 5.0;
Dd1 = (d1max-d1min)/(Nd1-1);
d1 = gda_cvec(d1min + Dd1*np.linspace(0,Nd1-1,Nd1));
 
# m1-axis
Nm1 = 51;
m1min = 0.0;
m1max = 5.0;
Dm1 = (m1max-m1min)/(Nm1-1);
m1 = gda_cvec(m1min + Dm1*np.linspace(0,Nm1-1,Nm1));
 
# setup for pdf
P1=np.zeros((Nd1,Nm1));
mbar1 = gda_cvec(2.5, 2.5);
sd1 = 0.5;
sd2 = 1.0;
C1 = np.diag( gda_cvec(sd1**2, sd2**2).ravel() );
CI1 = la.inv(C1);
DC1 = la.det(C1);
 
# normalization
norm1 = (1.0/(2.0*pi)) * (1.0/sqrt(DC1));
 
# tabulate distribution
for i in range(Nd1):
    for j in range(Nm1):
        x1 = gda_cvec(m1[i,0], m1[j,0]) - mbar1;
        r2 = np.matmul(x1.T,np.matmul(CI1,x1)); r2=r2[0,0];
        P1[i,j] = norm1*exp( -0.5*r2 );

# plot pdf
gda_draw(' ', P1 );
print("Caption: Probability density function p(m,d)");

gdapy06_10
No description has been provided for this image 
Caption: Probability density function p(m,d)
In [12]:
# gdapy06_11

# parameteric function passing thru distribution p(d1,m1)
 
print("gdapy06_11");
 
# d1-axis
Nd1 = 101;
d1min = 0.0;
d1max = 5.0;
Dd1 = (d1max-d1min)/(Nd1-1);
d1 = gda_cvec(d1min + Dd1*np.linspace(0,Nd1-1,Nd1));
 
# m1-axis
Nm1 = 101;
m1min = 0.0;
m1max = 5.0;
Dm1 = (m1max-m1min)/(Nm1-1);
m1 = gda_cvec(m1min + Dm1*np.linspace(0,Nm1-1,Nm1));
 
# setup for pdf
P1=np.zeros((Nd1,Nm1));
d1bar = 2.25;
m1bar = 2.08;
bar = gda_cvec(d1bar, m1bar);
sd1 = 0.5;
sm1 = 1.0;
C1 = np.diag( gda_cvec(sd1**2, sm1**2).ravel() );
CI1 = la.inv(C1);
DC1 = la.det(C1);
norm1 = (1.0/(2.0*pi)) * (1.0/sqrt(DC1));
 
# tablulate pdf
for i in range(Nm1):
    for j in range(Nd1):
        x1 = gda_cvec(d1[i,0], m1[j,0]) - bar;
        r2 = np.matmul(x1.T,np.matmul(CI1,x1)); r2=r2[0,0];
        P1[i,j] = norm1*exp( -0.5*r2 );

# s-axis for parametric curve
Ns = 51;
smin = 0.0;
smax = 5.0;
Ds = (smax-smin)/(Ns-1);
s = gda_cvec(smin + Ds*np.linspace(0,Ns-1,Ns));
 
# parameric curve
dp = 1.0+s-2.0*np.power(s/smax,2);
mp = np.copy(s);
 
# P on parameteric curve
ipd = (np.floor((dp-d1min)/Dd1)).astype(int);
ipm = (np.floor((mp-m1min)/Dm1)).astype(int);
# insure indices in range
k=np.where(ipd<0);
i=k[0];
if( not np.any(i) ):
    ipd[i,0]=0;

k=np.where(ipd>=Nd1);
i=k[0];
if( not np.any(i) ):
    ipd[i,0]=Nd1-1;

k=np.where(ipm<0);
i=k[0];
if( not np.any(i) ):
    ipm[i,0]=0;

k=np.where(ipm>=Nm1);
i=k[0];
if( not np.any(i) ):
    ipm[i,0]=Nm1-1;

Ps=np.zeros((Ns,1));
# evaluate P at indices
for i in range(Ns):
    Ps[i,0] = P1[ipd[i,0], ipm[i,0]];

# maximum along curve
Pmax = np.max(Ps);
ismax = np.argmax(Ps);
d1smax = dp[ismax,0];
m1smax = mp[ismax,0];

fig1 = plt.figure(1,figsize=(12,5));
ax1=plt.subplot(1,1,1);
mycmap=matplotlib.colormaps['jet'];
plt.axis( [d1min, d1max, m1min, m1max] );
left=d1min; right=d1max;
top=m1max; bottom=m1min;
dmax=np.max(P1);
ax1.invert_yaxis();
plt.imshow(np.flipud(P1), cmap=mycmap, vmin=0.0, vmax=dmax, extent=(left,right,bottom,top) );
plt.colorbar();
plt.plot( m1bar, d1bar, "wo", lw=3 );
plt.plot( gda_cvec(m1bar, m1bar), gda_cvec(d1max, d1max-0.1), "w-", lw=3 );
plt.plot( gda_cvec(m1min, m1min+0.1), gda_cvec(d1bar, d1bar), "w-", lw=3 );
plt.plot( mp, dp, 'w-', lw=3 );
plt.plot( m1smax, d1smax, 'ko', lw=4);
plt.plot( gda_cvec(m1min, m1smax), gda_cvec(d1smax, d1smax), "w:", lw=2);
plt.plot( gda_cvec(m1smax, m1smax), gda_cvec(d1max, d1smax), "w:", lw=2);
plt.xlabel("m");
plt.ylabel("d");
plt.show();
print("Caption: Probability density function p(d,m) (colors) with parametric curve");
print("f(d,m)=0 (white) superimposed.  The point of mximum probability");
print("along the curve is shown with a black circle.");


# plot pdf as a function of arclength along curve
fig1 = plt.figure(1,figsize=(12,5));
ax1=plt.subplot(1,1,1);
plt.axis( [smin, smax, 0, 0.3] );
plt.plot( s, Ps, "b-", lw=3 );
plt.xlabel("s");
plt.ylabel("p(s)");
plt.show();
print("Caption: Probability as a function of arc-length along the parametric curve.");
print("The probability has a single maximum.");

gdapy06_11
No description has been provided for this image 
Caption: Probability density function p(d,m) (colors) with parametric curve
f(d,m)=0 (white) superimposed.  The point of mximum probability
along the curve is shown with a black circle.
No description has been provided for this image 
Caption: Probability as a function of arc-length along the parametric curve.
The probability has a single maximum.
In [13]:
# gdapy06_12

# parameteric function passing thru distribution p(d1,m1)
 
print("gdapy06_12");
 
# d1-axis
Nd1 = 101;
d1min = 0.0;
d1max = 5.0;
Dd1 = (d1max-d1min)/(Nd1-1);
d1 = gda_cvec(d1min + Dd1*np.linspace(0,Nd1-1,Nd1));
 
# m1-axis
Nm1 = 101;
m1min = 0.0;
m1max = 5.0;
Dm1 = (m1max-m1min)/(Nm1-1);
m1 = gda_cvec(m1min + Dm1*np.linspace(0,Nm1-1,Nm1));
 
# setup for pdf
P1=np.zeros((Nd1,Nm1));
d1bar = 2.25;
m1bar = 2.08;
bar = gda_cvec(d1bar, m1bar);
sd1 = 1;
sm1 = 0.2;
C1 = np.diag( gda_cvec(sd1**2, sm1**2).ravel() );
CI1 = la.inv(C1);
DC1 = la.det(C1);
norm1 = (1.0/(2.0*pi)) * (1.0/sqrt(DC1));
 
# tablulate pdf
for i in range(Nm1):
    for j in range(Nd1):
        x1 = gda_cvec(d1[i,0], m1[j,0]) - bar;
        r2 = np.matmul(x1.T,np.matmul(CI1,x1)); r2=r2[0,0];
        P1[i,j] = norm1*exp( -0.5*r2 );

# s-axis for parametric curve
Ns = 51;
smin = 0.0;
smax = 5.0;
Ds = (smax-smin)/(Ns-1);
s = gda_cvec(smin + Ds*np.linspace(0,Ns-1,Ns));
 
# parameric curve
dp = 1.0+s-2.0*np.power(s/smax,2);
mp = np.copy(s);
 
# P on parameteric curve
ipd = (np.floor((dp-d1min)/Dd1)).astype(int);
ipm = (np.floor((mp-m1min)/Dm1)).astype(int);
# insure indices in range
k=np.where(ipd<0);
i=k[0];
if( not np.any(i) ):
    ipd[i,0]=0;

k=np.where(ipd>=Nd1);
i=k[0];
if( not np.any(i) ):
    ipd[i,0]=Nd1-1;

k=np.where(ipm<0);
i=k[0];
if( not np.any(i) ):
    ipm[i,0]=0;

k=np.where(ipm>=Nm1);
i=k[0];
if( not np.any(i) ):
    ipm[i,0]=Nm1-1;

Ps=np.zeros((Ns,1));
# evaluate P at indices
for i in range(Ns):
    Ps[i,0] = P1[ipd[i,0], ipm[i,0]];

# maximum along curve
Pmax = np.max(Ps);
ismax = np.argmax(Ps);
d1smax = dp[ismax,0];
m1smax = mp[ismax,0];

fig1 = plt.figure(1,figsize=(12,5));
ax1=plt.subplot(1,1,1);
mycmap=matplotlib.colormaps['jet'];
plt.axis( [d1min, d1max, m1min, m1max] );
left=d1min; right=d1max;
top=m1max; bottom=m1min;
dmax=np.max(P1);
ax1.invert_yaxis();
plt.imshow(np.flipud(P1), cmap=mycmap, vmin=0.0, vmax=dmax, extent=(left,right,bottom,top) );
plt.colorbar();
plt.plot( m1bar, d1bar, "wo", lw=3 );
plt.plot( gda_cvec(m1bar, m1bar), gda_cvec(d1max, d1max-0.1), "w-", lw=3 );
plt.plot( gda_cvec(m1min, m1min+0.1), gda_cvec(d1bar, d1bar), "w-", lw=3 );
plt.plot( mp, dp, 'w-', lw=3 );
plt.plot( m1smax, d1smax, 'ko', lw=4);
plt.plot( gda_cvec(m1min, m1smax), gda_cvec(d1smax, d1smax), "w:", lw=2);
plt.plot( gda_cvec(m1smax, m1smax), gda_cvec(d1max, d1smax), "w:", lw=2);
plt.xlabel("m");
plt.ylabel("d");
plt.show();
print("Caption: Probability density function p(d,m) (colors) with parametric curve");
print("f(d,m)=0 (white) superimposed.  In this case the prior model is very certain");
print("The point of mximum probability along the curve is shown with a black circle.");

# plot pdf as a function of arclength along curve
fig1 = plt.figure(1,figsize=(12,5));
ax1=plt.subplot(1,1,1);
plt.axis( [smin, smax, 0.0, 1.0] );
plt.plot( s, Ps, "b-", lw=3 );
plt.xlabel("s");
plt.ylabel("p(s)");
plt.show();
print("Caption: Probability as a function of arc-length along the parametric curve.");
print("The probability has a single maximum.");

gdapy06_12
No description has been provided for this image 
Caption: Probability density function p(d,m) (colors) with parametric curve
f(d,m)=0 (white) superimposed.  In this case the prior model is very certain
The point of mximum probability along the curve is shown with a black circle.
No description has been provided for this image 
Caption: Probability as a function of arc-length along the parametric curve.
The probability has a single maximum.
In [14]:
# gdapy06_13

# parameteric function passing thru distribution p(d1,m1)
 
print("gdapy06_13");
 
# d1-axis
Nd1 = 101;
d1min = 0.0;
d1max = 5.0;
Dd1 = (d1max-d1min)/(Nd1-1);
d1 = gda_cvec(d1min + Dd1*np.linspace(0,Nd1-1,Nd1));
 
# m1-axis
Nm1 = 101;
m1min = 0.0;
m1max = 5.0;
Dm1 = (m1max-m1min)/(Nm1-1);
m1 = gda_cvec(m1min + Dm1*np.linspace(0,Nm1-1,Nm1));
 
# setup for pdf
P1=np.zeros((Nd1,Nm1));
d1bar = 2.25;
m1bar = 2.08;
bar = gda_cvec(d1bar, m1bar);
sd1 = 0.2;
sm1 = 1.0;
C1 = np.diag( gda_cvec(sd1**2, sm1**2).ravel() );
CI1 = la.inv(C1);
DC1 = la.det(C1);
norm1 = (1.0/(2.0*pi)) * (1.0/sqrt(DC1));
 
# tablulate pdf
for i in range(Nm1):
    for j in range(Nd1):
        x1 = gda_cvec(d1[i,0], m1[j,0]) - bar;
        r2 = np.matmul(x1.T,np.matmul(CI1,x1)); r2=r2[0,0];
        P1[i,j] = norm1*exp( -0.5*r2 );

# s-axis for parametric curve
Ns = 51;
smin = 0.0;
smax = 5.0;
Ds = (smax-smin)/(Ns-1);
s = gda_cvec(smin + Ds*np.linspace(0,Ns-1,Ns));
 
# parameric curve
dp = 1.0+s-2.0*np.power(s/smax,2);
mp = np.copy(s);
 
# P on parameteric curve
ipd = (np.floor((dp-d1min)/Dd1)).astype(int);
ipm = (np.floor((mp-m1min)/Dm1)).astype(int);
# insure indices in range
k=np.where(ipd<0);
i=k[0];
if( not np.any(i) ):
    ipd[i,0]=0;

k=np.where(ipd>=Nd1);
i=k[0];
if( not np.any(i) ):
    ipd[i,0]=Nd1-1;

k=np.where(ipm<0);
i=k[0];
if( not np.any(i) ):
    ipm[i,0]=0;

k=np.where(ipm>=Nm1);
i=k[0];
if( not np.any(i) ):
    ipm[i,0]=Nm1-1;

Ps=np.zeros((Ns,1));
# evaluate P at indices
for i in range(Ns):
    Ps[i,0] = P1[ipd[i,0], ipm[i,0]];

# maximum along curve
Pmax = np.max(Ps);
ismax = np.argmax(Ps);
d1smax = dp[ismax,0];
m1smax = mp[ismax,0];

fig1 = plt.figure(1,figsize=(12,5));
ax1=plt.subplot(1,1,1);
mycmap=matplotlib.colormaps['jet'];
plt.axis( [d1min, d1max, m1min, m1max] );
left=d1min; right=d1max;
top=m1max; bottom=m1min;
dmax=np.max(P1);
ax1.invert_yaxis();
plt.imshow(np.flipud(P1), cmap=mycmap, vmin=0.0, vmax=dmax, extent=(left,right,bottom,top) );
plt.colorbar();
plt.plot( m1bar, d1bar, "wo", lw=3 );
plt.plot( gda_cvec(m1bar, m1bar), gda_cvec(d1max, d1max-0.1), "w-", lw=3 );
plt.plot( gda_cvec(m1min, m1min+0.1), gda_cvec(d1bar, d1bar), "w-", lw=3 );
plt.plot( mp, dp, 'w-', lw=3 );
plt.plot( m1smax, d1smax, 'ko', lw=4);
plt.plot( gda_cvec(m1min, m1smax), gda_cvec(d1smax, d1smax), "w:", lw=2);
plt.plot( gda_cvec(m1smax, m1smax), gda_cvec(d1max, d1smax), "w:", lw=2);
plt.xlabel("m");
plt.ylabel("d");
plt.show();
print("Caption: Probability density function p(d,m) (colors) with parametric curve");
print("f(d,m)=0 (white) superimposed.  In this case the prior model is very certain");
print("The point of mximum probability along the curve is shown with a black circle.");

# plot pdf as a function of arclength along curve
fig1 = plt.figure(1,figsize=(12,5));
ax1=plt.subplot(1,1,1);
plt.axis( [smin, smax, 0.0, 1.0] );
plt.plot( s, Ps, "b-", lw=3 );
plt.xlabel("s");
plt.ylabel("p(s)");
plt.show();
print("Caption: Probability as a function of arc-length along the parametric curve.");
print("The probability has a single maximum.");

gdapy06_13
No description has been provided for this image 
Caption: Probability density function p(d,m) (colors) with parametric curve
f(d,m)=0 (white) superimposed.  In this case the prior model is very certain
The point of mximum probability along the curve is shown with a black circle.
No description has been provided for this image 
Caption: Probability as a function of arc-length along the parametric curve.
The probability has a single maximum.
In [15]:
# gdapy06_14

# prior distribution pA(d1,m1) interacting with inexact
# theory pg(d1,m1) to produce total distribution pT(d1,m1)
 
print("gdapy06_14");

# d1-axis
Nd1 = 101;
d1min = 0.0;
d1max = 5.0;
Dd1 = (d1max-d1min)/(Nd1-1);
d1 = gda_cvec(d1min + Dd1*np.linspace(0,Nd1-1,Nd1));
 
# m1-axis
Nm1 = 101;
m1min = 0;
m1max = 5.0;
Dm1 = (m1max-m1min)/(Nm1-1);
m1 = gda_cvec( m1min + Dm1*np.linspace(0,Nm1-1,Nm1));
 
# setup for distribution P1=PA
P1=np.zeros((Nd1,Nm1));
d1bar = 2.25;
m1bar = 2.08;
bar = gda_cvec(d1bar, m1bar);
sd1 = 0.5;
sm1 = 1.0;
C1 = np.diag( gda_cvec(sd1**2, sm1**2).ravel() );
CI1 = la.inv(C1);
DC1 = la.det(C1);
# note not normaalized, so max is unity
for i in range(Nm1):
    for j in range(Nd1):
        x1 = gda_cvec(d1[i,0], m1[j,0]) - bar;
        r2 = np.matmul(x1.T,np.matmul(CI1,x1)); r2=r2[0,0];
        P1[i,j] = np.exp( -0.5*r2 );

# s-axis for parametric curve s
Ns = 51;
smin = 0.0;
smax = 5.0;
Ds = (smax-smin)/(Ns-1);
s = gda_cvec(smin + Ds*np.linspace(0,Ns-1,Ns));
 
# Pg follows parameric curve
dp = 1.0+s-2.0*np.power(s/smax,2);
mp = np.copy(s);
 
# P2=Pg distribution; not normalizable
P2 = np.zeros((Nd1,Nm1));
sg1 = 0.35;
sg2 = sg1**2;
for i in range(Nm1):
    for j in range(Nd1):
        r2 = np.power(d1[i,0]-dp,2) + np.power(m1[j,0]-mp,2);
        r2min=np.min(r2);
        P2[i,j] = exp(-r2min/sg2);

# P3=PT is product of PA and Pg
P3 = np.multiply(P1,P2);

# find max
Pi, Pj, P3max = gda_matrixmin(P3)
Pmaxm1 = m1min+Dm1*Pj;
Pmaxd1 = d1min+Dd1*Pi;

fig1 = plt.figure(1,figsize=(12,5));

# plot PA
ax1=plt.subplot(1,3,1);
mycmap=matplotlib.colormaps['jet'];
plt.axis( [d1min, d1max, m1min, m1max] );
left=d1min; right=d1max;
top=m1max; bottom=m1min;
dmax=np.max(P1);
ax1.invert_yaxis();
plt.imshow(np.flipud(P1), cmap=mycmap, vmin=0.0, vmax=dmax, extent=(left,right,bottom,top) );
plt.colorbar(shrink=0.5);
plt.plot( m1bar, d1bar, "wo", lw=3 );
plt.plot( gda_cvec(m1bar, m1bar), gda_cvec(d1max, d1max-0.1), "w-", lw=3 );
plt.plot( gda_cvec(m1min, m1min+0.1), gda_cvec(d1bar, d1bar), "w-", lw=3 );
plt.plot( mp, dp, "w:", lw=2 );
plt.xlabel("m");
plt.ylabel("d");

# plot Pg
ax1=plt.subplot(1,3,2);
mycmap=matplotlib.colormaps['jet'];
plt.axis( [d1min, d1max, m1min, m1max] );
left=d1min; right=d1max;
top=m1max; bottom=m1min;
dmax=np.max(P2);
ax1.invert_yaxis();
plt.imshow(np.flipud(P2), cmap=mycmap, vmin=0.0, vmax=dmax, extent=(left,right,bottom,top) );
plt.colorbar(shrink=0.5);
plt.plot( m1bar, d1bar, "wo", lw=3 );
plt.plot( gda_cvec(m1bar, m1bar), gda_cvec(d1max, d1max-0.1), "w-", lw=3 );
plt.plot( gda_cvec(m1min, m1min+0.1), gda_cvec(d1bar, d1bar), "w-", lw=3 );
plt.plot( mp, dp, "w:", lw=2 );
plt.xlabel("m");
plt.ylabel("d");

# plot PT
ax1=plt.subplot(1,3,3);
mycmap=matplotlib.colormaps['jet'];
plt.axis( [d1min, d1max, m1min, m1max] );
left=d1min; right=d1max;
top=m1max; bottom=m1min;
dmax=np.max(P3);
ax1.invert_yaxis();
plt.imshow(np.flipud(P3), cmap=mycmap, vmin=0.0, vmax=dmax, extent=(left,right,bottom,top) );
plt.colorbar(shrink=0.5);
plt.plot( m1bar, d1bar, "wo", lw=3 );
plt.plot( gda_cvec(m1bar, m1bar), gda_cvec(d1max, d1max-0.1), "w-", lw=3 );
plt.plot( gda_cvec(m1min, m1min+0.1), gda_cvec(d1bar, d1bar), "w-", lw=3 );
plt.plot( mp, dp, "w:", lw=2 );
plt.xlabel("m");
plt.ylabel("d");
plt.show();
print("Caption: (Left) The prior pdf pA combining observation and prior information");
print("(Middle) the pdf of the theory Pg (right) the combined pdf pT");

gdapy06_14
No description has been provided for this image 
Caption: (Left) The prior pdf pA combining observation and prior information
(Middle) the pdf of the theory Pg (right) the combined pdf pT
In [16]:
# gdapy06_15

# Two different cases of:
#
# prior distribution pA(d1,m1) interacting with inexact
# theory pg(d1,m1) to produce total distribution pT(d1,m1)
#
# once case with very exact theory, other with very inexact one
 
print("gdapy06_15");

# d1-axis
Nd1 = 101;
d1min = 0.0;
d1max = 5.0;
Dd1 = (d1max-d1min)/(Nd1-1);
d1 = gda_cvec(d1min + Dd1*np.linspace(0,Nd1-1,Nd1));
 
# m1-axis
Nm1 = 101;
m1min = 0;
m1max = 5.0;
Dm1 = (m1max-m1min)/(Nm1-1);
m1 = gda_cvec( m1min + Dm1*np.linspace(0,Nm1-1,Nm1));

# Case 1: nearly exact theory
 
# setup for distribution P1=PA
P1=np.zeros((Nd1,Nm1));
d1bar = 2.25;
m1bar = 2.08;
bar = gda_cvec(d1bar, m1bar);
sd1 = 0.5;
sm1 = 1.0;
C1 = np.diag( gda_cvec(sd1**2, sm1**2).ravel() );
CI1 = la.inv(C1);
DC1 = la.det(C1);
# note not normaalized, so max is unity
for i in range(Nm1):
    for j in range(Nd1):
        x1 = gda_cvec(d1[i,0], m1[j,0]) - bar;
        r2 = np.matmul(x1.T,np.matmul(CI1,x1)); r2=r2[0,0];
        P1[i,j] = np.exp( -0.5*r2 );

# s-axis for parametric curve s
Ns = 51;
smin = 0.0;
smax = 5.0;
Ds = (smax-smin)/(Ns-1);
s = gda_cvec(smin + Ds*np.linspace(0,Ns-1,Ns));
 
# Pg follows parameric curve
dp = 1.0+s-2.0*np.power(s/smax,2);
mp = np.copy(s);
 
# P2=Pg distribution; not normalizable
P2 = np.zeros((Nd1,Nm1));
sg1 = 0.1;
sg2 = sg1**2;
for i in range(Nm1):
    for j in range(Nd1):
        r2 = np.power(d1[i,0]-dp,2) + np.power(m1[j,0]-mp,2);
        r2min=np.min(r2);
        P2[i,j] = exp(-r2min/sg2);

# P3=PT is product of PA and Pg
P3 = np.multiply(P1,P2);

# find max
Pi, Pj, P3max = gda_matrixmin(P3)
Pmaxm1 = m1min+Dm1*Pj;
Pmaxd1 = d1min+Dd1*Pi;

fig1 = plt.figure(1,figsize=(12,5));

# plot PA
ax1=plt.subplot(1,3,1);
mycmap=matplotlib.colormaps['jet'];
plt.axis( [d1min, d1max, m1min, m1max] );
left=d1min; right=d1max;
top=m1max; bottom=m1min;
dmax=np.max(P1);
ax1.invert_yaxis();
plt.imshow(np.flipud(P1), cmap=mycmap, vmin=0.0, vmax=dmax, extent=(left,right,bottom,top) );
plt.colorbar(shrink=0.5);
plt.plot( m1bar, d1bar, "wo", lw=3 );
plt.plot( gda_cvec(m1bar, m1bar), gda_cvec(d1max, d1max-0.1), "w-", lw=3 );
plt.plot( gda_cvec(m1min, m1min+0.1), gda_cvec(d1bar, d1bar), "w-", lw=3 );
plt.plot( mp, dp, "w:", lw=2 );
plt.xlabel("m");
plt.ylabel("d");

# plot Pg
ax1=plt.subplot(1,3,2);
mycmap=matplotlib.colormaps['jet'];
plt.axis( [d1min, d1max, m1min, m1max] );
left=d1min; right=d1max;
top=m1max; bottom=m1min;
dmax=np.max(P2);
ax1.invert_yaxis();
plt.imshow(np.flipud(P2), cmap=mycmap, vmin=0.0, vmax=dmax, extent=(left,right,bottom,top) );
plt.colorbar(shrink=0.5);
plt.plot( m1bar, d1bar, "wo", lw=3 );
plt.plot( gda_cvec(m1bar, m1bar), gda_cvec(d1max, d1max-0.1), "w-", lw=3 );
plt.plot( gda_cvec(m1min, m1min+0.1), gda_cvec(d1bar, d1bar), "w-", lw=3 );
plt.plot( mp, dp, "w:", lw=2 );
plt.xlabel("m");
plt.ylabel("d");

# plot PT
ax1=plt.subplot(1,3,3);
mycmap=matplotlib.colormaps['jet'];
plt.axis( [d1min, d1max, m1min, m1max] );
left=d1min; right=d1max;
top=m1max; bottom=m1min;
dmax=np.max(P3);
ax1.invert_yaxis();
plt.imshow(np.flipud(P3), cmap=mycmap, vmin=0.0, vmax=dmax, extent=(left,right,bottom,top) );
plt.colorbar(shrink=0.5);
plt.plot( m1bar, d1bar, "wo", lw=3 );
plt.plot( gda_cvec(m1bar, m1bar), gda_cvec(d1max, d1max-0.1), "w-", lw=3 );
plt.plot( gda_cvec(m1min, m1min+0.1), gda_cvec(d1bar, d1bar), "w-", lw=3 );
plt.plot( mp, dp, "w:", lw=2 );
plt.xlabel("m");
plt.ylabel("d");
plt.show();
print("Caption: Prior pdf pA(d1,m1) (left) interacting with");
print("inexact theory pg(d1,m1) (middle) to produce total pdf pT(d1,m1)");
print("(right). In this case, the theory is very exact.");

# Case 2: very inexact theory
 
# setup for distribution P1=PA
P1=np.zeros((Nd1,Nm1));
d1bar = 2.25;
m1bar = 2.08;
bar = gda_cvec(d1bar, m1bar);
sd1 = 0.5;
sm1 = 1.0;
C1 = np.diag( gda_cvec(sd1**2, sm1**2).ravel() );
CI1 = la.inv(C1);
DC1 = la.det(C1);
# note not normaalized, so max is unity
for i in range(Nm1):
    for j in range(Nd1):
        x1 = gda_cvec(d1[i,0], m1[j,0]) - bar;
        r2 = np.matmul(x1.T,np.matmul(CI1,x1)); r2=r2[0,0];
        P1[i,j] = np.exp( -0.5*r2 );

# s-axis for parametric curve s
Ns = 51;
smin = 0.0;
smax = 5.0;
Ds = (smax-smin)/(Ns-1);
s = gda_cvec(smin + Ds*np.linspace(0,Ns-1,Ns));
 
# Pg follows parameric curve
dp = 1.0+s-2.0*np.power(s/smax,2);
mp = np.copy(s);
 
# P2=Pg distribution; not normalizable
P2 = np.zeros((Nd1,Nm1));
sg1 = 2.0;
sg2 = sg1**2;
for i in range(Nm1):
    for j in range(Nd1):
        r2 = np.power(d1[i,0]-dp,2) + np.power(m1[j,0]-mp,2);
        r2min=np.min(r2);
        P2[i,j] = exp(-r2min/sg2);

# P3=PT is product of PA and Pg
P3 = np.multiply(P1,P2);

# find max
Pi, Pj, P3max = gda_matrixmin(P3)
Pmaxm1 = m1min+Dm1*Pj;
Pmaxd1 = d1min+Dd1*Pi;

fig1 = plt.figure(1,figsize=(12,5));

# plot PA
ax1=plt.subplot(1,3,1);
mycmap=matplotlib.colormaps['jet'];
plt.axis( [d1min, d1max, m1min, m1max] );
left=d1min; right=d1max;
top=m1max; bottom=m1min;
dmax=np.max(P1);
ax1.invert_yaxis();
plt.imshow(np.flipud(P1), cmap=mycmap, vmin=0.0, vmax=dmax, extent=(left,right,bottom,top) );
plt.colorbar(shrink=0.5);
plt.plot( m1bar, d1bar, "wo", lw=3 );
plt.plot( gda_cvec(m1bar, m1bar), gda_cvec(d1max, d1max-0.1), "w-", lw=3 );
plt.plot( gda_cvec(m1min, m1min+0.1), gda_cvec(d1bar, d1bar), "w-", lw=3 );
plt.plot( mp, dp, "w:", lw=2 );
plt.xlabel("m");
plt.ylabel("d");

# plot Pg
ax1=plt.subplot(1,3,2);
mycmap=matplotlib.colormaps['jet'];
plt.axis( [d1min, d1max, m1min, m1max] );
left=d1min; right=d1max;
top=m1max; bottom=m1min;
dmax=np.max(P2);
ax1.invert_yaxis();
plt.imshow(np.flipud(P2), cmap=mycmap, vmin=0.0, vmax=dmax, extent=(left,right,bottom,top) );
plt.colorbar(shrink=0.5);
plt.plot( m1bar, d1bar, "wo", lw=3 );
plt.plot( gda_cvec(m1bar, m1bar), gda_cvec(d1max, d1max-0.1), "w-", lw=3 );
plt.plot( gda_cvec(m1min, m1min+0.1), gda_cvec(d1bar, d1bar), "w-", lw=3 );
plt.plot( mp, dp, "w:", lw=2 );
plt.xlabel("m");
plt.ylabel("d");

# plot PT
ax1=plt.subplot(1,3,3);
mycmap=matplotlib.colormaps['jet'];
plt.axis( [d1min, d1max, m1min, m1max] );
left=d1min; right=d1max;
top=m1max; bottom=m1min;
dmax=np.max(P3);
ax1.invert_yaxis();
plt.imshow(np.flipud(P3), cmap=mycmap, vmin=0.0, vmax=dmax, extent=(left,right,bottom,top) );
plt.colorbar(shrink=0.5);
plt.plot( m1bar, d1bar, "wo", lw=3 );
plt.plot( gda_cvec(m1bar, m1bar), gda_cvec(d1max, d1max-0.1), "w-", lw=3 );
plt.plot( gda_cvec(m1min, m1min+0.1), gda_cvec(d1bar, d1bar), "w-", lw=3 );
plt.plot( mp, dp, "w:", lw=2 );
plt.xlabel("m");
plt.ylabel("d");
plt.show();
print("Caption: Prior pdf pA(d1,m1) (left) interacting with");
print("inexact theory pg(d1,m1) (middle) to produce total pdf pT(d1,m1)");
print("(right). In this case, the theory is very inexact");

gdapy06_15
No description has been provided for this image 
Caption: Prior pdf pA(d1,m1) (left) interacting with
inexact theory pg(d1,m1) (middle) to produce total pdf pT(d1,m1)
(right). In this case, the theory is very exact.
No description has been provided for this image 
Caption: Prior pdf pA(d1,m1) (left) interacting with
inexact theory pg(d1,m1) (middle) to produce total pdf pT(d1,m1)
(right). In this case, the theory is very inexact
In [17]:
# gdapy06_16

# total distribution pT(m1,d1) and the projected
# distribution p(m) = (integral) pT(m1,d1) d(d1)
# highlighting the possibility that the maximum
# liklihood points of the two distributions are
# at two different values of m1

print("gdapy06_16");

# d1-axis
Nd1 = 101;
d1min = 0.0;
d1max = 5.0;
Dd1 = (d1max-d1min)/(Nd1-1);
d1 = gda_cvec( d1min + Dd1*np.linspace(0,Nd1-1,Nd1));
 
# m1-axis
Nm1 = 101;
m1min = 0.0;
m1max = 5.0;
Dm1 = (m1max-m1min)/(Nm1-1);
m1 = gda_cvec(m1min + Dm1*np.linspace(0,Nm1-1,Nm1));
 
# distribution P1 = pA
P1=np.zeros((Nd1,Nm1));
d1bar = 2.25;
m1bar = 2.08;
bar = gda_cvec(d1bar, m1bar);
sd1 = 0.5;
sm1 = 1.0;
C1 = np.diag( gda_cvec(sd1**2, sm1**2).ravel() );
CI1 = la.inv(C1);
DC1 = la.det(C1);
# note not normaalized, so max is unity
for i in range(Nm1):
    for j in range(Nd1):
        x1 = gda_cvec( d1[i,0], m1[j,0] ) - bar;
        r2 = np.matmul(x1.T,np.matmul(CI1,x1)); r2=r2[0,0];
        P1[i,j] = exp( -0.5*r2 );

# s-axis for parametric curve
Ns = 51;
smin = 0.0;
smax = 5.0;
Ds = (smax-smin)/(Ns-1);
s = gda_cvec(smin + Ds*np.linspace(0,Ns-1,Ns));
 
# parameric curve
dp = 1.0+s-2.0*(s/smax)**2;
mp = np.copy(s);
 
# P2 = Pg
P2 = np.zeros((Nd1,Nm1));
sg1 = 0.35;
sg2 = sg1**2;
for i in range(Nm1):
    for j in range(Nd1):
        r2 = np.power(d1[i,0]-dp,2) + np.power(m1[j,0]-mp,2);
        r2min=np.min(r2);
        P2[i,j] = exp( -r2min/sg2 );

#P3 = pT with pT=pA pg
P3 = np.multiply(P1,P2);

# find maximum
Pi, Pj, P3max = gda_matrixmin(-P3);
P3max=-P3max;
Pmaxm1 = m1min+Dm1*Pj;
Pmaxd1 = d1min+Dd1*Pi;
 
# integrate over data
Pm = np.sum(P3,axis=0);
Pmmax = np.max(Pm);
Pmi = np.argmax(Pm);
Pmm = m1min+Dm1*Pmi;
norm=Dm1*np.sum(Pm);
Pm = Pm/norm;
 
# plot pT(d1,m1)
fig1 = plt.figure(1,figsize=(12,5));
ax1=plt.subplot(1,1,1);
mycmap=matplotlib.colormaps['jet'];
plt.axis( [d1min, d1max, m1min, m1max] );
left=d1min; right=d1max;
top=m1max; bottom=m1min;
dmax=np.max(P3);
ax1.invert_yaxis();
# for testing orientation
#P3=np.zeros((Nd1,Nm1)); id1=2; im1=2; P3[id1,im1]=dmax;
#P3=np.zeros((Nd1,Nm1)); id1=Nd1-2; im1=2; P3[id1,im1]=dmax;
#P3=np.zeros((Nd1,Nm1)); id1=2; im1=Nm1-2; P3[id1,im1]=dmax;
plt.imshow(np.flipud(P3), cmap=mycmap, vmin=0.0, vmax=dmax, extent=(left,right,bottom,top) );
plt.colorbar(shrink=0.5);
plt.plot( mp, dp, "w:", lw=2 );
plt.plot( m1bar, d1bar, "wo", lw=3 );
plt.plot( Pmaxm1, Pmaxd1, "ko", lw=3 );
plt.plot( gda_cvec(m1bar, m1bar), gda_cvec(d1max, d1max-0.1), "w-", lw=3 );
plt.plot( gda_cvec(m1min, m1min+0.1), gda_cvec(d1bar, d1bar), "w-", lw=3 );
plt.plot( gda_cvec(Pmaxm1, Pmaxm1), gda_cvec(d1max, d1max-0.1), "k-", lw=3 );
plt.plot( gda_cvec(m1min, m1min+0.1), gda_cvec(Pmaxd1, Pmaxd1), "k-", lw=3 );
plt.xlabel("m");
plt.ylabel("d");
plt.show();
print("Caption: Total pdf pT(d,m) with prior value (white circle), maximum");
print("likelihood point (black cirle), and theory (dashed line) superimposed.");
        
# plot p(m1)
fig2 = plt.figure(2,figsize=(12,5));
ax1=plt.subplot(1,1,1);
plt.axis( [m1min, m1max, 0, 1] );
plt.plot( m1, Pm, "b-", lw=3 );
plt.plot( Pmm, np.max(Pm), "ko", lw=3 );
plt.xlabel("m");
plt.ylabel("p(m)");
plt.show();
print("Caption: Projected pdf p(m)=(integral) pT(d,m) dd (blue curve) and");
print("maximum likelihood point (black cirle).");

gdapy06_16
No description has been provided for this image 
Caption: Total pdf pT(d,m) with prior value (white circle), maximum
likelihood point (black cirle), and theory (dashed line) superimposed.
No description has been provided for this image 
Caption: Projected pdf p(m)=(integral) pT(d,m) dd (blue curve) and
maximum likelihood point (black cirle).
In [18]:
# gdapy06_17
 
# sample chi-squared analysis
# distribution of errors P (Phi), E and L
# for a problem of a band-limited timeseries
# with redundant noisy observations of the
# time series and with a priori smoothness
# information
 
print("gdapy06_17");

# time t-axis
M=101;
Dt=1;
t = gda_cvec( Dt*np.linspace(0,M-1,M));
tmax = t[M-1,0];
 
# make a wiggly curve m(t) by starting out with random noise
# and bandpass filtering it around a narrow range of angular
# frquencies centered on w0.  The bandpass filtering is
# accomplished by the gda_chebyshevfilt() function, which
# though not mentioned in the text (sorry about that) is pretty
# easy to use and has many other applications.  Its arguments are
# output_timeseries = gda_chebyshevfilt( input_timeseries, Dt, f1, f2 )
# where Dt is the sampling of the timeseries (say in s) and
# (f1,f2) are the (low,high) size of the frequency band (say in Hz)
n = 10;
A = 2.0;
w0 = n*pi/tmax;
Dw = w0/10.0;
nse = np.random.normal(loc=0.0,scale=1.0,size=(M,1));
mtrue, uu, vv = gda_chebyshevfilt( nse, Dt, (w0-Dw)/(2*pi), (w0+Dw)/(2*pi) );
mtrue = A*mtrue/np.max(np.abs(mtrue));
 
# components of data kernel, all normalize to produce
# data of about the same value
 
# data kernel is just 3 independent observations of same point
N = 3*M;
G = np.concatenate((np.identity(M),np.identity(M),np.identity(M)),axis=0);
dtrue = np.matmul(G,mtrue);
sigma_d = 0.1;
dobs = dtrue + np.random.normal(loc=0.0,scale=sigma_d,size=(N,1));
sqrtcovdi = np.identity(N)/sigma_d;
 
# prior information: second derivative close to zero
tr = gda_cvec( 1.0, np.zeros((M-3,1)) );
tc = gda_cvec( 1.0, -2.0, 1.0, np.zeros((M-3,1)) );
H = (1.0/Dt**2)*la.toeplitz( tr.ravel(), tc.ravel() );
h = np.zeros((M-2,1));
K = M-2;
mddtrue = np.matmul(H,mtrue);
sigma_mdd = (w0**2)*A/sqrt(2.0);
sigma_mdd2 = np.std(mddtrue);  # for test purposes
sqrtcovhi = np.identity(K)/sigma_mdd;

# generalized least squares equation and its solution
F = np.concatenate( (np.matmul(sqrtcovdi,G), np.matmul(sqrtcovhi,H)), axis=0 );
f = np.concatenate( (np.matmul(sqrtcovdi,dobs), np.matmul(sqrtcovhi,h)), axis=0);
FTF = np.matmul(F.T,F);
FTf = np.matmul(F.T,f);
mest = la.solve(FTF,FTf);
FTFinv = la.inv(FTF);
Cm = FTFinv;
 
mddest = np.matmul(H,mest);
sigma_mdd3 = np.std(mddest);  # for test purposes
 
# plot m
fig1 = plt.figure(1,figsize=(12,5));
ax1=plt.subplot(1,1,1);
plt.axis( [0, tmax, -1.1*A, 1.1*A] );
plt.plot( t, dobs[0:M,0], "b.", lw=2 );
plt.plot( t, dobs[M:2*M,0], "b.", lw=2 );
plt.plot( t, dobs[2*M:3*M], "b.", lw=2 );
plt.plot( t, mtrue, "k-", lw=4 );
plt.plot( t, mest, "r-", lw=2 );
plt.xlabel("t (s)");
plt.ylabel("m");
plt.show();
print("Caption: Bandlimited time series, m(t). True (black)");
print("estimated (red), redundant observations (blue).");

# chi-squared values
e = np.matmul(sqrtcovdi,(np.matmul(G,mest)-dobs)); # normalized prediction error
l = np.matmul(sqrtcovhi,(np.matmul(H,mest)-h)); # normalized prior error
E = np.matmul(e.T,e); E=E[0,0]; # Eobs
L = np.matmul(l.T,l); L=L[0,0]; # Lobs
P = E+L;  # PHIobs
vP = N+K-M; # degrees of freeedom in PHI
vE = vP*N/(N+K); # approx degrees of freeedom in E
vL = vP*K/(N+K); # approx degrees of freeedom in L
 
# Chi-squared test of the Null Hypothesis
print("   ");
print("Analytic 95% confidence tests");
if ( P>(vP-2.0*sqrt(2.0*vP)) and P<(vP+2.0*sqrt(2.0*vP)) ):
    print("P %.2f vP %.2f null hypothesis cannot be rejected" % (P, vP) );
else:
    print("P %.2f vP %.2f null hypothesis can be rejected" % (P, vP));
    
if ( E>(vE-2.0*sqrt(2.0*vE)) and E<(vE+2.0*sqrt(2.0*vE)) ):
    print("E %.2f vE %.2f null hypothesis cannot be rejected" % (E, vE) );
else:
    print("E %.2f vE %.2f null hypothesis can be rejected" % (E, vE) );
    
if ( L>(vL-2.0*sqrt(2.0*vL)) and L<(vL+2.0*sqrt(2.0*vL)) ):
    print("L %.2f vL %.2f null hypothesis cannot be rejected" % (L, vL) );
else:
    print("L %.2f vL %.2f null hypothesis can be rejected" % (L, vL) );
print("   ");
 
# Now do it lots of times to create empirical pdf's
# p(PHI), p(E) and p(L)
Nreal = 10000;
Pvec = np.zeros((Nreal,1));
Evec = np.zeros((Nreal,1));
Lvec = np.zeros((Nreal,1));
for ireal in range(Nreal):
    dobs = dtrue + np.random.normal(loc=0.0,scale=sigma_d,size=(N,1));
    f = np.concatenate( (np.matmul(sqrtcovdi,dobs), np.matmul(sqrtcovhi,h)), axis=0);
    FTf = np.matmul(F.T,f);
    mest = la.solve(FTF,FTf);
    ei = np.matmul(sqrtcovdi,(np.matmul(G,mest)-dobs)); # normalized prediction error
    li = np.matmul(sqrtcovhi,(np.matmul(H,mest)-h)); # normalized prior error
    Ei = np.matmul(ei.T,ei); Ei=Ei[0,0]; # Eobs
    Li = np.matmul(li.T,li); Li=Li[0,0]# Lobs
    Evec[ireal,0] = Ei;
    Lvec[ireal,0] = Li;
    Pvec[ireal,0] = Ei+Li;
Evec = np.sort(Evec,axis=0);
Lvec = np.sort(Lvec,axis=0);
Pvec = np.sort(Pvec,axis=0);

# empirical pdf's via histograms
Nb=201;
binmin = 0;
binmax = 2*vP;
h, e = np.histogram(Pvec,Nb,(binmin,binmax));  # create histogram
Nbin = len(h);    # lengths of histogram
Ne = len(e);      # length of edges
edges =  gda_cvec(e[0:Ne-1]);         # convert e to column-vector
bins = edges + 0.5*(edges[1,0]-edges[0,0]); # centers of bins
Db = bins[1,0]-bins[0,0]; # size of bins
Phist = gda_cvec(h);            # empirical pdf
h, e = np.histogram(Evec,Nb,(binmin,binmax));  # create histogram
Ehist = gda_cvec(h);            # empirical pdf
h, e = np.histogram(Lvec,Nb,(binmin,binmax));  # create histogram
Lhist = gda_cvec(h);            # empirical pdf
 
# p(P)
fig2 = plt.figure(2,figsize=(12,5));
ax1=plt.subplot(1,1,1);
top = 1.1*max(Phist);
plt.axis( [binmin, binmax, 0, top ] );
plt.plot( bins,Phist, "k-", lw=3 );
plt.plot( gda_cvec(vP, vP), gda_cvec(0, top/5), "r-", lw=2 );
plt.plot( gda_cvec(vP-2*sqrt(2*vP), vP-2*sqrt(2*vP)), gda_cvec(0.0, top/10), "b-", lw=4 );
plt.plot( gda_cvec(vP+2*sqrt(2*vP), vP+2*sqrt(2*vP)), gda_cvec(0.0, top/10), "b-", lw= 4 );
plt.xlabel("P");
plt.ylabel("counts");
plt.show();
 
# p(E)
fig3 = plt.figure(3,figsize=(12,5));
ax1=plt.subplot(1,1,1);
top = 1.1*max(Ehist);
plt.axis( [binmin, binmax, 0, top ] );
plt.plot( bins,Ehist, "k-", lw=3 );
plt.plot( gda_cvec(vE, vE), gda_cvec(0, top/5), "r-", lw=2 );
plt.plot( gda_cvec(vE-2*sqrt(2*vE), vE-2*sqrt(2*vE)), gda_cvec(0.0, top/10), "b-", lw=4 );
plt.plot( gda_cvec(vE+2*sqrt(2*vE), vE+2*sqrt(2*vE)), gda_cvec(0.0, top/10), "b-", lw= 4 );
plt.xlabel("E");
plt.ylabel("counts");
plt.show();
 
# p(L)
fig4 = plt.figure(4,figsize=(12,5));
ax1=plt.subplot(1,1,1);
top = 1.1*max(Lhist);
plt.axis( [binmin, binmax, 0, top ] );
plt.plot( bins,Lhist, "k-", lw=3 );
plt.plot( gda_cvec(vL, vL), gda_cvec(0, top/5), "r-", lw=2 );
plt.plot( gda_cvec(vL-2*sqrt(2*vL), vL-2*sqrt(2*vL)), gda_cvec(0.0, top/10), "b-", lw=4 );
plt.plot( gda_cvec(vL+2*sqrt(2*vL), vL+2*sqrt(2*vL)), gda_cvec(0.0, top/10), "b-", lw= 4 );
plt.xlabel("L");
plt.ylabel("counts");
plt.show();

print("Caption: Empirical pdfs for (top) generalized error Phi,");
print("(middle) prediction error E and (bottom) error in prior information");
print("Also shown are mean z(red) and 95 percent confidence bounds (blue).");

# empirical confidence limits based on the empirical pdf's
print("   ");
print("Numerical 95% confidence tests");

k=np.where(P<Pvec);
i=k[0];
if( not np.any(i) ):
    p=0.0;
else: 
    p = float(i[0])/float(Nreal);
if ( p>0.025 and p<0.975):
    print("P %.2f vP %.2f null hypothesis cannot be rejected" % (P, vP) );
else:
    print("P %.2f vP %.2f null hypothesis can be rejected" % (P, vP) );
    
k=np.where(E<Evec);
i=k[0];
if( not np.any(i) ):
    p=0.0;
else: 
    p = float(i[0])/float(Nreal);
if ( p>0.025 and p<0.975):
    print("E %.2f vE %.2f null hypothesis cannot be rejected" % (E, vE) );
else:
    print("E %.2f vE %.2f null hypothesis can be rejected" % (E, vE) );
    
k=np.where(L<Lvec);
i=k[0];
if( not np.any(i) ):
    p=0.0;
else: 
    p = float(i[0])/float(Nreal);
if ( p>0.025 and p<0.975):
    print("L %.2f vL %.2f null hypothesis cannot be rejected" % (L, vL) );
else:
    print("L %.2f vL %.2f null hypothesis can be rejected" % (L, vL) );
    

gdapy06_17
No description has been provided for this image 
Caption: Bandlimited time series, m(t). True (black)
estimated (red), redundant observations (blue).
   
Analytic 95% confidence tests
P 314.01 vP 301.00 null hypothesis cannot be rejected
E 227.54 vE 226.87 null hypothesis cannot be rejected
L 86.47 vL 74.13 null hypothesis cannot be rejected
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image 
Caption: Empirical pdfs for (top) generalized error Phi,
(middle) prediction error E and (bottom) error in prior information
Also shown are mean z(red) and 95 percent confidence bounds (blue).
   
Numerical 95% confidence tests
P 314.01 vP 301.00 null hypothesis cannot be rejected
E 227.54 vE 226.87 null hypothesis cannot be rejected
L 86.47 vL 74.13 null hypothesis cannot be rejected
In [19]:
# gdapy06_18
 
# example F-test to assess difference in fit of two models
# the two models are linear and cubic fit of the same d(z)
# dataset
 
print("gdapy06_18");

# auxially variable z
N = 11;
z = gda_cvec(np.linspace(0.0,N-1,N)/(N-1));
dtrue = z - 0.5*np.power(z,2);
 
# simulated data
sigmad=0.03;
dobs = dtrue + np.random.normal(loc=0.0, scale=sigmad,size=(N, 1) );
 
fig1 = plt.figure(1,figsize=(12,5));
ax1=plt.subplot(1,1,1);
plt.axis( [0, 1, -1, 1] );
plt.plot(z,dobs,"ro",lw=3);
 
# plot the observed data
# fit 1, straight line
M=2;
G=np.concatenate((np.ones((N,1)), z), axis=1);
GTG = np.matmul(G.T,G);
GTd = np.matmul(G.T,dobs);
mestA=la.solve(GTG,GTd);
dpreA = np.matmul(G,mestA);
 
# plot the predicted data for the linear fit
plt.plot(z,dpreA,"b-",lw=3);
plt.title("linear fit");
plt.xlabel("z");
plt.ylabel("d(z)");
plt.show();
print("Caption: Linear fit (blue curve) to data (red circles).");
print("  ");
# linear error
EA = np.matmul((dobs-dpreA).T,dobs-dpreA); EA=EA[0,0];
vA = N-M; # degrees of freedom
print("linear error %f, degrees of freedom %d" % (EA, vA) );

# plot the observed data
fig2 = plt.figure(2,figsize=(12,5));
ax1=plt.subplot(1,1,1);
plt.axis( [0, 1, -1, 1] );
plt.plot(z,dobs,"ro",lw=3);
 
# fit 2, cubic
M=4;
G=np.concatenate((np.ones((N,1)),z,np.power(z,2),np.power(z,3)),axis=1);
GTG = np.matmul(G.T,G);
GTd = np.matmul(G.T,dobs);
mestB=la.solve(GTG,GTd);
dpreB = np.matmul(G,mestB);
 
# plot the predicted data for the cubic fit
plt.plot(z,dpreB,"b-",lw=3);
plt.title("cubic fit");
plt.xlabel("z");
plt.ylabel("d(z)");
plt.show();
print("Caption: Cubic fit (blue curve) to data (red circles)");
print("   ");

# error of cubic fit
EB = np.matmul((dobs-dpreB).T,dobs-dpreB); EB=EB[0,0];
vB = N-M; # degrees of freedom
print("cubic error %f, degrees of freedom %d" % (EB,vB) );

Fobs = (EA/vA) / (EB/vB); # F-value
print("1/F %f F %f" % (1.0/Fobs, Fobs));

# probability value associated with F-value
if( Fobs<1.0 ):
    Fobs=1.0/Fobs;

P = 1.0 - (st.f.cdf(Fobs,vA,vB)-st.f.cdf(1.0/Fobs,vA,vB));
Pleft = st.f.cdf(1.0/Fobs,vA,vB);
Pright = 1.0-st.f.cdf(Fobs,vA,vB);
 
print("P(F<%f) = %f" % (1.0/Fobs, Pleft) );
print("P(F>%f) = %f" % (Fobs, Pright) );
print("P(F<%f or F>%f) = %f" % (1.0/Fobs, Fobs, P));
if( (Pleft+Pright)<0.05 ):
    print("Null Hypothesis can be rejected to 95% confidence");
else:
    print("Null Hypothesis cannot be rejected to 95% confidence");

gdapy06_18
No description has been provided for this image 
Caption: Linear fit (blue curve) to data (red circles).
  
linear error 0.025540, degrees of freedom 9
No description has been provided for this image 
Caption: Cubic fit (blue curve) to data (red circles)
   
cubic error 0.004077, degrees of freedom 7
1/F 0.205219 F 4.872853
P(F<0.205219) = 0.015787
P(F>4.872853) = 0.024329
P(F<0.205219 or F>4.872853) = 0.040116
Null Hypothesis can be rejected to 95% confidence
In [20]:
# gdapy06_19
# example of F test for problem with prior information
 
print("gdapy06_19");
 
# dense set of model parameters
M1=101; # must be odd
 
# densely sampled time t variable
Dt1=1.0;
t1 = gda_cvec(Dt1*np.linspace(0,M1-1,M1));
t1max=t1[M1-1,0];
 
# sparse set of model parameters
M2=int(floor((M1+1)/2));
 
# sparsely sampled time t variable
Dt2=2.0*Dt1;
t2 = gda_cvec(Dt2*np.linspace(0,M2-1,M2));
 
# matrix D2S takes dense to sparse model parameters by decimation
# m2 = D2S*m1
D2S = np.zeros((M2, M1));
j=0;
for i in range(M2):
    D2S[i,j]=1.0;
    j=j+2;

# matrix S2D takes sparse to dense model parameters by linear interpolation
# m1 = S2D*m2
S2D = np.zeros((M1, M2));
j=0;
for i in range(0,M1,2):
    S2D[i,j]=1;
    j=j+1;
j=0;
for i in range(1,M1-1,2):
    S2D[i,j]=0.5;
    S2D[i,j+1]=0.5;
    j=j+1;

# make a densely sampled wiggly curve m(t) by starting out with random noise
# and bandpass filtering it around a narrow range of angular
# frquencies centered on w0.  The bandpass filtering is
# accomplished by the gda_chebyshevfilt() function, which
# though not mentioned in the text (sorry about that) is pretty
# easy to use and has many other applications.  Its arguments are
# output_timeseries = gda_chebyshevfilt( input_timeseries, Dt, f1, f2 )
# where Dt is the sampling of the timeseries (say in s) and
# (f1,f2) are the (low,high) size of the frequency band (say in Hz)
n = 10;
A = 2.0;
w0 = n*pi/t1max;
Dw = w0/10.0;
nse = np.random.normal(loc=0.0,scale=1.0,size=(M1,1));
mtrue1,uu,vv = gda_chebyshevfilt( nse, Dt1, (w0-Dw)/(2.0*pi), (w0+Dw)/(2.0*pi) );
mtrue1 = A*mtrue1/np.max(np.abs(mtrue1));
 
# data kernel is just 3 independent observations of same point
N = 3*M1;
G1 = np.concatenate((np.identity(M1),np.identity(M1),np.identity(M1)),axis=0);
G2 = np.concatenate((S2D, S2D, S2D), axis=0 );
dtrue = np.matmul(G1,mtrue1);
 
# observed data is noisy
sigma_d = 0.1;
dobs = dtrue + np.random.normal(loc=0.0,scale=sigma_d,size=(N,1));
sqrtcovdi = np.identity(N)/sigma_d;
 
# prior information: second derivative close to zero
K1 = M1-2;
tr = gda_cvec(1.0, np.zeros((M1-3,1)) );
tc = gda_cvec(1.0, -2.0, 1.0, np.zeros((M1-3)) );
H1 = (1.0/(Dt1**2))*la.toeplitz( tr.ravel(), tc.ravel() );
h1 = np.zeros((M1-2,1));
K2 = M2-2;
tr = gda_cvec( 1.0, np.zeros((M2-3,1)) );
tc = gda_cvec( 1.0, -2.0, 1.0, np.zeros((M2-3,1)) );
H2 = (1/(Dt2**2))*la.toeplitz( tr, tc );
h2 = np.zeros((M2-2,1));
 
# prior covariance of model parameters
sigma_mdd = (w0**2)*A/sqrt(2.0);
sqrtcovhi1 = np.identity(K1)/sigma_mdd;
sqrtcovhi2 = np.identity(K2)/sigma_mdd;
 
# overall least squares equation and its solution dense problem
F1=np.concatenate((np.matmul(sqrtcovdi,G1),np.matmul(sqrtcovhi1,H1)),axis=0);
f1=np.concatenate((np.matmul(sqrtcovdi,dobs),np.matmul(sqrtcovhi1,h1)),axis=0);
FTF1 = np.matmul(F1.T,F1);
FTf1 = np.matmul(F1.T,f1);
mest1 = la.solve(FTF1,FTf1);
# sparse problem
F2=np.concatenate((np.matmul(sqrtcovdi,G2),np.matmul(sqrtcovhi2,H2)),axis=0);
f2=np.concatenate((np.matmul(sqrtcovdi,dobs),np.matmul(sqrtcovhi2,h2)),axis=0);
FTF2 = np.matmul(F2.T,F2);
FTf2 = np.matmul(F2.T,f2);
mest2 = la.solve(FTF2,FTf2);
 
# plot m
fig1 = plt.figure(1,figsize=(12,5));
ax1=plt.subplot(1,1,1);
plt.axis( [0, t1max, -1.1*A, 1.1*A] );
plt.plot( t1, dobs[0:M1,0], "b.", lw=2 );
plt.plot( t1, dobs[M1:2*M1], "b.", lw=2 );
plt.plot( t1, dobs[2*M1:3*M1], "b.", lw=2 );
plt.plot( t1, mtrue1, "k-", lw=8 );
plt.plot( t1, mest1, "r-", lw=5 );
plt.plot( t2, mest2, "g-", lw=3 );
plt.xlabel("t (s)");
plt.ylabel("m");
plt.show();
print("Caption: Bandlimited time series, true (black), estimated");
print("sparse solution (red), estimated dense solution (green), data (blue dots)");
print("  ");

# error associated with dense problem
e1 = np.matmul(sqrtcovdi ,np.matmul(G1,mest1)-dobs);
l1 = np.matmul(sqrtcovhi1,np.matmul(H1,mest1)-h1);
E1 = np.matmul(e1.T,e1); E1=E1[0,0];
L1 = np.matmul(l1.T,l1); L1=L1[0,0];
P1 = E1+L1;
vP1 = N+K1-M1;
vE1 = vP1*N/(N+K1);
vL1 = vP1*K1/(N+K1);
print("Fit 1: E1 %f vE1 %d L1 %f vL1 %d" % (E1,vE1,L1,vL1) );
print("N %d K1 %d M1 %d" % (N,K1,M1) );
 
# error associated with sparse problem
e2 = np.matmul(sqrtcovdi, np.matmul(G2,mest2)-dobs);
l2 = np.matmul(sqrtcovhi2,np.matmul(H2,mest2)-h2);
E2 = np.matmul(e2.T,e2); E2=E2[0,0];
L2 = np.matmul(l2.T,l2); L2=L2[0,0];
P2 = E2+L2;
vP2 = N+K2-M2;
vE2 = vP2*N/(N+K2);
vL2 = vP2*K2/(N+K2);
print("Fit 2: E2 %f vE2 %d L2 %f vL2 %d" % (E2,vE2,L2,vL2) );
print("N %d K2 %d M2 %d" % (N,K2,M2) );
    
# F ratio
vA=vP1;
vB=vP2;
PhiA = P1;
PhiB = P2;
Fobs = (PhiA/vA) / (PhiB/vB);
if( Fobs<1.0 ):
    Fobs=1.0/Fobs;
print("1/F %f F %f" % (1.0/Fobs, Fobs) );

Pval = 1.0 - abs(st.f.cdf(1.0/Fobs,vA,vB)-st.f.cdf(Fobs,vA,vB));
Pleft = st.f.cdf(1.0/Fobs,vA,vB);
Pright = 1.0-st.f.cdf(Fobs,vA,vB);
print("P(F<%f) = %f" % (1.0/Fobs, Pleft) );
print("P(F>%f) = %f" % (Fobs, Pright) );
print("P(F<%f or F>%f) = %f" % (1.0/Fobs, Fobs, Pval) );

if( (Pleft+Pright)<0.05 ):
    print("Null Hypothesis can be rejected to 95% confidence");
else:
    print("Null Hypothesis cannot be rejected to 95% confidence");
    
# now do a whole lot of the same problems
# with different realizations of observational noise
Nreal = 10000;
FFvec = np.zeros((Nreal,1));
for ireal in range(Nreal):
    # add new noise to data
    dobs = dtrue + np.random.normal(loc=0.0,scale=sigma_d,size=(N,1));
    
    # dense problem
    f1=np.concatenate((np.matmul(sqrtcovdi,dobs),np.matmul(sqrtcovhi1,h1)),axis=0);
    FTf1 = np.matmul(F1.T,f1);
    mest1 = la.solve(FTF1,FTf1);
    e1 = np.matmul(sqrtcovdi, np.matmul(G1,mest1)-dobs);
    l1 = np.matmul(sqrtcovhi1,np.matmul(H1,mest1)-h1);
    PP1 = np.matmul(e1.T,e1)+np.matmul(l1.T,l1); PP1=PP1[0,0];
   
    # add new noise to data
    dobs = dtrue + np.random.normal(loc=0.0,scale=sigma_d,size=(N,1));
    # sparse problem
    f2=np.concatenate((np.matmul(sqrtcovdi,dobs),np.matmul(sqrtcovhi2,h2)),axis=0);
    FTf2 = np.matmul(F2.T,f2);
    mest2 = la.solve(FTF2,FTf2);
    e2 = np.matmul(sqrtcovdi, np.matmul(G2,mest2)-dobs);
    l2 = np.matmul(sqrtcovhi2,np.matmul(H2,mest2)-h2);
    PP2 = np.matmul(e2.T,e2)+np.matmul(l2.T,l2); PP2=PP2[0,0];
    
    # tabulate F ratio
    FFvec[ireal,0] = (PP1/vP1) / (PP2/vP2);

FFvec = np.sort(FFvec);
FFmean = np.mean(FFvec);

# empirical pdf's via histograms
Nb=101;
binmin = 0.5;
binmax = 1.5;
h, e = np.histogram(FFvec,Nb,(binmin,binmax));  # create histogram
Nbin = len(h);    # lengths of histogram
Ne = len(e);      # length of edges
edges =  gda_cvec(e[0:Ne-1]);         # convert e to column-vector
bins = edges + 0.5*(edges[1,0]-edges[0,0]); # centers of bins
Db = bins[1,0]-bins[0,0]; # size of bins
FFhist = gda_cvec(h);            # empirical pdf
FFhist = FFhist/(Db*np.sum(FFhist));
 
# plot pdf
fig1 = plt.figure(1,figsize=(12,5));
ax1=plt.subplot(1,1,1);
top = 1.1*np.max(FFhist);
plt.axis( [binmin, binmax, 0.0, top ] );
plt.plot( bins,FFhist, "k-", lw=3 );
plt.plot( bins,st.f.pdf(bins,vP1,vP2), "r-", lw=3 );
plt.plot( gda_cvec(Fobs,Fobs), gda_cvec(0, top/5), "k-", lw=4 );
plt.plot( gda_cvec(FFmean,FFmean), gda_cvec(0, top/5.0), "g-", lw=4 );
plt.plot( gda_cvec(st.f.ppf(0.5,vP1,vP2),st.f.ppf(0.5,vP1,vP2)),gda_cvec(0.0, top/5.0), "r-", lw=4 );
plt.plot( gda_cvec(st.f.ppf(0.025,vP1,vP2),st.f.ppf(0.025,vP1,vP2)), gda_cvec(0.0, top/8.0), "b-", lw=4 );
plt.plot( gda_cvec(st.f.ppf(0.975,vP1,vP2),st.f.ppf(0.975,vP1,vP2)), gda_cvec(0.0, top/8.0), "b-", lw=4 );
plt.xlabel("F");
plt.ylabel("p(F)");
plt.show();
print("Caption: Analytic pdf p(F) (red curve), empirical pdf (black curve)");
print("observed F (black bar), empirical mean (green bar), true mean (red bar)");
print("95 percent confidence (blue bars).");

gdapy06_19
No description has been provided for this image 
Caption: Bandlimited time series, true (black), estimated
sparse solution (red), estimated dense solution (green), data (blue dots)
  
Fit 1: E1 219.345135 vE1 226 L1 86.133665 vL1 74
N 303 K1 99 M1 101
Fit 2: E2 275.286438 vE2 259 L2 35.663579 vL2 41
N 303 K2 49 M2 51
1/F 0.982405 F 1.017910
P(F<0.982405) = 0.438859
P(F>1.017910) = 0.438859
P(F<0.982405 or F>1.017910) = 0.877718
Null Hypothesis cannot be rejected to 95% confidence
No description has been provided for this image 
Caption: Analytic pdf p(F) (red curve), empirical pdf (black curve)
observed F (black bar), empirical mean (green bar), true mean (red bar)
95 percent confidence (blue bars).
In [ ]: