edapy_03

In [1]:

    # edapy_03_00 clear all variables and import vatious modules
    %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_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):
            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):
                    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):
                    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));

    # 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 argument.
    # 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;

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 [10]:

    # 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 [12]:

    # 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 [13]:

    # 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 [78]:

    # 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 [79]:

    # 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 [80]:

    # 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]; # 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 [81]:

    # 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 [88]:

    # 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;
    Ddq = 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( Dd*np.sum(p,axis=1) );

    # integrate along d1-axis to get p2, appoximating integral as sum
    p_2 = eda_cvec( Dd*np.sum(p,axis=0) );

    eda_draw( 'title p(d1,d2)', p, 'title p(d1)', p_1, 'title p(d2)', p_2);

In [15]:

    # 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 = Dd*Dd*np.sum(P);
    print('area %.4f' % A);

    eda_draw( 'title p(d1,d2)', P);

    area 1.0000

In [18]:

    # 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 [20]:

    # 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 [26]:

    # 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 [31]:

    # 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 [33]:

    # 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 );
            # 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 ); # 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 ); # 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 [38]:

    # 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));

    1.5 2.5
    cov +: true 0.0500 est 0.0500
    cov -: true -0.0500 est -0.0500
    cov 0: true 0.0000 est -0.0000

In [2]:

    # 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 [ ]: