In [1]:
# gdapy12_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:
# 2022/11/25 - Created by W. Menke
# 2023/07/27 - Tweaks by W. Menke
# 2024/05/23 -  Patches to accomodate new verions of numpy, scipy, matplotlib
#               patched dot product to return scalar value,
#               and fixed similar issues due to new unequivaence of scalar and 1x1 array, and
#               improved the output of gdapy12_07
# 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] = v[0:N,0];
        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;
In [4]:
# gdapy12_01
#
# Monte Carlo search
# applied to exemplary curve-fitting problem

# tunable parameters
# variance of data

print("gdapy12_01");
sd = 0.3;

# data are in a simple auxillary variable, x
N=40;
xmin=0;
xmax=1.0;
Dx=(xmax-xmin)/(N-1);
x = gda_cvec( np.linspace(0,Dx*(N-1),N) );

# two model parameters
M=2;

# true model parameters
mt = gda_cvec( 1.21, 1.54 );


# d(x)=g(x, m1, m2) with g = sin(w0*m(1)*x) + m(1)*m(2);
w0=20.0;
dtrue = np.sin(w0*mt[0,0]*x) + mt[0,0]*mt[1,0]; # true data

# observed data
sd2 = sd**2;
dobs = dtrue + np.random.normal(loc=0.0,scale=sd,size=(N,1)); 

# plot data
fig1 = plt.figure(1,figsize=(12,5));
ax1=plt.subplot(1,1,1);
plt.axis( [0, xmax, 0, 4 ] );
plt.plot(x,dtrue,"k-",lw=3);
plt.plot(x,dobs,"ro",lw=3);
plt.xlabel("x");
plt.ylabel("d");
plt.show();

# Define 2D grid (for plotting purposes only)
L = 101;
Dm = 0.02;
m1min=0.0;
m2min=0.0;
m1a = gda_cvec( m1min+Dm*np.linspace(0.0,L-1,L) );
m2a = gda_cvec( m2min+Dm*np.linspace(0.0,L-1,L) );
m1max = m1a[L-1,0];
m2max = m2a[L-1,0];

# tabulate total error E on grid (for plotting purposed only)
# Note m1 varies along rows of E and m2 along columns
E = np.zeros((L,L));
for j in range(L):
    for k in range(L):
        dpre = np.sin(w0*m1a[j,0]*x) + m1a[j,0]*m2a[k,0];
        e = dobs-dpre;
        myE=np.matmul(e.T,e); myE=myE[0,0];
        E[j,k] = myE/sd2;
Emax = np.max(E);

# Monte Carlo search
# save best solutions, so far
Niter = 10000;
mbest = np.zeros((M,Niter));
Ebest = np.zeros((Niter,1));
iterbest = np.zeros((Niter,1));
Nsaved = 0;
# starting solution
m0 = gda_cvec(0.1,0.1);  # start in upepr-left hand corner of plot
# predicted data
dpre0 = np.sin(w0*m0[0,0]*x) + m0[0,0]*m0[1,0];
# prediction error
e0 = dobs-dpre0;
E0 = np.matmul(e0.T,e0)/sd2; E0=E0[0,0];
# save this solution first solution, the best so far
mbest[0:M,Nsaved:Nsaved+1] = m0;
Ebest[Nsaved,0] = E0;
iterbest[Nsaved,0] = 0;
Nsaved=Nsaved+1;
for myiter in range(1,Niter):
    # randomly select trial solution
    m1trial = np.random.uniform(low=m1min,high=m1max);
    m2trial = np.random.uniform(low=m2min,high=m2max);
    mtrial = gda_cvec( m1trial, m2trial );
    # predict data
    dtrial = np.sin(w0*m1trial*x) + m1trial*m2trial;
    # prediction error
    etrial = dobs-dtrial;
    Etrial = np.matmul(etrial.T,etrial)/sd2; Etrial=Etrial[0,0];
    # accept trial solution if error decreases
    if( Etrial<E0 ):
        mbest[0:M,Nsaved:Nsaved+1] = mtrial;
        Ebest[Nsaved,0] = Etrial;
        iterbest[Nsaved,0] = myiter;
        Nsaved=Nsaved+1;
        # reset best-solution so far to this one
        E0=Etrial;
        m0 = mtrial;

# throw away unused part of arrays
mbest = np.copy( mbest[0:2,0:Nsaved] );
Ebest = np.copy( Ebest[0:Nsaved,0:1] );
iterbest = np.copy( iterbest[0:Nsaved,0:1] );

# plot error surface
fig1 = plt.figure(3,figsize=(12,5));
ax1=plt.subplot(1,1,1);
mycmap=matplotlib.colormaps['jet'];
plt.axis( [m2min, m2max, m1min, m1max] );
left=m2min;
right=m2max;
bottom=m1max;
top=m1min;
plt.imshow( E, cmap=mycmap, vmin=0, vmax=Emax, extent=(left,right,bottom,top) );
plt.plot( mbest[1,0:Nsaved], mbest[0,0:Nsaved], "w*", lw=2);
plt.plot( mbest[1,0:Nsaved], mbest[0,0:Nsaved], "w-", lw=2);
plt.colorbar();
plt.xlabel('m2');
plt.ylabel('m1');
plt.gca().invert_yaxis();
plt.show();

# plot progress
fig1 = plt.figure(1,figsize=(12,5));
ax1=plt.subplot(3,1,1);
K = np.max(iterbest);
plt.axis( [0, K, 0, 4 ] );
plt.plot( iterbest[:,0], mbest[0,:],"k-",lw=1);
plt.plot( gda_cvec(0,K),gda_cvec(mt[0,0],mt[0,0]),"r:",lw=1);
plt.plot( iterbest[:,0], mbest[0,:],"k*",lw=3);
plt.xlabel("iteration, i");
plt.ylabel("m1");
ax1=plt.subplot(3,1,2);
plt.axis( [0, K, 0, 4 ] );
plt.plot( iterbest[:,0], mbest[1,:],"k-",lw=1);
plt.plot( gda_cvec(0,K),gda_cvec(mt[1,0],mt[1,0]),"r:",lw=1);
plt.plot( iterbest[:,0], mbest[1,:],"k*",lw=3);
plt.xlabel("iteration, i");
plt.ylabel("m2");
ax1=plt.subplot(3,1,3);
plt.axis( [0, K, np.min(np.log10(Ebest)), np.max(np.log10(Ebest)) ] );
plt.plot( iterbest[:,0], np.log10(Ebest[:,0]),"k-",lw=1);
plt.plot( iterbest[:,0], np.log10(Ebest[:,0]),"k*",lw=3);
plt.xlabel("iteration, i");
plt.ylabel("log10(E)");
plt.show();

print("%d iterations, %d accepted" % (Niter,Nsaved) );
gdapy12_01
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
10000 iterations, 10 accepted
In [7]:
# gdapy12_02
#
# Simulated Annealing search
# applied to exemplary curve-fitting problem

print("gdapy12_02");

# tunable parameters
# variance of data
sd = 0.3;
# neighborhood parameter
sq = 0.3;

# data are in a simple auxillary variable, x
N=40;
xmin=0;
xmax=1.0;
Dx=(xmax-xmin)/(N-1);
x = gda_cvec( np.linspace(0,Dx*(N-1),N) );

# two model parameters
M=2;

# true model parameters
mt = gda_cvec( 1.21, 1.54 );

# d(x)=g(x, m1, m2) with g = sin(w0*m(1)*x) + m(1)*m(2);
w0=20.0;
dtrue = np.sin(w0*mt[0,0]*x) + mt[0,0]*mt[1,0]; # true data

# observed data
sd2 = sd**2;
dobs = dtrue + np.random.normal(loc=0.0,scale=sd,size=(N,1)); 

# plot data
fig1 = plt.figure(1,figsize=(12,5));
ax1=plt.subplot(1,1,1);
plt.axis( [0, xmax, 0, 4 ] );
plt.plot(x,dtrue,"k-",lw=3);
plt.plot(x,dobs,"ro",lw=3);
plt.xlabel("x");
plt.ylabel("d");
plt.show();

# Define 2D grid (for plotting purposes only)
L = 101;
Dm = 0.02;
m1min=0.0;
m2min=0.0;
m1a = gda_cvec( m1min+Dm*np.linspace(0.0,L-1,L) );
m2a = gda_cvec( m2min+Dm*np.linspace(0.0,L-1,L) );
m1max = m1a[L-1,0];
m2max = m2a[L-1,0];

# tabulate total error psi on grid.
# Note m1 varies along rows of E and m2 along columns
E = np.zeros((L,L));
for j in range(L):
    for k in range(L):
        dpre = np.sin(w0*m1a[j,0]*x) + m1a[j,0]*m2a[k,0];
        e = dobs-dpre;
        myE = np.matmul(e.T,e); myE=myE[0,0];
        E[j,k] = myE/sd2;
Emax = np.max(E);

# Simulted annealing
# save all solutions and their error
Niter = 2000;
msave = np.zeros((M,Niter));
Esave = np.zeros((Niter,1));
n = np.arange(Niter);

# temperature 
Tstart = 100.0;
Tend = 0.001;

T = gda_cvec(np.exp(np.linspace(log(Tstart),log(Tend),Niter)));

# starting solution
m0 = gda_cvec(0.2,0.2);  # start in upepr-left hand corner of plot
# predicted data
d0 = np.sin(w0*m0[0,0]*x) + m0[0,0]*m0[1,0];
# prediction error
e0 = dobs-d0;
E0 = np.matmul(e0.T,e0); E0=E0[0,0];
# likelihood
logpm0 = -E0/(T[0,0]*sd2);
# save this solution first solution, the best so far
msave[0:M,0:1] = m0;
Esave[0,0] = E0;
Naccept = 0;
for myiter in range(1,Niter):
    # randomly select trial solution
    mp = np.random.normal(loc=m0,scale=sq,size=(M,1));
    # predict data
    dp = np.sin(w0*mp[0,0]*x) + mp[0,0]*mp[1,0];
    # prediction error
    ep = dobs-dp;
    Ep = np.matmul(ep.T,ep); Ep=Ep[0,0]
    logpmp = -0.5*Ep/(T[myiter,0]*sd2);
    logalpha = logpmp - logpm0;
    # standard metropolis acceptance test
    if( logalpha >= 0.0 ):
        m0 = mp;
        E0 = Ep;
        logpm0 = logpmp;
        Naccept=Naccept+1;
    else:
        alpha = exp(logalpha);
        r = np.random.uniform(low=0.0,high=1.0);
        if( alpha > r ):
            m0 = mp;
            E0 = Ep;
            logpm0 = logpmp;
            Naccept=Naccept+1;
        else:
            mp = m0;
            Ep = E0;
            logpmp = logpm0;
    msave[0:M,myiter:myiter+1] = mp;
    Esave[myiter,0] = Ep;

# plot error surface
fig1 = plt.figure(3,figsize=(12,5));
ax1=plt.subplot(1,1,1);
mycmap=matplotlib.colormaps['jet'];
plt.axis( [m2min, m2max, m1min, m1max] );
left=m2min;
right=m2max;
bottom=m1max;
top=m1min;
plt.imshow( E, cmap=mycmap, vmin=0, vmax=Emax, extent=(left,right,bottom,top) );
plt.plot( msave[1,0:Niter], msave[0,0:Niter], "w*", lw=2);
plt.plot( msave[1,0:Niter], msave[0,0:Niter], "w-", lw=2);
plt.plot( mt[1,0], mt[0,0], "g*", lw=2);
plt.plot( msave[1,Niter-1], msave[0,Niter-1], "r*", lw=2);
plt.colorbar();
plt.xlabel('m2');
plt.ylabel('m1');
plt.gca().invert_yaxis();
plt.show();

# plot progress
fig1 = plt.figure(1,figsize=(12,5));
ax1=plt.subplot(4,1,1);
plt.axis( [0, Niter, 0, 4 ] );
plt.plot( n, msave[0,:],"k-",lw=1);
plt.plot( gda_cvec(0,Niter),gda_cvec(mt[0,0],mt[0,0]),"r:",lw=1);
plt.xlabel("iteration, i");
plt.ylabel("m1");
ax1=plt.subplot(4,1,2);
plt.axis( [0, Niter, 0, 4 ] );
plt.plot( n, msave[1,:],"k-",lw=1);
plt.plot( gda_cvec(0,Niter),gda_cvec(mt[1,0],mt[1,0]),"r:",lw=1);
plt.xlabel("iteration, i");
plt.ylabel("m2");
ax1=plt.subplot(4,1,3);
plt.axis( [0, Niter, np.min(np.log10(Esave)), np.max(np.log10(Esave)) ] );
plt.plot( n, np.log10(Esave[:,0]),"k-",lw=1);
plt.xlabel("iteration, i");
plt.ylabel("log10(E)");
ax1=plt.subplot(4,1,4);
plt.axis( [0, Niter, np.min(np.log(T)), np.max(np.log(T)) ] );
plt.plot( n, np.log10(T[:,0]),"k-",lw=1);
plt.xlabel("iteration, i");
plt.ylabel("log10(T)");
plt.show();

print("%d iterations, %d accepts" % (Niter,Naccept) );
gdapy12_02
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
2000 iterations, 49 accepts
In [8]:
# gdapy12_03
#
# Using the Metropolis-Hastibgs algorith to sample the likelihood
# associated with the exemplary curve-fitting problem

print("gdapy12_03");

# tunable parameters
# variance of data
sd = 0.3;
# variance of prior information
sm = 0.5;
# neighborhood parameter
sn = 0.05;

# data are in a simple auxillary variable, x
N=40;
xmin=0;
xmax=1.0;
Dx=(xmax-xmin)/(N-1);
x = gda_cvec( np.linspace(0,Dx*(N-1),N) );

# two model parameters
M=2;

# true model parameters
mt = gda_cvec( 1.21, 1.54 );

# prior model parameters and their variance
mA = gda_cvec( 1.43, 1.50 );
sm2 = sm**2;

# d(x)=g(x, m1, m2) with g = sin(w0*m(1)*x) + m(1)*m(2);
w0=20.0;
dtrue = np.sin(w0*mt[0,0]*x) + mt[0,0]*mt[1,0]; # true data

# observed data
sd2 = sd**2;
dobs = dtrue + np.random.normal(loc=0.0,scale=sd,size=(N,1)); 

# plot data
fig1 = plt.figure(1,figsize=(12,5));
ax1=plt.subplot(1,1,1);
plt.axis( [0, xmax, 0, 4 ] );
plt.plot(x,dtrue,"k-",lw=3);
plt.plot(x,dobs,"ro",lw=3);
plt.xlabel("x");
plt.ylabel("d");
plt.show();

# Define 2D grid
L = 101;
Dm = 0.02;
m1min=0.0;
m2min=0.0;
m1a = gda_cvec( m1min+Dm*np.linspace(0.0,L-1,L) );
m2a = gda_cvec( m2min+Dm*np.linspace(0.0,L-1,L) );
m1max = m1a[L-1,0];
m2max = m2a[L-1,0];

# tabulate total error psi on grid.
# Note m1 varies along rows of E and m2 along columns
Psi = np.zeros((L,L));
for j in range(L):
    for k in range(L):
        dpre = np.sin(w0*m1a[j,0]*x) + m1a[j,0]*m2a[k,0];
        e = dobs-dpre;
        l = mA-gda_cvec(m1a[j,0], m2a[k,0]);
        myE = np.matmul(e.T,e); myE=myE[0,0];
        myL = np.matmul(l.T,l); myL=myL[0,0];
        Psi[j,k] = myE/sd2+myL/sm2;
Psimax = np.max(Psi);

# Metropolis Hastings
Nburnt=100000;  # burn in
Niter=1000000+Nburnt;
m = np.zeros((M,Niter));
m1c = (m1min + m1max)/2.0;
m2c = (m2min + m2max)/2.0;
mold = gda_cvec( m1c, m2c );
m[0:M,0:1] = mold;
dold = np.sin(w0*mold[0,0]*x) + mold[0,0]*mold[1,0];
e = dobs-dold;
Eold = np.matmul(e.T,e);
l = mA-mold;
Lold = np.matmul(l.T,l);
logpold = -0.5*Eold/sd2-0.5*Lold/sm2;
for myiter in range(1,Niter):

    # perturbation
    Dm = np.random.normal(loc=0.0,scale=sn,size=(M,1));
    mnew = mold + Dm;
    
    # prior error
    l = mA-mnew;
    Lnew = np.matmul(l.T,l); Lnew=Lnew[0,0];
    
    # prediction
    dnew = np.sin(w0*mnew[0,0]*x) + mnew[0,0]*mnew[1,0];
    e = dobs-dnew;
    Enew = np.matmul(e.T,e); Enew=Enew[0,0];
    logpnew = -0.5*Enew/sd2-0.5*Lnew/sm2;
    
    # acceptance test
    
    if( (logpnew-logpold) >= 0.0):
        m[0:M,myiter:myiter+1] = mnew;
        mold = mnew;
        Eold = Enew;
        Lold = Lnew;
        logpold = logpnew;
    else:
        alpha = exp( logpnew - logpold );
        r = np.random.uniform(low=0.0,high=1.0);
        if( alpha > r ):
            m[0:M,myiter:myiter+1] = mnew;
            mold = mnew;
            Eold = Enew;
            Lold = Lnew;
            logpold = logpnew;
        else:
            m[0:M,myiter:myiter+1] = mold;

# keep only values past burn in
m = np.copy( m[0:M,Nburnt:Niter]);
i, Niter = np.shape(m);
# estimate is the mean
mest = gda_cvec( np.mean(m,axis=1) );

# sample covariance
Cm = np.cov( m );
# 95% confidence of the mean
Dm95 = gda_cvec( 2.0*np.sqrt(np.diag(Cm)/Niter) );
            
# plot error surface
fig2 = plt.figure(2,figsize=(12,5));
ax1=plt.subplot(1,1,1);
mycmap=matplotlib.colormaps['jet'];
plt.axis( [m2min, m2max, m1min, m1max] );
left=m2min;
right=m2max;
bottom=m1max;
top=m1min;
plt.imshow( Psi, cmap=mycmap, vmin=0, vmax=Psimax, extent=(left,right,bottom,top) );
plt.plot( mest[1,0], mest[0,0], "go", lw=3);
plt.plot( gda_cvec(mest[1,0]-Dm95[1,0],mest[1,0]+Dm95[1,0]),
          gda_cvec(mest[0,0],mest[0,0]), "g-", lw=1);
plt.plot( gda_cvec(mest[1,0],mest[1,0]),
          gda_cvec(mest[0,0]-Dm95[0,0],mest[0,0]+Dm95[0,0]), "g-", lw=1);
plt.plot( mA[1,0], mA[0,0], "y*", lw=3);
plt.colorbar();
plt.xlabel('m2');
plt.ylabel('m1');
plt.gca().invert_yaxis();
plt.show();

# plot error surface
fig3 = plt.figure(3,figsize=(12,5));
ax1=plt.subplot(1,1,1);
mycmap=matplotlib.colormaps['jet'];
plt.axis( [m2min, m2max, m1min, m1max] );
left=m2min;
right=m2max;
bottom=m1max;
top=m1min;
plt.imshow( Psi, cmap=mycmap, vmin=0, vmax=Psimax, extent=(left,right,bottom,top) );
plt.plot( m[1,0:Niter], m[0,0:Niter], "w.", lw=1)
plt.plot( mA[1,0], mA[0,0], "y*", lw=3);
plt.colorbar();
plt.xlabel('m2');
plt.ylabel('m1');
plt.gca().invert_yaxis();
plt.show();

# histogram of predicted data
NHx=512;
NHv=128;
xH=gda_cvec(np.linspace(0,xmax,NHx))
vmin=0.0;
vmax=4.0;
Dv = (vmax-vmin)/(NHv-1);
H=np.zeros((NHv,NHx));
n = np.arange(NHx);
for myiter in range(Niter):
    dH = np.sin(w0*m[0,myiter]*xH) + m[0,myiter]*m[1,myiter];
    # painful!
    j = ((dH-vmin)/Dv).astype(int);  # convert data to integer array index
    r = np.where((j>=0)&(j<NHv));    # exclude indices that are out of bound
    k = r[0];
    jg = j[k].ravel();  # good data value indices
    ig = n[k].ravel();  # corresponding x position indices
    H[jg,ig]=H[jg,ig]+1; # accumulate histogram

# plot data
fig1 = plt.figure(1,figsize=(12,5));
ax1=plt.subplot(1,1,1);
plt.axis( [0, xmax, vmin, vmax ] );
left=0.0;
right=xmax;
bottom=vmin;
top=vmax;
Hmin=0.0;
Hmax = np.max(H);
plt.imshow( H, cmap=mycmap, vmin=Hmin, vmax=Hmax, origin='lower', extent=(left,right,bottom,top), aspect='auto' );
plt.plot(x,dobs,"wo",lw=3);
plt.xlabel("x");
plt.ylabel("d");
plt.colorbar();
plt.show();
gdapy12_03
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
In [9]:
# gdapy12_04
#
# Using the Metropolis-Hastibgs algorith to sample the likelihood
# associated with the exemplary Laplace transform problem

print("gdapy12_04");

# data are in a simple auxillary variable, z

# z-axis
M=11;
zmin=0;
zmax=1.0;
Dz=(zmax-zmin)/(M-1);
z = gda_cvec( np.linspace(0,Dz*(M-1),M) );

# true model
mtrue = np.ones((M,1))+np.sin(2*pi*z/zmax);

# data kernel, exponentially decaying
N=M;
G = np.zeros((N,M));
c = 6.0*z; # decay rates
for i in range(M):
    G[i:i+1,0:M] = np.exp(-c[i,0]*z ).T;

# prior model parameters and their variance
mA = np.ones((M,1));
sm = 0.1;
sm2 = sm**2;

# d = Gm
dtrue = np.matmul(G,mtrue);

# observed data
sd=0.01;
sd2 = sd**2;
dobs = dtrue + np.random.normal(loc=0.0,scale=sd,size=(N,1)); 

gda_draw( "title G", G, "title m", mtrue, "=", "title dobs", dobs );

# plot data
fig1 = plt.figure(1,figsize=(12,5));
ax1=plt.subplot(1,1,1);
plt.axis( [0, zmax, 0, np.max(dobs) ] );
plt.plot(z,dtrue,"k-",lw=3);
plt.plot(z,dobs,"ro",lw=3);
plt.xlabel("z");
plt.ylabel("d");
plt.show();
           
# neighborhood parameter
sn = 0.2;
sn = 0.05;

# Metropolis Hastings
Niter=1000000;
Nburnt=1000;  # burn in
m = np.zeros((M,Niter));
mold = np.random.uniform(low=0.8,high=1.0,size=(M,1));
m[0:M,0:1] = mold;
dold = np.matmul(G,mold);
e = dobs-dold;
Eold = np.matmul(e.T,e); Eold=Eold[0,0];
l = mA-mold;
Lold = np.matmul(l.T,l); Lold=Lold[0,0];
logpold = -0.5*Eold/sd2-0.5*Lold/sm2;
for myiter in range(1,Niter):

    # perturbation
    Dm = np.random.normal(loc=0.0,scale=sn,size=(M,1));
    mnew = mold + Dm;
    
    # prior error
    l = mA-mnew;
    Lnew = np.matmul(l.T,l); Lnew=Lnew[0,0];
    
    # prediction
    dnew = np.matmul(G,mnew);
    e = dobs-dnew;
    Enew = np.matmul(e.T,e); Enew=Enew[0,0];
    logpnew = -0.5*Enew/sd2-0.5*Lnew/sm2;
    
    # acceptance test
    Dlogp = logpnew - logpold;
    if( Dlogp>0.0 ): # prevent overflow
        alpha=1.1;
    else:
        alpha = exp( Dlogp );
    if( alpha > 1.0):
        m[0:M,myiter:myiter+1] = mnew;
        mold = mnew;
        Eold = Enew;
        Lold = Lnew;
        logpold = logpnew;
    else:
        r = np.random.uniform(low=0.0,high=1.0);
        if( alpha > r ):
            m[0:M,myiter:myiter+1] = mnew;
            mold = mnew;
            Eold = Enew;
            Lold = Lnew;
            logpold = logpnew;
        else:
            m[0:M,myiter:myiter+1] = mold;

# keep only values past burn in
m = np.copy( m[0:M,Nburnt:Niter]);
i, Niter = np.shape(m);
# estimate is the mean
mest = gda_cvec( np.mean(m,axis=1) );

# sample covariance
Cm = np.cov( m );
# 2 * std dev
Dm95 = 2.0*gda_cvec(np.sqrt(np.diag(Cm)));

# plot model
fig1 = plt.figure(1,figsize=(12,5));
ax1=plt.subplot(1,1,1);
plt.axis( [0, zmax, 0, 3 ] );
plt.plot(z,mtrue,"k-",lw=3);
plt.plot(z,mest,"r-",lw=2);
plt.plot( z, mest-Dm95, "g-", lw=1);
plt.plot( z, mest+Dm95, "g-", lw=1);
plt.xlabel("z");
plt.ylabel("m");
plt.show();

# histogram
Nbins = 101;
binmin=0.0;
binmax=3.0;
Dbin = (binmax-binmin)/(Nbins-1);
bins = gda_cvec(np.linspace(binmin,binmax,Nbins));
H = np.zeros((Nbins,M));
for iter in range(Niter):
    for i in range(M):
        m0 = m[i,iter];
        k = int((m0-binmin)/Dbin);
        if( (k<0) or (k>=Nbins) ):
            continue;
        H[k,i]=H[k,i]+1.0;

# plot error surface
fig3 = plt.figure(3,figsize=(12,5));
ax1=plt.subplot(1,1,1);
mycmap=matplotlib.colormaps['jet'];
plt.axis( [0, zmax, binmin, binmax] );
left=0;
right=zmax;
bottom=binmax;
top=binmin;
Hroot = np.sqrt(H);
Hmax = np.max(Hroot);
plt.imshow( Hroot, cmap=mycmap, vmin=0, vmax=Hmax,
           extent=(left,right,bottom,top), aspect='auto' );
plt.plot( z, mtrue, "y-", lw=3);
plt.plot( z, mest, "w-", lw=3);
plt.colorbar();
plt.xlabel('z');
plt.ylabel('m');
plt.show();

count=np.sum(np.logical_and(np.greater(m[2,:],m[1,:]),np.greater(m[2,:],m[3,:])));
print("%.6f percent have m2>m1 and m2>m3" % (100*count/Niter) );
gdapy12_04
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
65.856056 percent have m2>m1 and m2>m3
In [23]:
# gdapy12_05
#
# Using the extended Metropolis-Hastibgs algorith to sample the a
# trans-dimensional Normal pdf  with two data types:
#
# p(m) = [1d Normal with probability beta] + [2d Normal with probability (1-beta)]

print("gdapy12_05");

beta = 0.2; # probaility of dim 1

# dim-1 pdf is normal with this variance and normalization facto
sm1 = 0.5;
sm12 = sm1**2;
N1 = beta/sqrt(2.0*pi*sm12);
logN1 = log(N1);

# dim-2 pdf is normal with this variance and normalization factor
sm2 = 1.0;
sm22 = sm2**2;
N2 = (1.0-beta)/(2.0*pi*sm22);
logN2 = log(N2);

# Metropolis Hastings
Niter=1000000;
Nburnt=1000;  # burn in
doswitch=500;  # switch dimensions

# stores links in chain
m = np.zeros((3,Niter)); # m[2,.] 1 or 2 signifies dimension

dimold = 2;  # starting dimenion
dimnew = dimold;

# assign first link in chain
if( dimold==1 ):
    mold = np.random.uniform(low=-1.0,high=1.0,size=(1,1));
    logpold = logN1-0.5*(mold[0,0]**2)/sm12;
    m[0,0] = mold[0,0];
    m[0,1] = 0.0;
    m[0,2] = dimold;
elif( dimold==2 ):
    mold = gda_cvec(np.random.uniform(low=-1.0,high=1.0,size=(2,1)));
    logpold = logN2-0.5*(mold[0,0]**2+mold[1,0]**2)/sm22;
    m[0,0] = mold[0,0];
    m[0,1] = mold[1,0];
    m[0,2] = dimold;

# random verctors u are normal with these variance and normlizations
su1 = 0.2;
su12 = su1**2;
Nu1 = 1/(sqrt(2.0*pi)*su1);
logNu1 = log(Nu1);

su2 = 0.2;
su22 = su2**2;
Nu2 = 1/(sqrt(2.0*pi)*su2);
logNu2 = log(Nu2);

su3 = 1.0;
su32 = su3**2;
Nu3 = 1/(sqrt(2.0*pi)*su3);
logNu3 = log(Nu3);

# Transdimensional Metropolis Hastings 
for myiter in range(1,Niter):
    
    # switch dimensions every doswitch iterations
    if( (myiter%doswitch)==0 ):
        if( dimnew == 1 ):
            dimnew = 2;
        elif ( dimnew == 2 ):
            dimnew = 1;

    # three random vectors u1, u2, u3 with transformation
    # for the dim-1 to dim-2 transformation
    # m1' = [ 1 1 0 0 ] [m1]        (m1' = m1 + u1)
    # m2' = [ 0 0 0 1 ] [u1]        (m2' = u3)
    # u1' = [ 0 1 0 0 ] [u2]        (u1' = u1)
    # u2' = [ 0 0 1 0 ] [u3]        (u2' = u2)
    # which has a Jacobian of unity
    u1 = np.random.normal(loc=0.0,scale=su1); # g1
    u2 = np.random.normal(loc=0.0,scale=su2); # g2
    u3 = np.random.normal(loc=0.0,scale=su3); # g3
       
    if( (dimold==1) and (dimnew==1) ):
        mnew = np.zeros((1,1));
        mnew[0,0] = mold[0,0] + u1;
        logpnew = logN1-0.5*(mnew[0,0]**2)/sm12;
        logalpha=logpnew-logpold;
    elif( (dimold==1) and (dimnew==2) ):
        # pnew = p12'  g1'  g2'
        # pold = p1    g1   g2  g3
        # but since g1=g1' and g2=g2'
        # so pnew/pold = p12' / (p1 g3)
        mnew = np.zeros((2,1));
        logg3 = logNu3-0.5*(u3**2)/su32;
        mnew[0,0]= mold[0,0] + u1;
        mnew[1,0]= u3;
        logpnew = logN2-0.5*(mnew[0,0]**2+mnew[1,0]**2)/sm22;
        logalpha=logpnew-(logpold+logg3);
    elif( (dimold==2) and (dimnew==2) ):
        mnew = np.zeros((2,1));
        mnew[0,0]= mold[0,0] + u1;
        mnew[1,0]= mold[1,0] + u2;
        logpnew = logN2-0.5*(mnew[0,0]**2+mnew[1,0]**2)/sm22;
        logalpha=logpnew-logpold;
    elif( (dimold==2) and (dimnew==1) ):
        mnew = np.zeros((1,1));
        mnew[0,0] = mold[0,0] - u1;
        logpnew = logN1-0.5*(mnew[0,0]**2)/sm12;
        up3 = mold[1,0];
        loggp3 = logNu3-0.5*(up3**2)/su32;
        logalpha=(loggp3+logpnew)-logpold;

    if( logalpha>0.0 ):
        alpha=1.1;
    else:
        alpha=exp(logalpha);
    if( alpha > 1.0):
        dimold = dimnew;
        mold = mnew;
        logpold = logpnew;
    else:
        r = np.random.uniform(low=0.0,high=1.0);
        if( alpha > r ):
            dimold = dimnew;
            mold = mnew;
            logpold = logpnew;
        else:
            dimnew = dimold;
            mnew = mold;
            logpnew = logpold;
    if( dimnew == 1 ):
        m[0,myiter]=mnew[0,0];
        m[1,myiter]=0.0;
        m[2,myiter]=dimnew;
    elif( dimnew == 2 ):
        m[0,myiter]=mnew[0,0];
        m[1,myiter]=mnew[1,0];
        m[2,myiter]=dimnew;

# keep only values past burn in
m = np.copy( m[0:3,Nburnt:Niter]);
i, Niter = np.shape(m);
# estimate is the mean

Nbins = 41;
binmin = -3.0;
binmax = 3.0;
Dbin = (binmax-binmin)/(Nbins-1);
hbin = gda_cvec(np.linspace(binmin,binmax,Nbins));

H1 = np.zeros((Nbins,1));
H1count=0;
for myiter in range(Niter):
    if( m[2,myiter]==1 ):
        m1 = m[0,myiter];
        k = int((m1-binmin)/Dbin);
        if( (k<0) or (k>=Nbins) ):
            continue;
        H1[k,0] = H1[k,0]+1;
        H1count = H1count+1;

H2 = np.zeros((Nbins,Nbins));
H2count=0;
for myiter in range(Niter):
    if( m[2,myiter]==2 ):
        m1 = m[0,myiter];
        m2 = m[1,myiter];
        k1 = int((m1-binmin)/Dbin);
        k2 = int((m2-binmin)/Dbin);
        if( (k1<0) or (k1>=Nbins) or (k2<0) or (k2>=Nbins) ):
            continue;
        H2[k1,k2] = H2[k1,k2]+1;
        H2count = H2count+1;

print('estimated');
gda_draw(H1,H2);

print("beta: true: %.4f  est: %.4f" % (beta,H1count/(H1count+H2count)));

# true H1 and H2

H1t = np.zeros((Nbins,1));
for i in range(Nbins):
    H1t[i,0] = N1*exp(-0.5*(hbin[i,0]**2)/sm12);

H2t = np.zeros((Nbins,Nbins));
for i in range(Nbins):
    for j in range(Nbins):
        H2t[i,j] = N2*exp(-0.5*(hbin[i,0]**2+hbin[j,0]**2)/sm22);

print('true');
gda_draw(H1t,H2t);


H1n = (H1/(H1count+H2count))/(Dbin);
H2n = (H2/(H1count+H2count))/(Dbin**2);

fig1=plt.figure();
ax1 = plt.subplot(1,1,1);
plt.plot(hbin+Dbin/2,H1n,'k-'); # nis-registered by half a bin
plt.plot(hbin,H1t,'r-');
plt.xlabel('m');
plt.ylabel('p1(m)');
plt.show();
gdapy12_05
estimated
No description has been provided for this image
beta: true: 0.2000  est: 0.1908
true
No description has been provided for this image
No description has been provided for this image
In [31]:
# gdapy12_06
#
# Using the extended Metropolis-Hastibgs algorith to sample the a
# trans-dimensional pdf that simultageously fits data to a sinusoid
# or a bell corve
#
# p(m|d) = [sinusoid, 1 model parameter] + [bell curve, 2 model parameters]

print("gdapy12_06");

# tunable parameter:
# true dimension, can be either 1 or 2
dimtrue = 2;

# x-axis
Nx=41;
xmin = 0.0;
xmax = 1.0;
xmid = (xmax-xmin)/2.0;
Dx = (xmax-xmin)/(Nx-1);
x = gda_cvec(np.linspace(xmin,xmax,Nx));

# model 1, sine wave of amplitude m1
if( dimtrue == 1):
    mtrue = 1.0;
    dtrue = mtrue*np.sin(pi*(x-xmin)/(xmax-xmin));
    smnt2 = 0.2**2;
    dnottrue = mtrue*np.exp(-0.5*np.power(x-xmid,2)/smnt2);
elif( dimtrue == 2): # model 1, gaussian amplitude mp1 and variance mp2
    mptrue = gda_cvec( 1.0, 0.2**2 );
    dtrue = mptrue[0,0]*np.exp(-0.5*np.power(x-xmid,2)/mptrue[1,0]);
    dnottrue = mptrue[0,0]*np.sin(pi*(x-xmin)/(xmax-xmin));

# prior variance of data
sd = 0.2;
sd2 = sd**2;
dobs = dtrue + np.random.normal(loc=0.0,scale=sd,size=(Nx,1));

# plot data
fig1 = plt.figure(1,figsize=(12,5));
ax1=plt.subplot(1,1,1);
plt.axis( [xmin, xmax, 0, 2 ] );
plt.plot(x,dtrue,"k-",lw=3);
plt.plot(x,dnottrue,"g-",lw=3);
plt.plot(x,dobs,"r.",lw=2);
plt.xlabel("z");
plt.ylabel("d");
plt.show();

# dim=1
# p(m|d) = p(d|m) pA(m) / p(d)

# dim=2
# pp(mp1,mp2|d) = pp(d|mp1,mp2) ppA(m1,m2) / p(d)

# note that p(d) cancels for all ratios, and can be ignored
# note that normalization p(d|m) and pp(d|mp1,mp2) is the same
#       and so cancels for all rations and can be ignored

# pA(m) assumed to be normal
# pA(m) = NA exp(-0.5*(m-mA)**2/sA2) with NA = 1/sqrt(2*pi*sA2)

# ppA(mp1,mp2) asumed to be normal
# ppA(mp1,mp2) = NpA exp(-0.5*((mp1-mpA1)**2+(mp2-mpA2)**2)/spA2)
#        with NpA = 1/(2*pi*spA2)

# prior information
# dim=1
mA = 0.9;
sA = 1.0;
sA2 = sA**2;
NA = 1/sqrt(2*pi*sA2);
logNA = log(NA);
# dim=2
mpA = gda_cvec(0.9,0.4**2);
spA = 1.0;
spA2 = spA**2;
NpA = 1/(2*pi*spA2);
logNpA = log(NpA);

# Metropolis Hastings
Niter=1000000;
Nburnt=1000;  # burn in
doswitch=50;  # switch dimensions

# stores links in chain
m = np.zeros((3,Niter)); # m[2,.] 1 or 2 signifies dimension

dimold = 2;  # starting dimension
dimnew = dimold;

# assign first link in chain
if( dimold==1 ):
    mold = np.zeros((1,1));
    mold[0,0] = np.random.uniform(low=0.5,high=1.5);
    dpre = mold*np.sin(pi*(x-xmin)/(xmax-xmin));
    e = dobs-dpre;
    E = np.matmul(e.T,e); E=E[0,0];
    l = mold - mA
    L = l[0,0]**2;
    psi = E/sd2 + L/sA2;
    logpold = logNA-0.5*psi;
    m[0,0] = mold[0,0];
    m[0,1] = 0.0;
    m[0,2] = dimold;
elif( dimold==2 ):
    mold = np.zeros((2,1));
    mold[0,0] = np.random.uniform(low=0.5,high=1.5);
    mold[1,0] = np.random.uniform(low=0.1,high=1.0);
    dpre = mold[0,0]*np.exp(-0.5*np.power(x-xmid,2)/mold[1,0]);
    e = dobs-dpre;
    E = np.matmul(e.T,e); E=E[0,0];
    l = mold - mpA
    L = np.matmul(l.T,l); L=L[0,0];
    psi = E/sd2 + L/spA2;
    logpold = logNpA-0.5*psi;
    m[0,0] = mold[0,0];
    m[0,1] = mold[1,0];
    m[0,2] = dimold;

# random verctors u are normal with these variance and normlizations
su1 = 0.2;
su12 = su1**2;
Nu1 = 1/(sqrt(2.0*pi)*su1);
logNu1 = log(Nu1);

su2 = 0.2;
su22 = su2**2;
Nu2 = 1/(sqrt(2.0*pi)*su2);
logNu2 = log(Nu2);

su3 = 1.0;
su32 = su3**2;
Nu3 = 1/(sqrt(2.0*pi)*su3);
logNu3 = log(Nu3);

# patch to avoid variace going negative
cutoff=0.001;

# Transdimensional Metropolis Hastings 
for myiter in range(1,Niter):
    
    # switch dimensions every doswitch iterations
    if( (myiter%doswitch)==0 ):
        if( dimnew == 1 ):
            dimnew = 2;
        elif ( dimnew == 2 ):
            dimnew = 1;

    # three random vectors u1, u2, u3 with transformation
    # for the dim-1 to dim-2 transformation
    # m1' = [ 1 1 0 0 ] [m1]        (m1' = m1 + u1)
    # m2' = [ 0 0 0 1 ] [u1]        (m2' = u3)
    # u1' = [ 0 1 0 0 ] [u2]        (u1' = u1)
    # u2' = [ 0 0 1 0 ] [u3]        (u2' = u2)
    # which has a Jacobian of unity
    u1 = np.random.normal(loc=0.0,scale=su1); # g1
    u2 = np.random.normal(loc=0.0,scale=su2); # g2
    u3 = np.random.normal(loc=0.0,scale=su3); # g3
       
    if( (dimold==1) and (dimnew==1) ):
        mnew = np.zeros((1,1));
        mnew[0,0] = mold[0,0] + u1;
        dpre = mnew[0,0]*np.sin(pi*(x-xmin)/(xmax-xmin));
        e = dobs-dpre;
        E = np.matmul(e.T,e); E=E[0,0];
        l = mold - mA
        L = l[0,0]**2;
        psi = E/sd2 + L/sA2;
        logpnew = logNA-0.5*psi;
        logalpha = logpnew-logpold;
    elif( (dimold==1) and (dimnew==2) ):
        # pnew = p12'  g1'  g2'
        # pold = p1    g1   g2  g3
        # but since g1=g1' and g2=g2'
        # so pnew/pold = p12' / (p1 g3)
        logg3 = logNu3 - 0.5*(u3**2)/su32;
        mnew = np.zeros((2,1));
        mnew[0,0]= mold[0,0] + u1;
        mnew[1,0]= u3;
        if( mnew[1,0]<cutoff ):
            mnew[1,0]=cutoff;
        dpre = mnew[0,0]*np.exp(-0.5*np.power(x-xmid,2)/mnew[1,0]);
        e = dobs-dpre;
        E = np.matmul(e.T,e); E=E[0,0];
        l = mnew - mpA
        L = np.matmul(l.T,l); L=L[0,0];
        psi = E/sd2 + L/spA2;
        logpnew = logNpA-0.5*psi;
        logalpha = logpnew - (logpold + logg3);
    elif( (dimold==2) and (dimnew==2) ):
        mnew = np.zeros((2,1));
        mnew[0,0]= mold[0,0] + u1;
        mnew[1,0]= mold[1,0] + u2;
        if( mnew[1,0]<cutoff ):
            mnew[1,0]=cutoff;
        dpre = mnew[0,0]*np.exp(-0.5*np.power(x-xmid,2)/mnew[1,0]);
        e = dobs-dpre;
        E = np.matmul(e.T,e); E=E[0,0];
        l = mnew - mpA
        L = np.matmul(l.T,l); L=L[0,0];
        psi = E/sd2 + L/spA2;
        logpnew = logNpA - 0.5*psi;
        logalpha = logpnew - logpold;
    elif( (dimold==2) and (dimnew==1) ):
        mnew = np.zeros((1,1));
        mnew[0,0] = mold[0,0] - u1;
        u3p = mold[1,0];
        dpre = mnew[0,0]*np.sin(pi*(x-xmin)/(xmax-xmin));
        e = dobs-dpre;
        E = np.matmul(e.T,e); E=E[0,0];
        l = mnew - mA
        L = l[0,0]**2;
        psi = E/sd2 + L/sA2;
        logpnew = logNA - 0.5*psi;
        loggp3 = logNu3 - 0.5*(up3**2)/su32;
        logalpha = (loggp3+logpnew)-logpold;
    
    if( logalpha > 0.0):
        dimold = dimnew;
        mold = mnew;
        logpold = logpnew;
    else:
        alpha = exp(logalpha);
        r = np.random.uniform(low=0.0,high=1.0);
        if( alpha > r ):
            dimold = dimnew;
            mold = mnew;
            logpold = logpnew;
        else:
            dimnew = dimold;
            mnew = mold;
            logpnew = logpold;
    if( dimnew == 1 ):
        m[0,myiter]=mnew[0,0];
        m[1,myiter]=0.0;
        m[2,myiter]=dimnew;
    elif( dimnew == 2 ):
        m[0,myiter]=mnew[0,0];
        m[1,myiter]=mnew[1,0];
        m[2,myiter]=dimnew;

# keep only values past burn in
m = np.copy( m[0:3,Nburnt:Niter]);
i, Niter = np.shape(m);

if( dimtrue == 1):
    print("true dim=1, m = %.4f" % (mtrue));
elif( dimtrue == 2): # model 1, gaussian amplitude mp1 and variance mp2
    print("true dim=2, m = %.4f %.4f" % (mptrue[0,0],mptrue[1,0]));
    
Nbins = 41;
binmin = 0.0;
binmax = 2.0;
Dbin = (binmax-binmin)/(Nbins-1);
hbin = gda_cvec(np.linspace(binmin,binmax,Nbins));

H1 = np.zeros((Nbins,1));
H1count=0;
m1sum = 0.0;
for myiter in range(Niter):
    if( m[2,myiter]==1 ):
        m1 = m[0,myiter];
        m1sum = m1sum+m1;
        k = int((m1-binmin)/Dbin);
        if( (k<0) or (k>=Nbins) ):
            continue;
        H1[k,0] = H1[k,0]+1;
        H1count = H1count+1;
if( H1count>0 ):
    m1est = m1sum/H1count;
    print("dim=1 mest = %.4f" % (m1est) );
else:
    print("no dim=1 members");

Nbins = 41;
bin1min = 0.0;
bin1max = 2.0;
bin2min = 0.0;
bin2max = 0.1;
Dbin1 = (bin1max-bin1min)/(Nbins-1);
Dbin2 = (bin2max-bin2min)/(Nbins-1);
hbin1 = gda_cvec(np.linspace(bin1min,bin1max,Nbins));
hbin2 = gda_cvec(np.linspace(bin2min,bin2max,Nbins));
H2 = np.zeros((Nbins,Nbins));
H2count=0;
m1sum = 0.0;
m2sum = 0.0;
for myiter in range(Niter):
    if( m[2,myiter]==2 ):
        m1 = m[0,myiter];
        m2 = m[1,myiter];
        m1sum = m1sum + m1;
        m2sum = m2sum + m2;
        k1 = int((m1-bin1min)/Dbin1);
        k2 = int((m2-bin2min)/Dbin2);
        if( (k1<0) or (k1>=Nbins) or (k2<0) or (k2>=Nbins) ):
            continue;
        H2[k1,k2] = H2[k1,k2]+1;
        H2count = H2count+1;
        
if( H2count>0 ):
    m1est = m1sum/H2count;
    m2est = m2sum/H2count;
    print("dim=2 mest = %.4f %.4f" % (m1est,m2est) );
else:
    print("no dim=2 members");

print("percent of dim=1 %.4f" % (100.0*H1count/(H1count+H2count)));
gda_draw(H1,H2);
gdapy12_06
No description has been provided for this image
true dim=2, m = 1.0000 0.0400
dim=1 mest = 0.9241
dim=2 mest = 1.1128 0.0406
percent of dim=1 11.1011
No description has been provided for this image
In [67]:
# gdapy12_07
#
# Using the extended Metropolis-Hastibgs algorith to sample the
# pdf associated with the Laplace Transform problem, in which the
# model is trans-dimensional, with 1, 2, 4 ... layers.

# The script used a fun transformation and its inverse for layered media
# that presumes the nummber of layers is a power of two.
# forward tranformation:
#   Each old layer is split into N = 2 or 4 or 8 ... new layers,
#      with N=Mnew/Mold and with new layers having same value os old,
#      plus random variation
#   [mp; up] = T [m,u]  where m are layer values, u are random variation
# invere transformation:
#   Each group of N old Layers are grouped into one new layer
#      with new layer having same  value as the mean of the old group,
#      plus random variation
#   [m; u] = Ti [mp,up]  where m are layer values, u are random variation
# definitions Mold, Mnew and N=Mnew/Mold, K=Mnew-Mold, J=K/Mold
# the matrix T is 2*Mnew by 2*Mnew and is broken into several blocks
# Block A, in the upper left is Mnew by Mnew, and is broken into 2 sub-blocks
#     Block AL on the upper left is Mnew by Mold
#         divided into (Mold groups of N rows) and Mold columns
#         each group has a column of 1's, moving one place to right
#     Block AR on the upper right, is Mnew by K
#         it is divided into Mold by Mold sub-blocks that are N by J
#         only the diagonal sub-blocks are populated
#         the top row of each diagonal sub-block is populated by -1's
#         the rest have 1's along their diagonal,starting flush left
# Block B, on the upper right is Mnew by Mnew
#      it is the identity matrix
# Block C, on the lower left is Mnew by Mnew and is unpopulated
# Block D, on the lower left is Mnew by Mnew and is an identify matrix
# for no change of dimension, there is only a forward transformation
# with random noise being added to each model parameter

print("gdapy12_07");

# tunable parameters  (beware: they are carefully tuned!)
# sigma of data = sigmad_factor * datamax
sigmad_factor = 0.005;
sigmad_factor = 0.006;
# sigma of prior model
sigmam = 5.0;
# random variation of models
sigmau = 0.01;

# for error checking
maxE = 1.0E-6;
Ecount=0;

# this codebuilds and checks the transformation
dmax = 4;           # maximum dimension
Mmax = 2**dmax;       # maximum number of layers

# transformation and its invere (only dimension changes considered here)
T = np.zeros((dmax+1,dmax+1,2*Mmax,2*Mmax));
Jac = np.zeros((dmax+1,dmax+1));

for oldpower in range (dmax+1):
    Mold = 2**oldpower;
    for newpower in range(oldpower+1,dmax+1):
        Mnew = 2**newpower;

        N = int(Mnew/Mold);
        K = Mnew-Mold;
        J = int(K/Mold);
        TA = np.zeros((Mnew,Mnew));
        TB = np.identity(Mnew);
        TC = np.zeros((Mnew,Mnew));
        TD = np.identity(Mnew);
        # AL block
        for i in range(Mold):
            for j in range(N):
                TA[i*N+j,i]=1.0;
        # AR block
        for i in range(Mold):
                # populate top row with -1's
                for k in range(J):
                    TA[i*N,Mold+i*J+k]=-1.0
                for k in range(min(N,J)):
                    TA[i*N+k+1,Mold+i*J+k]=1.0
                    
        To = np.concatenate((np.concatenate((TA,TB),axis=1),
                               np.concatenate((TC,TD),axis=1)), axis=0 );
        Toi = la.inv(To);
        
        # T is always 2*Mmax by 2*Mmax
        # initialized to diagional matrix
        T[oldpower,newpower,0:2*Mmax,0:2*Mmax] = np.identity(2*Mmax);
        T[newpower,oldpower,0:2*Mmax,0:2*Mmax] = np.identity(2*Mmax);
        # then set upper-left block
        T[oldpower,newpower,0:2*Mnew,0:2*Mnew] = To;
        T[newpower,oldpower,0:2*Mnew,0:2*Mnew] = Toi;
        
        # Jacobian
        Jac[newpower,oldpower] = abs(la.det(To));
        Jac[oldpower,newpower] = abs(la.det(Toi));
         
        # check crital part of forward transformation, that it preserve value
        m = np.random.uniform(low=0.0,high=1.0,size=(Mold,1));
        u = np.zeros((K+Mnew,1));
        mu = np.concatenate((m,u),axis=0);   
        mup = np.matmul(To,mu);
        E = 0.0;
        for i in range(Mold):
            for j in range(N):
                e = m[i,0] - mup[i*N+j,0];
                E = E+e**2;
        if( E>maxE ):
            Ecount=Ecount+1;
                              
        # check critcal part of inverse transformation, that it preserve mean value
        mp = np.random.uniform(low=0.0,high=1.0,size=(Mnew,1));
        up = np.zeros((Mnew,1));
        mup = np.concatenate((mp,up),axis=0);
        mu = np.matmul(Toi,mup);
        E = 0.0;
        for i in range(Mold):
            s=0.0;
            for j in range(N):
                s = s + mup[i*N+j,0];
            e = mu[i,0] - s/N;
            E = E+e**2;
        if( E>maxE ):
            Ecount=Ecount+1;
            
        # check that the mp-up part of Ti has at least 1 nonzero value per row
        # so that the resulting m's have random variation
        for i in range(Mold):
            Nonzero=0;
            for j in range(Mnew,2*Mnew):
                if (abs(Toi[i,j]) > maxE ):
                    Nonzero=Nonzero+1;
            if (Nonzero==0):
                Ecount=Ecount+1;
                
# no change in dimension
for pof2 in range(dmax+1):
    To = np.zeros((2*Mmax,2*Mmax));
    M = 2**pof2;
    for i in range(M):
        To[i,i] = 1.0;
        To[i,i+M] = 1.0;
    for i in range(M,2*Mmax):
        To[i,i] = 1.0;
    Jo = abs(la.det(To));
    T[pof2,pof2,0:2*Mmax,0:2*Mmax] = To;
    Jac[pof2,pof2] = Jo;

print("Transformations build for dmax=%d Mmax=%d: number of errors: %d" % (dmax, Mmax, Ecount));
print(" ");

dold = 2;
dnew = 2;
Mold = 2**dold;
Mnew = 2**dnew;
print("Example: dold=%d dnew=%d Mold=%d Mnew=%d" % (dold, dnew, Mold, Mnew));
print("Jacobian: %.4f" % (Jac[dold, dnew]));
print("Transformation:");
print( T[dold,dnew,0:2*Mnew,0:2*Mnew] );
print(" ");

dold = 1;
dnew = 2;
Mold = 2**dold;
Mnew = 2**dnew;
print("Example: dold=%d dnew=%d Mold=%d Mnew=%d" % (dold, dnew, Mold, Mnew));
print("Jacobian: %.4f" % (Jac[dold, dnew]));
print("Transformation:")
print( T[dold,dnew,0:2*Mnew,0:2*Mnew] );
print("Inverse Transformation:")
print( T[dnew,dold,0:2*Mnew,0:2*Mnew] );
print("Transformation * Inverse Transformation")
P = np.matmul( T[dold,dnew,0:2*Mnew,0:2*Mnew], T[dnew,dold,0:2*Mnew,0:2*Mnew]);
print(P);
print(' ');

dold = 2;
dnew = 1;
Mold = 2**dold;
Mnew = 2**dnew;
print("Example: dold=%d dnew=%d Mold=%d Mnew=%d" % (dold, dnew, Mold, Mnew));
print("Jacobian: %.4f" % (Jac[dold, dnew]));
print("Transformation:")
print( T[dold,dnew,0:2*Mnew,0:2*Mnew] );
print("Inverse Transformation:")
print( T[dnew,dold,0:2*Mnew,0:2*Mnew] );
print("Transformation * Inverse Transformation")
P = np.matmul( T[dold,dnew,0:2*Mnew,0:2*Mnew], T[dnew,dold,0:2*Mnew,0:2*Mnew]);
print(P);
print(' ');

# I think this easier to understand if I use labels "s" for small "b" for big
# then the vectors is [ms, ub] where ms is length Ms and the other [mb, us]
# where mb is length Mb.  Also, let K = Mb-Ms.

