

In [1]:

    
from numpy import pi, log
from numpy.linalg import inv, det
import numpy as np



In [2]:

    
v1 = np.array([ 8, 11, 16, 18,  6,  4, 20, 25, 9, 13])
v2 = np.array([10, 14, 16, 15, 20,  4, 18, 22])

$$ \hat{\mu}_1 = \sum_{j=1}^n w_{1,j}/n \\ \hat{\sigma}_{1,1} = \sum_{j=1}^n (w_{1,j} - \hat{\mu}_1)^2/n \\ \hat{\mu}_2 = \bar{w}_2 + \hat{\beta}(\hat{\mu}_1 - \bar{w}_1) \\ \hat{\sigma}_{2,2} = s_{2,2} + \hat{\beta}^2(\hat{\sigma}_{1,1} - s_{1,1}) \\ \hat{\sigma}_{1,2} = \hat{\beta}\hat{\sigma}_{1,1} \\ \hat{\beta} = s_{1,2}/s_{1,1} $$

where

$$ \bar{w}_i = \sum_{j=1}^m w_{i,j}/m $$

and

$$ s_{h,i} = \sum_{j=1}^m (w_{h,j} - \bar{y}_h)(w_{i,j} - \bar{y}_i)/m $$



In [3]:

    
v1.mean()









    Out[3]:





13.0



In [4]:

    
sigma11 = v1.std()**2
sigma11









    Out[4]:





40.20000000000001



In [5]:

    
v2.mean()









    Out[5]:





14.875



In [6]:

    
mu = np.array([[v1.mean(), v2.mean()]])
print mu.shape
print mu









    



(1L, 2L)
[[ 13.     14.875]]



In [9]:

    
s11 = 0
for _ in range(len(v1)):
    s11 += (v1[_] - v1.mean())*(v1[_] - v1.mean()) / len(v1)
s11









    Out[9]:





40.199999999999996



In [10]:

    
s12 = 0
for _ in range(len(v2)):
    s12 += (v1[_] - v1.mean())*(v2[_] - v2.mean()) / len(v2)
s12









    Out[10]:





24.9375



In [11]:

    
s22 = 0
for _ in range(len(v2)):
    s22 += (v2[_] - v2.mean())*(v2[_] - v2.mean()) / len(v2)
s22









    Out[11]:





28.859375



In [12]:

    
beta = s12/s11
beta









    Out[12]:





0.62033582089552242



In [13]:

    
sigma12 = beta * v1.std()**2
sigma12









    Out[13]:





24.937500000000007



In [14]:

    
sigma22 = s22 + beta*(sigma11 - s11)
sigma22









    Out[14]:





28.859375000000007



In [15]:

    
n = len(v1)



In [16]:

    
m1 = 0
m2 = 2
m = n - m1 - m2
print 'm: %d; m1: %d; m2: %d; n: %d' % (m, m1, m2, n)









    



m: 8; m1: 0; m2: 2; n: 10



In [17]:

    
covM = np.array([[sigma11, sigma12], [sigma12, sigma22]])
covM









    Out[17]:





array([[ 40.2     ,  24.9375  ],
       [ 24.9375  ,  28.859375]])



In [18]:

    
half = 1./2.



In [19]:

    
V = np.concatenate([v1[:m].reshape(-1,1), v2.reshape(-1,1)], axis=1)
V









    Out[19]:





array([[ 8, 10],
       [11, 14],
       [16, 16],
       [18, 15],
       [ 6, 20],
       [ 4,  4],
       [20, 18],
       [25, 22]])



In [20]:

    
ss = 0
for _ in range(m):
    ss += V[_].dot(inv(covM)).dot(V[_])
print ss









    



80.8610604326



In [21]:

    
ss2 = 0
for _ in range(m, n):
    ss2 += (v1[_] - mu[0,0])**2
ss2









    Out[21]:





16.0



In [22]:

    
l = -(m + half*(m1+m2)*log(2. * pi)) - half*m*log(det(covM)) \
-half*ss - half*(m1*log(sigma22) + m2*log(sigma11)) \
-half*(1./sigma11 * ss2)



In [23]:

    
print l * -2









    



158.629410662



In [35]:

    
def ss(x, y):
    return ((x-x.mean())*(y-y.mean())).sum()



In [34]:

    
w1 = v1[:m].mean()
w2 = v2.mean()
print 'w1 %.6f - w2 %.6f' % (w1, w2)









    



w1 13.500000 - w2 14.875000



In [40]:

    
ss(v1[:m], v1[:m])/m









    Out[40]:





48.0



In [41]:

    
ss(v2,v2)/m









    Out[41]:





28.859375



In [42]:

    
ss(v1[:m], v2)/m









    Out[42]:





24.9375



In [ ]: