-
Notifications
You must be signed in to change notification settings - Fork 0
/
MEGA_Numpy.py
65 lines (47 loc) · 1.41 KB
/
MEGA_Numpy.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
# -*- coding: utf-8 -*-
"""
Created on 12/02/21
@author: Cedric Beaulac
LVM-MEGA : Implementation for MEGA1 and MEGA2
"""
import numpy as np
##################################
# frobnorm : Frobenius norm of the matrix A
##################################
def frobnorm(A):
A2 = np.dot(A,A)
Trace = np.trace(A2)
Frob = np.sqrt(Trace)
return Frob
##################################
# frobvect : Frobenius norm of the vector a
##################################
def frobvect(a):
Frob = np.sqrt(np.sum(np.power(a,2)))
return Frob
#################################
# MEGA: Returns MEGA1 and MEGA2 (using Frob norm)
# Inputs: data (Data set), mean (sample of E[x|z]), var (sample of V[x|z])
# Output : MEGA1 (MEGA for first moment) and MEGA2 (MEGA for second moment)
#################################
def MEGA(data,mean,var):
#Initialize MEGA
MEGA = np.zeros(2)
nz = np.shape(mean)[0]
d = np.shape(mean)[1]
#First moment MEGA
#DE
xbar = np.mean(data,0)
#FME
EzEx = np.mean(mean,0)
MEGA[0] = frobvect(EzEx-xbar)
# Second moment MEGA
# LHS (DE)
S2 = np.cov(np.transpose(data))
xbar2 = np.outer(xbar,np.transpose(xbar))
LHS = S2+xbar2
# RHS (FME)
E2 = np.matmul(mean.reshape(-1,np.shape(mean)[1],1),mean.reshape(-1,1,np.shape(mean)[1]))
RHS = np.mean(E2+var,0)
MEGA[1] = frobnorm(LHS-RHS)
return MEGA