# Case 1: Small to Big (dold<dnew). The you start with [ms, ub] where you supply ub
# by generating 2*Mb-Ms random numbers.  Then you compute [mb, us] via the trasnformation
# [mb, us] = T[dold,dnew,0:2*Mnew,0:2*Mnew] [ms, ub]
# The acceptance function is
#
#         p(mb) Jsb                            gb(ub_K+1)  ...  gb(ub_2*Mmax)   
# Asb  =  --------------------------------     ----------------------------
#         p(ms) gs(us_1) ... gs(us_K)          gs(us_K+1)  ...  gs(us_2*Mmax)
#
# where the right hand fraction is unity, as gs(us_i) is literlly the
# same function as gb(ub_i).
#
# Case 2: Big to Small (dnew<dold). The you start with [mb, us] where you supply us
# by generating Mb random numbers.  Then you compute [mb, us] via the trasnformation
# [mb, us] = T[dnew,dold,0:2*Mnew,0:2*Mnew] [ms, ub]
# The acceptance function is
#
#         p(ms) gs(us_1) ... gs(us_K) Jbs      gs(us_K+1)  ...  gs(us_2*Mmax)   
# Abs =   --------------------------------     ----------------------------
#         p(mb)                                gb(ub_K+1)  ...  gb(ub_2*Mmax)
#
# where the right hand fraction is unity.    Note that the two A's are exactly
# reciprocals, as Jsb=1/Jbs.
#
# Case 3: Same dimension (d=dnew=dold). The you start with [m0, u0] where you supply u0
# by generating M=Mb=Ms random numbers.  Then you compute [mp, up] via the trasnformation
# [mp, up] = T[d,d,0:2*M,0:2*M] [m0, u0]
# The acceptance function is
#
#         p(mp)  
# A =    -----
#         p(m0)
#
# as the Jacobian is unity and all the g's cancel..

# prior for number of dimensions
i = gda_cvec(np.linspace(0,dmax,dmax+1));
cD = 0.2;
PD = np.exp(-cD*i);
PD = PD/np.sum(PD);
logPD = np.log(PD);

# synthetic data
G = np.zeros((Mmax,Mmax));
cmax = Mmax/20.0;
c = gda_cvec(np.linspace(0,cmax,Mmax));
Dz = 1.0;
zmax=Dz*(Mmax-1);
z = gda_cvec(np.linspace(0,zmax,Mmax));
for i in range(Mmax):
    r = np.exp(-c[i,0]*z).T;
    G[i:i+1,0:Mmax] = r;
dtrue = 1;
Mtrue = 2**dtrue;
mtrue = gda_cvec(0.5, 1.5);
print("True dimensions %d" % (Mtrue));
for i in range(Mtrue):
    print("mtrue[%d] = %.2f" % (i,mtrue[i,0]) );
mtrueext = np.concatenate( (mtrue, np.zeros((Mmax-Mtrue,1))), axis=0 );
To = T[dtrue,dmax,0:Mmax,0:Mmax];
mint = np.matmul(To,mtrueext);
N = Mmax;
datatrue = np.matmul(G,mint);
datamax = np.max(np.abs(datatrue));
sigmad = sigmad_factor*datamax
sigmad2 = sigmad**2;
logsigmad = log(sigmad);
logtwopi = log(2.0*pi);
dataobs = datatrue + np.random.normal(loc=0.0,scale=sigmad,size=(Mmax,1));
gda_draw("title G",G,"title m",mint,"=","title d", datatrue);

# prior model
sigmam2 = sigmam**2;
logsigmam = log(sigmam);
mpri = np.ones((Mtrue,1));
mpriext = np.concatenate( (mpri, np.zeros((Mmax-Mtrue,1))), axis=0 );

# random variation of models
sigmau2 = sigmau**2;
logsigmau = log(sigmau);

# starting model
dold = 0;
Mold = 2**dold;
mold = np.ones((Mold,1));
moldext = np.concatenate( (mold, np.zeros((Mmax-Mold,1))), axis=0 );
To = T[dold,dmax,0:Mmax,0:Mmax];
moldint = np.matmul(To,moldext);
dataold = np.matmul(G,moldint);
eold = dataobs - dataold;
Eold = np.matmul(eold.T,eold); Eold=Eold[0,0];
logpEold = -logtwopi-Mmax*logsigmad-0.5*Eold/sigmad**2;
lold = mpriext-moldext ;
Lold = np.matmul(lold.T,lold)/(Mmax/Mold);  Lold=Lold[0,0]; # (Mmax/Mold) accounts for duplicates
logpLold = -logtwopi-Mold*logsigmam-0.5*Lold/sigmam**2;

# Transdimensional Metropolis Hastings 
Nburnt=100000;  # burn in
Niter=1000000+Nburnt; # length of chain
doswitch=50;  # switch model types every doswitch iterations
# stores links in chain
msave = np.zeros((Mmax,Niter));
dsave = np.zeros((Niter,1),dtype=int);
# save starting model
msave[0:Mmax,0:1] = moldint;  # save interpolated models
dsave[0,0] = dold; # save dimension
for myiter in range(1,Niter):

    # switch dimensions every doswitch iterations
    if( (myiter%doswitch)==0 ):
        dnew = np.random.randint(low=0, high=dmax+1);
    else:
        dnew = dold;

    # new model parameters
    Mnew=2**dnew;

    # new model from random perturbatio of old
    uold=np.random.normal(loc=0.0,scale=sigmau,size=(2*Mmax-Mold,1));
    muold = np.concatenate( (mold, uold), axis=0 );
    To = T[dold,dnew,0:2*Mmax,0:2*Mmax] ;
    Jo = Jac[dold,dnew];
    munew = np.matmul(To,muold);
    mnew = munew[0:Mnew,0:1];
    unew = munew[Mnew:2*Mmax,0:1];
    if( dnew>dold ):
        loggstuff = log(Jo);
        K=Mnew-Mold;
        for k in range(K):
            # in numerator, so minus
            loggstuff = loggstuff - logtwopi - logsigmau - 0.5*(uold[k,0]**2)/sigmau2;
    elif( dold>dnew ):
        loggstuff = log(Jo);
        K=Mold-Mnew;
        for k in range(K):
            # in denominator, so positive
            loggstuff = loggstuff + logtwopi + logsigmau + 0.5*(unew[k,0]**2)/sigmau2;
    else:
        loggstuff = 0.0;
    if( dnew < dmax ):
        mnewext = np.concatenate( (mnew, np.zeros((Mmax-Mnew,1))), axis=0 );
        To = T[dnew,dmax,0:Mmax,0:Mmax];
        mnewint = np.matmul(To,mnewext);
    else:
        mnewext = mnew;
        mnewint = mnew;
    datanew = np.matmul(G,mnewint);
    enew = dataobs - datanew;
    Enew = np.matmul(enew.T,enew); Enew=Enew[0,0];
    logpEnew = -logtwopi-Mmax*logsigmad-0.5*Enew/sigmad**2;
    lnew = mpriext-mnewext;
    Lnew = np.matmul(lnew.T,lnew)/(Mmax/Mnew); Lnew=Lnew[0,0]; # (Mmax/Mnew) accounts for duplicates
    logpLnew = -logtwopi-Mnew*logsigmam-0.5*Lnew/sigmam**2;
    logalpha = (logpEnew + logpLnew + loggstuff + logPD[dnew,0]) - (logpEold + logpLold + logPD[dold,0]);

    # error check on reversibility; the loggstuff must be equal and opposite
    # (the reveribility condition must be correct for the method to work)
    RCHECK=0;
    RCHECKITER=151;
    if(RCHECK and (myiter==RCHECKITER)):
        # swap old and new
        loggstuffsave = loggstuff;
        logalphasave = logalpha;
        # save old
        doldsave = dold;
        Moldsave = Mold;
        moldsave = mold;
        uoldsave = uold;
        muoldsave = muold;
        moldextsave = moldext;
        moldintsave = moldint;
        dataoldsave = dataold;
        logpEoldsave = logpEold;
        logpLoldsave = logpLold;
        # old becomes new
        dold = dnew;
        Mold = Mnew;
        mold = mnew;
        uold = unew;
        muold = munew;
        moldext = mnewext;
        moldint = mnewint;
        dataold = datanew;
        logpEold = logpEnew;
        logpLold = logpLnew; 
        # new becomes old
        dnew = doldsave;
        Mnew = Moldsave;
        mnew = moldsave;
        unew = uoldsave;
        munew = muoldsave;
        mnewext = moldextsave;
        mnewint = moldintsave;
        datanew = dataoldsave;
        logpEnew = logpEoldsave;
        logpLnew = logpLoldsave;
        To = T[dold,dnew,0:2*Mmax,0:2*Mmax] ;
        Jo = Jac[dold,dnew];
        if( dnew>dold ):
            loggstuff = log(Jo);
            K=Mnew-Mold;
            for k in range(K):
                # in numerator, so minus
                loggstuff = loggstuff - logtwopi - logsigmau - 0.5*(uold[k,0]**2)/sigmau2;
        elif( dold>dnew ):
            loggstuff = log(Jo);
            K=Mold-Mnew;
            for k in range(K):
                # in denominator, so positive
                loggstuff = loggstuff + logtwopi + logsigmau + 0.5*(unew[k,0]**2)/sigmau2;
        else:
            loggstuff = 0.0;
        if( dnew < dmax ):
            mnewext = np.concatenate( (mnew, np.zeros((Mmax-Mnew,1))), axis=0 );
            To = T[dnew,dmax,0:Mmax,0:Mmax];
            mnewint = np.matmul(To,mnewext);
        else:
            mnewext = mnew;
            mnewint = mnew;
        datanew = np.matmul(G,mnewint);
        enew = dataobs - datanew;
        Enew = np.matmul(enew.T,enew); Enew=Enew[0,0];
        logpEnew = -logtwopi-Mmax*logsigmad-0.5*Enew/sigmad**2;
        lnew = mpriext-mnewext;
        Lnew = np.matmul(lnew.T,lnew)/(Mmax/Mnew); Lnew=Lnew[0,0]# (Mmax/Mnew) accounts for duplicates
        logpLnew = -logtwopi-Mnew*logsigmam-0.5*Lnew/sigmam**2;
        logalpha = (logpEnew + logpLnew + loggstuff + logPD[dnew,0]) - (logpEold + logpLold + logPD[dold,0]);
        print("reversibility check 1: should be zero: ",loggstuff+loggstuffsave);
        print("reversibility check 2: should be zero: ",logalpha+logalphasave);
        # this will only work once, so terminate script
        raise Exception('Exit Early Due to Error Check Being Activted with RCHECK=1');
        
    # standard acceptance test
    if( logalpha > 0.0):
        dold = dnew;
        Mold = Mnew
        mold = mnew;
        uold = unew;
        muold = munew;
        moldint = mnewint;
        Eold = Enew;
        Lold = Lnew;
        logpEold = logpEnew;
        logpLold = logpLnew;
    else:
        alpha = exp(logalpha);
        r = np.random.uniform(low=0.0,high=1.0);
        if( alpha > r ):
            dold = dnew;
            Mold = Mnew
            mold = mnew;
            uold = unew;
            muold = munew;
            moldint = mnewint;
            moldext = mnewint
            Eold = Enew;
            Lold = Lnew;
            logpEold = logpEnew;
            logpLold = logpLnew;
        else:
            dnew = dold;
            Mnew = Mold
            mnew = mold;
            unew = uold;
            munew = muold;
            mnewint = moldint;
            mnewext = moldext;
            Enew = Eold;
            Lnew = Lold;
            logpEnew = logpEold;
            logpLnew = logpLold;

    msave[0:Mmax,myiter:myiter+1] = moldint;  # save interpolated models
    dsave[myiter,0] = dold; # save dimension
    
# delete burn-in
msave = np.copy( msave[0:Mmax,Nburnt:Niter] );
dsave = np.copy( dsave[Nburnt:Niter,0:1] );
Niter, i = np.shape(dsave);

Hdim = np.zeros((dmax+1,1));
for i in range(Niter):
    j = dsave[i,0];
    Hdim[j,0] = Hdim[j,0]+1;
n=gda_cvec(np.linspace(0,dmax,dmax+1));
Hdim = Hdim/np.sum(Hdim);

print("Histogram of model dimensions");
print(Hdim);
print(' ');

dbest = np.argmax(Hdim);
Mbest = 2**dbest;
print("Best-fit dimensions %d" % (Mtrue));
print(' ');

print("mean interpolated model");
mest= np.mean(msave,axis=1);
print(mest);
print(' ');

mbest = np.zeros((Mbest,1));
Nbest=0;
r = int(Mmax/Mbest);
for myiter in range(Niter):
    myd = dsave[myiter,0]
    if( myd != dbest ):
        continue;
    Nbest = Nbest + 1;
    mbest = mbest + msave[0:Mmax:r,myiter:myiter+1];
mbest = mbest/Nbest;
print("mean best-dimensional model");
print(mbest);
print(' ');

# plot data
fig1 = plt.figure(1,figsize=(12,5));
ax1=plt.subplot(1,1,1);
plt.axis( [0, dmax, 0, 1 ] );
plt.plot(n,PD,"k-",lw=3);
plt.plot(n,Hdim,"g-",lw=3);
plt.xlabel("dimensions, n");
plt.ylabel("P(n)");
plt.show();

# histogram of layer probabilities
Nv = 1001;
vmin = 0.0;
vmax = 2.0;
Dv = (vmax-vmin)/(Nv-1)
Hv = np.zeros((Nv,Mmax));
for myiter in range(Niter):
    mint = msave[0:Mmax,myiter:myiter+1];
    for i in range(Mmax):
        j = int((mint[i,0]-vmin)/Dv);
        if( (j<0) or (j>=Nv) ):
            continue;
        Hv[j,i] = Hv[j,i]+1;

# plot model probabilities
fig2 = plt.figure(3,figsize=(12,5));
ax1=plt.subplot(1,1,1);
mycmap=matplotlib.colormaps['jet'];
plt.axis( [0, zmax, vmin, vmax] );
left=0;
right=zmax;
bottom=vmin;
top=vmax;
Hmax = np.max(Hv);
plt.imshow( Hv, cmap=mycmap, vmin=0, vmax=Hmax,
           extent=(left,right,top,bottom), aspect='auto' );
plt.plot( z, mest, "w-", lw=3);
plt.colorbar();
plt.xlabel('z');
plt.ylabel('m(z)');
plt.show();
        
gdapy12_07
Transformations build for dmax=4 Mmax=16: number of errors: 0
 
Example: dold=2 dnew=2 Mold=4 Mnew=4
Jacobian: 1.0000
Transformation:
[[1. 0. 0. 0. 1. 0. 0. 0.]
 [0. 1. 0. 0. 0. 1. 0. 0.]
 [0. 0. 1. 0. 0. 0. 1. 0.]
 [0. 0. 0. 1. 0. 0. 0. 1.]
 [0. 0. 0. 0. 1. 0. 0. 0.]
 [0. 0. 0. 0. 0. 1. 0. 0.]
 [0. 0. 0. 0. 0. 0. 1. 0.]
 [0. 0. 0. 0. 0. 0. 0. 1.]]
 
Example: dold=1 dnew=2 Mold=2 Mnew=4
Jacobian: 0.2500
Transformation:
[[ 1.  0. -1.  0.  1.  0.  0.  0.]
 [ 1.  0.  1.  0.  0.  1.  0.  0.]
 [ 0.  1.  0. -1.  0.  0.  1.  0.]
 [ 0.  1.  0.  1.  0.  0.  0.  1.]
 [ 0.  0.  0.  0.  1.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  1.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  1.  0.]
 [ 0.  0.  0.  0.  0.  0.  0.  1.]]
Inverse Transformation:
[[ 0.5  0.5  0.   0.  -0.5 -0.5  0.  -0. ]
 [ 0.   0.   0.5  0.5  0.   0.  -0.5 -0.5]
 [-0.5  0.5  0.   0.   0.5 -0.5  0.  -0. ]
 [ 0.   0.  -0.5  0.5  0.   0.   0.5 -0.5]
 [ 0.   0.   0.   0.   1.   0.   0.  -0. ]
 [ 0.   0.   0.   0.   0.   1.   0.  -0. ]
 [ 0.   0.   0.   0.   0.   0.   1.  -0. ]
 [ 0.   0.   0.   0.   0.   0.   0.   1. ]]
Transformation * Inverse Transformation
[[1. 0. 0. 0. 0. 0. 0. 0.]
 [0. 1. 0. 0. 0. 0. 0. 0.]
 [0. 0. 1. 0. 0. 0. 0. 0.]
 [0. 0. 0. 1. 0. 0. 0. 0.]
 [0. 0. 0. 0. 1. 0. 0. 0.]
 [0. 0. 0. 0. 0. 1. 0. 0.]
 [0. 0. 0. 0. 0. 0. 1. 0.]
 [0. 0. 0. 0. 0. 0. 0. 1.]]
 
Example: dold=2 dnew=1 Mold=4 Mnew=2
Jacobian: 4.0000
Transformation:
[[ 0.5  0.5  0.   0. ]
 [ 0.   0.   0.5  0.5]
 [-0.5  0.5  0.   0. ]
 [ 0.   0.  -0.5  0.5]]
Inverse Transformation:
[[ 1.  0. -1.  0.]
 [ 1.  0.  1.  0.]
 [ 0.  1.  0. -1.]
 [ 0.  1.  0.  1.]]
Transformation * Inverse Transformation
[[1. 0. 0. 0.]
 [0. 1. 0. 0.]
 [0. 0. 1. 0.]
 [0. 0. 0. 1.]]
 
True dimensions 2
mtrue[0] = 0.50
mtrue[1] = 1.50
No description has been provided for this image
Histogram of model dimensions
[[0.18185]
 [0.27885]
 [0.21645]
 [0.17795]
 [0.1449 ]]
 
Best-fit dimensions 2
 
mean interpolated model
[0.60248536 0.60692625 0.62925189 0.63026336 0.72231027 0.72261172
 0.72640253 0.72679177 1.19199353 1.19241832 1.19342743 1.19330794
 1.2032806  1.20299979 1.20257721 1.20278107]
 
mean best-dimensional model
[[0.57174234]
 [1.34588434]]
 
No description has been provided for this image
No description has been provided for this image
In [66]:
print(dbest);
1
8
275350
[[0.57761439]
 [1.37422419]]
In [ ]: