# Graphical model experiment¶

The 8th KIAS CAC Summer School

Machine Learning Practice III - Graphical models example

Author: Yung-Kyun Noh, Ph.D.

#### We show the application of graphical model for Gaussian data. We consider each class has the casuality graph,¶

## Data generation¶

\begin{eqnarray} \left(\begin{array}{c} X_1 \\ X_2 \end{array}\right) \sim \mathcal{N}(\mu_{12}, \Sigma_{12}), \quad \left(\begin{array}{c} X_7 \\ X_8 \\ X_9 \\ X_{10} \end{array}\right) \sim \mathcal{N}(\mu_{7890}, \Sigma_{7890}), \end{eqnarray}\begin{eqnarray} \left(\begin{array}{c} X_3 \\ X_4 \\ X_5 \end{array}\right) &=& W_{345,12} \left(\begin{array}{c} X_1 \\ X_2 \end{array}\right) + \gamma\left(\begin{array}{c} \epsilon_3 \\ \epsilon_4 \\ \epsilon_5 \end{array}\right) \\ X_6 &=& W_{6,345} \left(\begin{array}{c} X_3 \\ X_4 \\ X_5 \end{array}\right) + W_{6,7890} \left(\begin{array}{c} X_7 \\ X_8 \\ X_9 \\ X_{10} \end{array}\right) + \gamma \epsilon_6 \end{eqnarray}\begin{eqnarray} \epsilon_3, \epsilon_4, \epsilon_5, \epsilon_6 \sim \mathcal{N}(0, 1) \end{eqnarray}

## Parameter estimation¶

\begin{eqnarray} P(X_1,\ldots,X_{10}) &=& P(X_1,X_2)P(X_3,X_4,X_5|X_1,X_2)P(X_7,X_8,X_9,X_{10})P(X_6|X_3,X_4,X_5,X_7,\ldots,X_{10})\\ &=& P(X_1,X_2)\frac{P(X_1,X_2,X_3,X_4,X_5)}{P(X_1,X_2)}P(X_7,X_8,X_9,X_{10})\frac{P(X_3,\ldots,X_{10})}{P(X_3,X_4,X_5,X_7,\ldots,X_{10})} \\ &=& P(X_1,X_2)\frac{P(X_1,X_2,X_3,X_4,X_5)}{P(X_1,X_2)}P(X_7,X_8,X_9,X_{10})\frac{P(X_3,\ldots,X_{10})}{P(X_3,X_4,X_5)P(X_7,\ldots,X_{10})} \\ &=& P(X_1,X_2,X_3,X_4,X_5)\frac{P(X_3,X_4,X_5,X_6,X_7,X_8,X_9,X_{10})}{P(X_3,X_4,X_5)} \end{eqnarray}
In [1]:

    import numpy as np

    %pylab inline

    Populating the interactive namespace from numpy and matplotlib

### Make true density parameters¶

In [2]:

    dim = 10;

    # True density
    trueMean_1 = np.zeros(dim)
    trueMean_2 = np.concatenate(([1],np.zeros(dim-1)))
    print trueMean_1
    print trueMean_2

    # covariance matrix of X1, X2 (root nodes)
    curDim = 2
    trueCov12_1 = np.random.randn(curDim,curDim)
    trueCov12_1 = np.dot(trueCov12_1.T, trueCov12_1)
    trueCov12_2 = np.random.randn(curDim,curDim)
    trueCov12_2 = trueCov12_1 + .001*np.dot(trueCov12_2.T, trueCov12_2)

    print trueCov12_1
    print trueCov12_2

    # covariance matrix of X7, X8, X9, X10 (root nodes)
    curDim = 4;
    trueCov7890_1 = np.random.randn(curDim,curDim)
    trueCov7890_1 = np.dot(trueCov7890_1.T, trueCov7890_1)
    trueCov7890_2 = np.random.randn(curDim,curDim)
    trueCov7890_2 = trueCov7890_1 + .1*np.dot(trueCov7890_2.T, trueCov7890_2)

    print trueCov7890_1
    print trueCov7890_2

    [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
    [ 1.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
    [[ 0.66990231  1.05063391]
     [ 1.05063391  2.23574801]]
    [[ 0.67841322  1.04711639]
     [ 1.04711639  2.2385393 ]]
    [[ 2.95261466  0.89918369 -0.71586731  0.00907099]
     [ 0.89918369  2.37577021  0.96099359 -0.29608912]
     [-0.71586731  0.96099359  1.18659954  1.00897774]
     [ 0.00907099 -0.29608912  1.00897774  4.05881863]]
    [[ 3.16863259  0.94613626 -0.65645078  0.05682842]
     [ 0.94613626  2.5625525   1.17490649 -0.49851749]
     [-0.65645078  1.17490649  1.48813987  0.73875608]
     [ 0.05682842 -0.49851749  0.73875608  4.36524807]]

### Generate data according to the true parameters¶

In [3]:

    # data generation
    datanum = 11
    tstdatanum = 1000

In [4]:

    # root node X1, X2
    curDim = 2

    X12_1 = np.random.multivariate_normal(trueMean_1[0:2], trueCov12_1, datanum)
    tstX12_1 = np.random.multivariate_normal(trueMean_1[0:2], trueCov12_1, tstdatanum)

    X12_2 = np.random.multivariate_normal(trueMean_2[0:2], trueCov12_2, datanum)
    tstX12_2 = np.random.multivariate_normal(trueMean_2[0:2], trueCov12_2, tstdatanum)

    # root node X7, X8, X9, X10
    curDim = 4
    X7890_1 = np.random.multivariate_normal(trueMean_1[6:10], trueCov7890_1, datanum)
    tstX7890_1 = np.random.multivariate_normal(trueMean_1[6:10], trueCov7890_1, tstdatanum)
    X7890_2 = np.random.multivariate_normal(trueMean_2[6:10], trueCov7890_2, datanum)
    tstX7890_2 = np.random.multivariate_normal(trueMean_2[6:10], trueCov7890_2, tstdatanum)

In [5]:

    # Non-root node X3, X4, X5
    curDim = 3
    curW = np.random.randn(2,curDim)
    #print curW
    X345_1 = np.dot(X12_1,curW) + .1*np.random.randn(datanum,curDim) # non-isotropic noise can be added
    #print X345_1
    tstX345_1 = np.dot(tstX12_1,curW) + .1*np.random.randn(tstdatanum,curDim) # non-isotropic noise can be added
    curW = curW + .1*np.random.randn(2,curDim)
    X345_2 = np.dot(X12_2,curW) + .1*np.random.randn(datanum,curDim)   # non-isotropic noise can be added
    tstX345_2 = np.dot(tstX12_2,curW) + .1*np.random.randn(tstdatanum,curDim)  # non-isotropic noise can be added

    # Non-root node X6
    curDim = 1
    curW = np.random.randn(7,curDim)
    #print np.hstack([X345_1,X7890_1])
    X6_1 = np.dot(np.hstack([X345_1,X7890_1]),curW) + .1*np.random.randn(datanum,curDim) # non-isotropic noise can be added
    tstX6_1 = np.dot(np.hstack([tstX345_1,tstX7890_1]),curW) + .1*np.random.randn(tstdatanum,curDim) # non-isotropic noise can be added
    curW = curW + .1*np.random.randn(7,curDim)
    X6_2 = np.dot(np.hstack([X345_2,X7890_2]),curW) + .1*np.random.randn(datanum,curDim) # non-isotropic noise can be added
    tstX6_2 = np.dot(np.hstack([tstX345_2,tstX7890_2]),curW) + .1*np.random.randn(tstdatanum,curDim) # non-isotropic noise can be added

In [6]:

    X1 = np.hstack([X12_1,X345_1,X6_1,X7890_1])
    X2 = np.hstack([X12_2,X345_2,X6_2,X7890_2])
    tstX1 = np.hstack([tstX12_1,tstX345_1,tstX6_1,tstX7890_1])
    tstX2 = np.hstack([tstX12_2,tstX345_2,tstX6_2,tstX7890_2])

In [7]:

    # Classification with joint density function
    # X1, X2 \in R^{dim x datanum}
    estM1 = sum(X1,axis=0)/datanum
    estM2 = sum(X2,axis=0)/datanum

    estCov1 = np.dot(X1.T,X1)/datanum - np.outer(estM1, estM1)
    estCov2 = np.dot(X2.T,X2)/datanum - np.outer(estM2, estM2)

In [8]:

    def getLogGaussian(X, mu, Sig):
        ''' evalute log probability of Gaussian distribution
        given the Gaussian distribution model (mean and covariance).'''
        # X: datanum x dim
        nData, nDim = X.shape
        detSig = np.linalg.det(Sig)
        invSig = np.linalg.inv(Sig)
        logPs = np.zeros(nData)

        const = - nDim / (2.) * np.log(2*np.pi) - 1 / (2.) * np.log(detSig) # sacalar
        for i in range(nData):
            x = X[i,:]
            expval = - 1/2. * (x - mu).reshape(1,-1).dot(invSig).dot((x - mu).reshape(-1,1)) # 2d array
            logPs[i] =  const + expval
        return logPs

## Full covariance matrix¶

In [9]:

    probC1tst1 = getLogGaussian(tstX1, estM1, estCov1)   # array
    probC2tst1 = getLogGaussian(tstX1, estM2, estCov2)
    probC1tst2 = getLogGaussian(tstX2, estM1, estCov1)
    probC2tst2 = getLogGaussian(tstX2, estM2, estCov2)
    print sum(probC1tst1 - probC2tst1 > 0)
    print sum(probC2tst2 - probC1tst2 > 0)
    accuracy = (sum(probC1tst1 - probC2tst1 > 0) +
        sum(probC2tst2 - probC1tst2 > 0))/(2.*tstdatanum)
    print accuracy

    658
    934
    0.796

## From graphical model¶

In [10]:

    # Graphical model
    # P(X1,...,X10)
    #        = P(X1,X2)*P(X3,...,X5|X1,X2)*
    #          P(X7,...,X10)*P(X6|X3,...X5,X7,...X10)
    #
    #        = P(X1,...,X5)*P(X3,...,X10)/P(X3,X4,X5)
    #

    estM1_12345 = sum(X1[:,0:5],axis=0)/datanum
    estM2_12345 = sum(X2[:,0:5],axis=0)/datanum
    estM1_34567890 = sum(X1[:,2:10],axis=0)/datanum
    estM2_34567890 = sum(X2[:,2:10],axis=0)/datanum
    estM1_345 = sum(X1[:,2:5],axis=0)/datanum
    estM2_345 = sum(X2[:,2:5],axis=0)/datanum

    estCov1_12345 = np.dot(X1[:,0:5].T, X1[:,0:5])/datanum - np.outer(estM1_12345, estM1_12345)
    estCov2_12345 = np.dot(X2[:,0:5].T, X2[:,0:5])/datanum - np.outer(estM2_12345, estM2_12345)
    estCov1_34567890 = np.dot(X1[:,2:10].T, X1[:,2:10])/datanum - np.outer(estM1_34567890, estM1_34567890)
    estCov2_34567890 = np.dot(X2[:,2:10].T, X2[:,2:10])/datanum - np.outer(estM2_34567890, estM2_34567890)
    estCov1_345 = np.dot(X1[:,2:5].T, X1[:,2:5])/datanum - np.outer(estM1_345, estM1_345)
    estCov2_345 = np.dot(X2[:,2:5].T, X2[:,2:5])/datanum - np.outer(estM2_345, estM2_345)

    probC1tst1 = (exp(getLogGaussian(tstX1[:,0:5], estM1_12345, estCov1_12345))*
        exp(getLogGaussian(tstX1[:,2:10], estM1_34567890, estCov1_34567890))/
        exp(getLogGaussian(tstX1[:,2:5], estM1_345, estCov1_345)))
    probC2tst1 = (exp(getLogGaussian(tstX1[:,0:5], estM2_12345, estCov2_12345))*
        exp(getLogGaussian(tstX1[:,2:10], estM2_34567890, estCov2_34567890))/
        exp(getLogGaussian(tstX1[:,2:5], estM2_345, estCov2_345)))
    probC1tst2 = (exp(getLogGaussian(tstX2[:,0:5], estM1_12345, estCov1_12345))*
        exp(getLogGaussian(tstX2[:,2:10], estM1_34567890, estCov1_34567890))/
        exp(getLogGaussian(tstX2[:,2:5], estM1_345, estCov1_345)))
    probC2tst2 = (exp(getLogGaussian(tstX2[:,0:5], estM2_12345, estCov2_12345))*
        exp(getLogGaussian(tstX2[:,2:10], estM2_34567890, estCov2_34567890))/
        exp(getLogGaussian(tstX2[:,2:5], estM2_345, estCov2_345)))
    GraphicalModelAccuracy = (sum(probC1tst1 - probC2tst1 > 0) +
        sum(probC2tst2 - probC1tst2 > 0))/(2.*tstdatanum)
    print GraphicalModelAccuracy

    0.8955

In [11]:

    print accuracy
    print GraphicalModelAccuracy

    0.796
    0.8955