In [1]:
# edapy_03_00 clear all variables and import vatious modules
# History
# 2022/06/26 - checked with most recent Python & etc. software
# 2023/07/13 - fixed some issues, incl colormap, loglof plt, eda_cvec()
# 2024/05/26 - Patches to accomodate new verions of numpy, scipy, matplotlib
# patched dot product to return scalar value
%reset -f
import os
from datetime import date
from math import exp, pi, sin, sqrt, floor, ceil
import numpy as np
import scipy.linalg as la
from matplotlib import pyplot as plt
import matplotlib.cm as cm
from matplotlib.colors import ListedColormap
# eda_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 eda_draw(*argv):
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;
# 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 eda_cvec(*argv):
t = int;
Nt = 0;
for a in argv:
v = eda_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 = eda_cvec1(a);
N,M = np.shape(v);
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=eda_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=eda_cvec(v1,v2,...) concatenates
# many argiments.
def eda_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("eda_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(v, np.int32):
pass;
elif isinstance(vi,float):
t=float;
elif isinstance(vi,complex):
t=complex;
break;
else:
raise TypeError("eda_cvec: list contains unsupported type %s" % type(vi));
w = np.zeros((r,1),dtype=t);
w[:,0] = v;
return w;
elif isinstance(v, tuple):
r = len(v);
t = int;
for vi in v:
if isinstance(vi,int) or isinstance(v, np.int32):
pass;
elif isinstance(vi,float):
t=float;
elif isinstance(vi,complex):
t=complex;
break;
else:
raise TypeError("eda_cvec: tuple contains unsupported type %s" % type(vi));
w = np.zeros((r,1),dtype=t);
w[:,0] = v;
return w;
else:
raise TypeError("eda_cvec: %s not supported" % type(v));
In [2]:
# edapy_03_01: mode calculated from p(d)
# define an exemplary p(d)
N=15;
Dd = 1.0;
d = eda_cvec( Dd*np.linspace(0, N-1, N) + 0.5*Dd );
p = eda_cvec( [5.00, 15.00, 10.00, 10.00, 10.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00] );
# normalise so integral is unity
p = p / (Dd*np.sum(p));
# mode is d for which p is maximum
# argmax returns the index k at which p is maximum
k = np.argmax(p);
# the mode is the d with index k
dmode = d[k,0];
print('mode', dmode );
# plot the p.d.f with the mode as a dotted line
fig1 = plt.figure(1,figsize=(10,4));
ax1 = plt.subplot(1,1,1);
plt.axis([np.min(d), np.max(d), 0.0, np.max(p)]);
plt.plot(d, p, 'k-');
plt.plot([dmode,dmode],[0,np.max(p)/10],'k:');
plt.xlabel('d');
plt.ylabel('p(d)');
mode 1.5
In [3]:
# edapy_03_02: median calculated from p(d)
# define an exemplary p(d)
N=15;
Dd = 1.0;
d = eda_cvec( Dd*np.linspace(0, N-1, N) + 0.5*Dd );
p = eda_cvec( [5.00, 15.00, 10.00, 10.00, 10.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00] );
# normalise so integral is unity
p = p / (Dd*np.sum(p));
# P is the cumulative probability
# approximate integral as summation
P = eda_cvec( Dd*np.cumsum(p) );
# median is d at the 50% area point of P
# search for first occurrence of P>=0.5
for i in range(N):
if P[i,0] >= 0.5:
dmedian = d[i,0];
break;
print('median', dmedian );
# plot the p.d.f with the median as a dotted line
fig1 = plt.figure(1,figsize=(10,4));
ax1 = plt.subplot(1,1,1);
plt.axis([np.min(d), np.max(d), 0.0, np.max(p)]);
plt.plot(d, p, 'k-');
plt.plot([dmedian,dmedian],[0,np.max(p)/10],'k:');
plt.xlabel('d');
plt.ylabel('p(d)');
median 4.5
In [4]:
# edaps_03_03: mean calculated from p(d)
# define an exemplary p(d)
N=15;
Dd = 1.0;
d = eda_cvec( Dd*np.linspace(0, N-1, N) + 0.5*Dd );
p = eda_cvec( [5.00, 15.00, 10.00, 10.00, 10.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00] );
# normalise so integral is unity
p = p / (Dd*np.sum(p));
# perform the integal for the mean
dmean = Dd*np.sum(np.multiply(p,d));
print('mean %.2f' % dmean );
# plot the p.d.f with the mean as a dotted line
fig1 = plt.figure(1,figsize=(10,4));
ax1 = plt.subplot(1,1,1);
plt.axis([np.min(d), np.max(d), 0.0, np.max(p)]);
plt.plot(d, p, 'k-'); # plot data d1
plt.plot([dmean,dmean],[0,np.max(p)/10],'k:');
plt.xlabel('d');
plt.ylabel('p(d)');
mean 6.30
In [5]:
# edapy_03_04: variance calculated from p(d)
# define an exemplary p(d)
N=15;
Dd = 1.0;
d = eda_cvec( Dd*np.linspace(0, N-1, N) + 0.5*Dd );
p = eda_cvec( [5.00, 15.00, 10.00, 10.00, 10.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00] );
# normalise so integral is unity
p = p / (Dd*np.sum(p));
# mean = (integral) (d-dbar)^2 p(d) dd
# approximate integral as sum
dmean = Dd*np.sum(np.multiply(p,d));
# perform the integral for the variance
q = np.power(d-dmean,2);
sigma2 = Dd*np.sum(np.multiply(q,p));
sigma = sqrt(sigma2);
print('sigma %.4f' % sigma );
# plot the p.d.f with the mean as a solid line and +/- 1 sigma as dotted line
fig1 = plt.figure(1,figsize=(10,4));
ax1 = plt.subplot(1,1,1);
plt.axis([np.min(d), np.max(d), 0.0, np.max(p)]);
plt.plot(d, p, 'k-'); # plot data d1
plt.plot([dmean,dmean],[0,np.max(p)/10],'k-'); # plot zero line
plt.plot([dmean-sigma,dmean-sigma],[0,np.max(p)/10],'k:'); # plot zero line
plt.plot([dmean+sigma,dmean+sigma],[0,np.max(p)/10],'k:'); # plot zero line
plt.xlabel('d');
plt.ylabel('p(d)');
sigma 4.3197
In [6]:
# edapy_03_05: visualation of process of computing variance
# The plot had three vectors, p(d), q(d) and p(d)q(d)
# depicted as grey-shaded coun vectors, side-by-side
# with mean and +/- one sigma levels marked
# set up d axis
Nd = 101;
d = eda_cvec(np.linspace(-30.0,30.0, Nd));
Dd = d[2,0]-d[1,0];
q = np.power(d,2); # q is d-squared
# probability at fixed sigma
db = 0; # mean
sd = 5; # sqrt of variance
sd2 = sd**2; # variance
p = (1/(sqrt(2*pi)*sd)) * np.exp( -0.5*np.power(d-db,2)/sd2 ); # Normal pdf
# variance calcuation (just to check)
# variance = (integral) (d-dbar)^2 p(d) dd = (integral) q p(d) dd with q=(d-dbar)^2
# approximate integral with sum
qp = np.multiply(q,p);
sigma2 =Dd * np.sum(qp);
sigma = sqrt(sigma2);
# find nearest element of array to d=dbar-sigma, d=dbar and d=dbar+sigma
v = np.abs(d-(db-sd));
ilo = np.argmin(v);
v = np.abs(d-db);
imi = np.argmin(v);
v = np.abs(d-(db+sd));
ihi = np.argmin(v);
# make a marker vector with non-zero elements corresponding to mean, mean+/-sigma
mk = np.zeros((Nd,1));
mk[ilo,0] = 1;
mk[imi,0] = 1;
mk[ihi,0] = 1;
eda_draw(mk,'title p', p, 'title q', q, 'title qp', qp, mk);
In [7]:
# edapy_03_06: Normal p.d.f.'s with different means
# list of means
Ndbar=5;
dbar = eda_cvec( np.linspace(10.0,30.0, Ndbar) );
# set up d-axis
Nd = 101;
d = eda_cvec( np.linspace(0,40, Nd) );
# probability at fixed sigma
sd = 5; # sqrt of variacne
sd2 = sd**2; # variance
titles = (); # titles for plot is built incrementally; starts as an empty tuple
p = np.zeros((Nd,Ndbar)); # array, each column of which is a pdf
for i in range(Ndbar): # loop over exemplary means
db = dbar[i,0]; # mean
g = (1/(sqrt(2*pi)*sd)) * np.exp( -0.5*np.power(d-db,2)/sd2 ); # Normal pdf
p[:,i] = g.ravel();
t = "title %.0f" % db; # title of plot
titles = titles + (t,); # append title
eda_draw(titles[0], p[:,0], titles[1], p[:,1], titles[2], p[:,2], titles[3], p[:,3], titles[4], p[:,4]);
In [9]:
# edapy_03_07: Normal p.d.f.'s with different variances
# list of sqrt-variances
Nsigma=5;
sigma = eda_cvec( [2.5, 5, 20, 20, 40] );
Nd = 101;
d = eda_cvec( np.linspace(0.0,40.0, Nd) );
# probability at fixed sigma
p = np.zeros((Nd,Nsigma)); # array, each column of which is a pdf
db = 20; # mean
titles = (); # build list of titles incrementally, starting with an empty tuple
for i in range(Nsigma):
sd = sigma[i,0]; # sqrt of variance
sd2 = sd**2; # variance
g = (1/(sqrt(2*pi)*sd)) * np.exp( -0.5*np.power(d-db,2)/sd2 );
p[:,i]=g.ravel();
t = "title %.0f" % sd
titles = titles + (t,);
eda_draw(titles[0], p[:,0], titles[1], p[:,1], titles[2], p[:,2], titles[3], p[:,3], titles[4], p[:,4]);
In [10]:
# edapy_03_08: transformation of variables m = d^2
# compare original p.d.f. is uniform; tat is, p(d)=1 on 0<d<1
# the transformation of variables is m = d^2
# with transformed p.d.f. is p(m)=0.5*m^(-0.5) on 0<m<1
Nd = 101;
d = eda_cvec( np.linspace(0.01,1.0, Nd) );
Dd = d[2,0] - d[1,0];
# uniform pdf p(d)
pd = np.ones((Nd,1))/Dd;
# set up m-axis
Nm = 101;
m = eda_cvec( np.linspace(0.01,1, Nm) );
Dm = m[2,0] - m[1,0];
# transformed, non-uniform pdf p(m)
pm = 0.5*np.reciprocal(np.sqrt(m));
eda_draw( 'title p(d)', pd, 'title p(m)', pm);
In [11]:
# edapy_03_09: The joint p.d.f. p(d1,d2) and the two univariate
# p.d.f.'s p(d1) and p(d2) formed from it
# set up the d1 and d2 axes
L = 40;
Dd1 = 1.0;
Dd2 = 1.0;
d1 = eda_cvec( Dd1*np.linspace(0,L-1,L) );
d2 = eda_cvec( Dd2*np.linspace(0,L-1,L) );
# 2D normal p.f.d. setup
dbar1 = 10.0;
dbar2 = 20.0;
sigma1 = 5.0;
sigma2 = 5.0;
# 2-d Normal p.d.f.
norm1 = (1/(sqrt(2*pi)*sigma1));
norm2 = (1/(sqrt(2*pi)*sigma2));
p1 = norm1 * np.exp( - np.power(d1-dbar1,2) / (2*sigma1*sigma1) );
p2 = norm2 * np.exp( - np.power(d2-dbar2,2) / (2*sigma2*sigma2) );
# note use of tensor product to construct p(d1,d2)=p(d1)p(d2)
p = np.matmul( p1, p2.T );
# integrate along d2-axis to get p2, appoximating integral as sum
p_1 = eda_cvec( Dd2*np.sum(p,axis=1) );
# integrate along d1-axis to get p2, appoximating integral as sum
p_2 = eda_cvec( Dd1*np.sum(p,axis=0) );
eda_draw( 'title p(d1,d2)', p, 'title p(d1)', p_1, 'title p(d2)', p_2);
In [12]:
# edapy_03_10 The uniform p.d.f. p(d1,d2)=constant
# set up the d1 and d2 axes
Nd1 = 40;
Nd2 = 40;
Dd1 = 1.0;
Dd2 = 1.0;
d1 = eda_cvec( Dd1*np.linspace(0,Nd1-1,Nd1) );
d2 = eda_cvec( Dd2*np.linspace(0,Nd2-1,Nd2) );
# make a uniform p.d.f. and normalize it
P = np.ones((Nd1, Nd2));
# calculate the area under the p.d.f.
norm = Dd1*Dd2*np.sum(P);
# normalize the p.d.f. by dividing by the area
P = P/norm;
# area of pdf (should be 1)
A = Dd1*Dd2*np.sum(P);
print('area %.4f' % A);
eda_draw( 'title p(d1,d2)', P);
area 1.0000
In [13]:
# edapy_03_11: construct a joint Normal p.f.d.
# and compute its mean and variance
# set up d1 axis
Nd1 = 101;
d1 = eda_cvec( np.linspace(0.0,5.0, Nd1) );
Dd1 = d1[1,0] - d1[0,0];
# set up d2 axis
Nd2 = 101;
d2 = eda_cvec( np.linspace(0,5,Nd2) );
Dd2 = d2[1,0] - d2[0,0];
# 2D normal p.f.d. setup
dbar1 = 1.5;
dbar2 = 2.5;
sigma1 = 0.25;
sigma2 = 0.25;
# since d1 and d2 are uncorrelated, the 2-d Normal
# p.d.f. p(d1,d2) is product of p(d1) and p(d2)
norm1=(1/(sqrt(2*pi)*sigma1)); # normalization of p(d1)
norm2=(1/(sqrt(2*pi)*sigma2)); # normalization of p(d2)
# p(d1)
p1=norm1*np.exp(-np.power(d1-dbar1,2)/(2*sigma1*sigma1));
# p(d2)
p2=norm2*np.exp(-np.power(d2-dbar2,2)/(2*sigma2*sigma2));
P=np.matmul(p1,p2.T); # p(d1,d2) = p(d1) p(d2)
eda_draw(P);
# integrate along d2-axis to get p1
p_1 = eda_cvec( Dd2*np.sum(P,axis=1) );
# integrate along d1-axis to get p2
p_2 = eda_cvec( Dd1*np.sum(P,axis=0) );
# areas (should all be 1) (approximate integrals as sums)
A = Dd1*Dd2*np.sum(P);
A1 = Dd1*np.sum(p_1);
A2 = Dd2*np.sum(p_2);
# vectors of ones
w1 = np.ones((Nd1,1));
w2 = np.ones((Nd2,1));
# mean along 1 axis, computed from p
D1 = np.matmul( d1, w2.T );
dbar_1 = Dd1*Dd2*np.sum( np.multiply(D1, P) );
# mean along 1 axis, computed from p_1
dbar_1a = Dd1*np.sum( np.multiply(d1, p_1) );
# mean along 2 axis, computed from p
D2 = np.matmul( w1, d2.T );
dbar_2 = Dd1*Dd2*np.sum( np.multiply(D2, P) );
# mean along 2 axis, computed from p_2
dbar_2a = Dd2*np.sum( np.multiply(d2, p_2) );
# eda03_13: variances
# variance along 1 axis, computed from p(d1,d2)
# use a tensor product to create a 2D version of the quantity (d1-d1bar)**2
D1s = np.matmul( np.power(d1-dbar1,2), w2.T );
# compute the variance with the usual integral, approximated as a sum
sigma_1_squared = Dd1*Dd2*np.sum( np.multiply(D1s, P) );
sigma_1 = sqrt(sigma_1_squared);
# variance along 1 axis, computed from p(d1) via the usual integral, approximated as a sum
sigma_1a_squared = Dd1*np.sum( np.multiply( np.power(d1-dbar_1,2), p_1) );
sigma_1a = sqrt(sigma_1_squared);
# variance along 2 axis, computed from p(d1,d2)
# use a tensor product to create a 2D version of the quantity (d2-d2bar)**2
D2s = np.matmul( w1, np.power(d2-dbar2,2).T );
# compute the variance with the usual integral, approximated as a sum
sigma2_squared = Dd1*Dd2*np.sum( np.multiply(D2s, P) );
sigma_2 = sqrt(sigma2_squared);
# variance along 2 axis, computed from p(d2) via the usual integral, approximated as a sum
sigma_2a_squared = Dd2*np.sum( np.multiply( np.power(d2-dbar_2,2), p_2) );
sigma_2a = sqrt(sigma_2a_squared);
print('area %.4f %.4f %.4f' % (A, A1, A2));
print('dbar1 %.4f %.4f %.4f' % (dbar1, dbar_1, dbar_1a));
print('dbar2 %.4f %.4f %.4f' % (dbar2, dbar_2, dbar_2a));
print('sigma1 %.4f %.4f %.4f' % (sigma1, sigma_1, sigma_1a));
print('sigma2 %.4f %.4f %.4f' % (sigma2, sigma_2, sigma_2a));
area 1.0000 1.0000 1.0000 dbar1 1.5000 1.5000 1.5000 dbar2 2.5000 2.5000 2.5000 sigma1 0.2500 0.2500 0.2500 sigma2 0.2500 0.2500 0.2500
In [14]:
# edapy_03_12: example of a conditional p.d.f.
# set up d1-axis
Nd1 = 101;
d1 = eda_cvec( np.linspace(0.0,5.0, Nd1) );
Dd1 = d1[1,0] - d1[0,0];
# set up d2-axis
Nd2 = 101;
d2 = eda_cvec( np.linspace(0.0,5.0,Nd2) );
Dd2 = d2[1,0] - d2[0,0];
# 2D normal p.f.d. setup
dbar1 = 1.5;
dbar2 = 3.5;
sigma1 = 0.75;
sigma2 = 0.75;
# since d1 and d2 are uncorrelated, the 2-d Normal
# p.d.f. p(d1,d2) is product of p(d1) and p(d2)
norm1 = (1/(sqrt(2*pi)*sigma1)); # p(d1) normalization
norm2 = (1/(sqrt(2*pi)*sigma2)); # p(d2) normalization
p1 = norm1 * np.exp( - np.power(d1-dbar1,2) / (2*sigma1*sigma1) ); # p(d1)
p2 = norm2 * np.exp( - np.power(d2-dbar2,2) / (2*sigma2*sigma2) ); # p(d2)
P = np.matmul( p1, p2.T ); # use tensor product to make p(d1,d2) = p(d1)p(d2)
# integrate p(d1,d2) along d2-axis to get p2
p_1 = eda_cvec( Dd2*np.sum(P,axis=1) );
# integrate p(d1,d2) along d1-axis to get p2
p_2 = eda_cvec( Dd1*np.sum(P,axis=0) );
# conditional p.d.f.: p(d1|d2) = p(d1,d2)/p2
p1g2 = np.divide( P, np.matmul( np.ones((Nd1,1)) ,p_2.T));
# conditional p.d.f.: p(d2|d1) = p(d1,d2)/p2
p2g1 = np.divide( P, np.matmul(p_1,np.ones((Nd2,1)).T));
eda_draw( 'title p(d1,d2)', P, 'title p(d1|d2)', p1g2, 'title p(d2|d1)', p2g1);
In [15]:
# edama_03_13: covariace calculated from p.d.f. p(d1,d2)
# set up d1-axis
Nd1 = 101;
d1 = eda_cvec( np.linspace(-2.5,2.5, Nd1) );
Dd1 = d1[1,0] - d1[0,0];
# set up d2-axis
Nd2 = 101;
d2 = eda_cvec( np.linspace(-2.5,2.5,Nd2) );
Dd2 = d2[1,0] - d2[0,0];
# create a highly and positively correlated p.d.f.
P = np.zeros((Nd1,Nd2));
for i in range(Nd1):
for j in range(Nd2):
d_lo = -1.0+d1[i,0];
d_hi = 1.0+d1[i,0];
r2 = 2.0 + d1[i,0]**2 + d2[j,0]**2;
P[i,j] = (d2[j,0]>=d_lo)*(d2[j,0]<=d_hi)/r2;
norm = Dd1*Dd2*np.sum(P);
P = P/norm;
# compute means
d1w = np.matmul(d1,np.ones((1,Nd2)));
wd2 = np.matmul(np.ones((Nd1,1)),d2.T);
d1bar = Dd1*Dd2*np.sum(np.multiply(d1w,P));
d2bar = Dd1*Dd2*np.sum(np.multiply(wd2,P));
# make the alternating sign function
S = np.matmul((d1-d1bar),(d2-d2bar).T);
# perform covariance integral
sigma12 = (Dd1*Dd2)*np.sum(np.multiply(S,P));
print("d1bar %.3f d2bar %.3f covariance %.3f" % (d1bar, d2bar, sigma12) );
eda_draw('title s(d1,d2)', S, 'title p(d1,d2)', P, 'title s(d1,d2) p(d1,d2)', np.multiply(S,P));
d1bar 0.002 d2bar 0.002 covariance 0.910
In [16]:
# edapy_03_14: Constructing Normal p.d.f. without for loop
# set up d1-axis
Nd1 = 101;
d1 = eda_cvec( np.linspace(0.0,5.0, Nd1) );
Dd1 = d1[1,0] - d1[0,0];
# set up d2-axis
Nd2 = 101;
d2 = eda_cvec( np.linspace(0.0,5.0,Nd2) );
Dd2 = d2[1,0] - d2[0,0];
# 2D normal p.f.d. setup, positively correlated
dbar1 = 1.5;
dbar2 = 2.5;
sigma1 = 0.25;
sigma2 = 0.25;
covar = 0.05; # positive covariance
covpos = covar;
C = np.array( [ [sigma1**2, covar], [covar, sigma2**2]] ); # covariance matrix
Ci = la.inv(C); # its inverse
sqrtdet = sqrt(la.det(C)); # square root of its determinant
norm=(1/(2*pi*sqrtdet)); # normalization
# vectors of ones
w1 = np.ones((Nd1,1));
w2 = np.ones((Nd2,1));
# Normal p.d.f., positive covariance
dd1 = d1-dbar1; # d1 minus its mean
dd2 = d2-dbar2; # d2 minus its mean
dd12 = np.power(dd1,2); # (d1 minus its mean) squared
dd22 = np.power(dd2,2); # (d2 minus its mean) squared
# 2D Normal pdf p(d1,d2) using tensor products
# my sense its that this expression is way to complicate to easy check
# and that an impementation using for loops, though computationally inefficient,
# would be clearer (see 03_16, below)
P=norm*np.exp(-0.5*Ci[0,0]*np.matmul(dd12,w2.T)
-0.5*Ci[1,1]*np.matmul(w1,dd22.T)
-Ci[0,1]*np.matmul(dd1,dd2.T));
pp=P;
# 2D normal p.f.d. setup, negatively correlated
dbar1 = 1.5;
dbar2 = 2.5;
sigma1 = 0.25;
sigma2 = 0.25;
covar = -0.05; # positive covariance
covpos = covar;
C = np.array( [ [sigma1**2, covar], [covar, sigma2**2]] ); # covariance matrix
Ci = la.inv(C); # its inverse
sqrtdet = sqrt(la.det(C)); # square root of its determinant
norm=(1/(2*pi*sqrtdet)); # normalization
# Normal p.d.f., negative covariance
dd1 = d1-dbar1; # d1 minus its mean
dd2 = d2-dbar2; # d2 minus its mean
dd12 = np.power(dd1,2); # (d1 minus its mean) squared
dd22 = np.power(dd2,2); # (d2 minus its mean) squared
# 2D Normal pdf p(d1,d2) using tensor products
# my sense its that this expression is way to complicate to easy check
# and that an impementation using for loops, though computationally inefficient,
# would be clearer (see 03_16, below)
pn=norm*np.exp(-0.5*Ci[0,0]*np.matmul(dd12,w2.T)
-0.5*Ci[1,1]*np.matmul(w1,dd22.T)
-Ci[0,1]*np.matmul(dd1,dd2.T));
# 2D normal p.f.d. setup, uncorrelated
dbar1 = 1.5;
dbar2 = 2.5;
sigma1 = 0.25;
sigma2 = 0.25;
covar = 0.0; # zero covariance
covpos = covar;
C = np.array( [ [sigma1**2, covar], [covar, sigma2**2]] ); # covariance matrix
Ci = la.inv(C); # its inverse
sqrtdet = sqrt(la.det(C)); # square root of its determinant
norm=(1/(2*pi*sqrtdet)); # normalization
# Normal p.d.f., uncorrelated
dd1 = d1-dbar1; # d1 minus its mean
dd2 = d2-dbar2; # d2 minus its mean
dd12 = np.power(dd1,2); # (d1 minus its mean) squared
dd22 = np.power(dd2,2); # (d2 minus its mean) squared
# 2D Normal pdf p(d1,d2) using tensor products
# my sense its that this expression is way to complicate to easy check
# and that an impementation using for loops, though computationally inefficient,
# would be clearer (see 03_16, below)
pu=norm*np.exp(-0.5*Ci[0,0]*np.matmul(dd12,w2.T)
-0.5*Ci[1,1]*np.matmul(w1,dd22.T)
-Ci[0,1]*np.matmul(dd1,dd2.T));
eda_draw( 'title cov +', pp, 'title cov -', pn, 'title uncor', pu);
In [18]:
# edapy_03_15: Constructing Normal p.d.f. with for loops
# set up d1-axis
Nd1 = 101;
d1 = eda_cvec( np.linspace(0.0,5.0, Nd1) );
Dd1 = d1[1,0] - d1[0,0];
# set up d2-axis
Nd2 = 101;
d2 = eda_cvec( np.linspace(0.0,5.0,Nd2) );
Dd2 = d2[1,0] - d2[0,0];
# Part 1: p(d1,d2) with positive covariance
# 2D normal p.f.d. setup
dbar1 = 1.5;
dbar2 = 2.5;
sigma1 = 0.25;
sigma2 = 0.25;
covar = 0.05;
covpos = covar;
Cd = [ [sigma1**2, covar], [covar, sigma2**2]];
Cdinv = la.inv(Cd);
sqrtdet = sqrt(la.det(Cd));
# 2-d Normal p.d.f.
P = np.zeros((Nd1,Nd2));
dd = np.zeros((2,1));
for i in range(Nd1):
for j in range(Nd2):
# vector of data minus their means
dd[:,0] = [ d1[i,0] - dbar1, d2[j,0] - dbar2 ];
# ddT inv(C) dd
E = np.matmul( np.matmul(dd.T,Cdinv), dd ); E=E[0,0];
# p(d1,d2)
P[i,j] = (1/(2*pi*sqrtdet)) * exp( -0.5*E );
ppos = np.copy(P);
# part 2: p(d1,d2) with negative covariance
# 2D normal p.f.d. setup
dbar1 = 1.5;
dbar2 = 2.5;
sigma1 = 0.25;
sigma2 = 0.25;
covar = -0.05;
covneg = covar;
Cd = [ [sigma1**2, covar], [covar, sigma2**2]]; # covariance matrix
Cdinv = la.inv(Cd); # its inverse
sqrtdet = sqrt(la.det(Cd)); # square root of it determinant
P = np.zeros((Nd1,Nd2));
dd = np.zeros((2,1));
for i in range(Nd1):
for j in range(Nd2):
dd[:,0] = [ d1[i,0] - dbar1, d2[j,0] - dbar2 ]; # vector of data minus their means
E = np.matmul( np.matmul(dd.T,Cdinv), dd ); E=E[0,0]; # ddT inv(C) dd
P[i,j] = (1/(2*pi*sqrtdet)) * exp( -0.5*E ); # p(d1,d2)
pneg = np.copy(P);
# Part 3: p(d1,d2) with zero covariance
# 2D normal p.f.d. setup
dbar1 = 1.5;
dbar2 = 2.5;
sigma1 = 0.25;
sigma2 = 0.25;
covar = 0.0;
covzero = covar;
Cd = [ [sigma1**2, covar], [covar, sigma2**2]]; # covariance matrix
Cdinv = la.inv(Cd); # its inverse
sqrtdet = sqrt(la.det(Cd)); # square root of it determinant
# 2-d Normal p.d.f.
P = np.zeros((Nd1,Nd2));
dd = np.zeros((2,1));
for i in range(Nd1):
for j in range(Nd2):
dd[:,0] = [ d1[i,0] - dbar1, d2[j,0] - dbar2 ]; # vector of data minus their means
E = np.matmul( np.matmul(dd.T,Cdinv), dd ); E=E[0,0]; # ddT inv(C) dd
P[i,j] = (1/(2*pi*sqrtdet)) * exp( -0.5*E ); # p(d1,d2)
pzero=np.copy(P);
eda_draw( 'title cov +', ppos, 'title cov -', pneg, 'title cov 0', pzero);
In [17]:
# edapy_03_16: compute covariance of p(d1,d2)
# warning: you must have run the previous cell before this one will work
S = np.matmul( (d1-dbar1), (d2-dbar2).T ); # use tensor product to make (d1-dbar1)(d2-dbar2)
covpos_est = Dd1*Dd2*np.sum( np.multiply(S, ppos) ); # covariance intergal, + covariance case
covneg_est = Dd1*Dd2*np.sum( np.multiply(S, pneg) ); # covariance intergal, - covariance case
covzero_est = Dd1*Dd2*np.sum( np.multiply(S, pzero) ); # covariance intergal, 0 covariance case
print( 'cov +: true %.4f est %.4f' % (covpos, covpos_est) );
print( 'cov -: true %.4f est %.4f' % (covneg, covneg_est) );
print( 'cov 0: true %.4f est %.4f' % (covzero, covzero_est) );
eda_draw( 'title p(d1,d2)', ppos, 'title s(d1,d2)', S, 'title sp', np.multiply(S,P));
cov +: true 0.0500 est 0.0500 cov -: true -0.0500 est -0.0500 cov 0: true 0.0000 est -0.0000
In [19]:
# edapy_03_17: error propagation from data to model parameters
# set up d1-axis
Nd1 = 101;
d1 = eda_cvec( np.linspace(0.0,5.0, Nd1) );
Dd1 = d1[1,0] - d1[0,0];
# set up d2-axis
Nd2 = 101;
d2 = eda_cvec( np.linspace(0.0,5.0,Nd2) );
Dd2 = d2[1,0] - d2[0,0];
# 2D normal p.f.d. setup, uncorrelated
dbar1 = 1.0;
dbar2 = 4.0;
sigma1 = 0.25;
sigma2 = 0.25;
covar = 0.0; # uncorrelated
Cd = np.array( [ [sigma1**2, covar], [covar, sigma2**2]] ); # covariance matrix
Cdi = la.inv(Cd); # its inverse
sqrtdet = sqrt(la.det(Cd)); # square root of its determinant
norm=(1/(2*pi*sqrtdet)); # normalization
# vectors of ones
w1 = np.ones((Nd1,1));
w2 = np.ones((Nd2,1));
# Normal p.d.f., positive covariance
dd1 = d1-dbar1; # d1 minus its mean
dd2 = d2-dbar2; # d2 minus its mean
dd12 = np.power(dd1,2); # (d1 minus its mean) squared
dd22 = np.power(dd2,2); # (d2 minus its mean) squared
pd=norm*np.exp(-0.5*Cdi[0,0]*np.matmul(dd12,w2.T)-0.5*Cdi[1,1]*np.matmul(w1,dd22.T)-Cdi[0,1]*np.matmul(dd1,dd2.T));
# model is related to data by m=Md
M = np.array( [ [-1, 1], [-2, 1] ] );
# error propagation
mbar = np.matmul(M, eda_cvec( dbar1, dbar2) );
mbar1 = mbar[0,0];
mbar2 = mbar[1,0];
Cm = np.matmul(M,np.matmul(Cd,M.transpose()));
Cmi = la.inv(Cm);
sqrtdet = sqrt(la.det(Cm)); # square root of its determinant
norm=(1/(2*pi*sqrtdet)); # normalization
# print results
print('mbar1=%.2f mbar2=%.2f' % (mbar1,mbar2) );
print('varm1=%.2f varm2=%.2f covm=%.2f' % (Cm[0,0], Cm[1,1], Cm[0,1]) );
# set up m1-axis
Nm1 = 101;
m1 = eda_cvec( np.linspace(0.0,5.0, Nm1) );
Dm1 = m1[1,0] - m1[0,0];
# set up m2-axis
Nm2 = 101;
m2 = eda_cvec( np.linspace(0.0,5.0,Nm2) );
Dm2 = m2[1,0] - m2[0,0];
# vectors of ones
w1 = np.ones((Nm1,1));
w2 = np.ones((Nm2,1));
# Normal p.d.f.
dm1 = m1-mbar1; # d1 minus its mean
dm2 = m2-mbar2; # d2 minus its mean
dm12 = np.power(dm1,2); # (d1 minus its mean) squared
dm22 = np.power(dm2,2); # (d2 minus its mean) squared
pm=norm*np.exp(-0.5*Cmi[0,0]*np.matmul(dm12,w2.T)-0.5*Cmi[1,1]*np.matmul(w1,dm22.T)-Cmi[0,1]*np.matmul(dm1,dm2.T));
eda_draw( 'title p(d)', pd, 'title p(m)', pm);
mbar1=3.00 mbar2=2.00 varm1=0.12 varm2=0.31 covm=0.19
In [ ]: