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
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
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
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
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
beta: true: 0.2000 est: 0.1908 true
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
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
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
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]]
In [66]:
print(dbest);
1 8 275350 [[0.57761439] [1.37422419]]
In [ ]: