def Estep(Y, A, mus, sigmas, decay_coef):
	assert(len(Y[0]) == len(A))
	N, D = Y.shape
	D, K = A.shape
	assert((sigmas.shape[0] == D) and (sigmas.shape[1] == 1))
	assert((mus.shape[0] == D) and (mus.shape[1] == 1))
	EX = np.zeros([N, D])
	EXZ = np.zeros([N, D, K])#this is a 3D tensor. 
	EX2 = np.zeros([N, D])
	EZ = np.zeros([N, K])
	EZZT = np.zeros([N, K, K])
	for i in range(N):
		#compute P(Z, X_0 | Y) following three step formula. 
		#1. compute P(Z, X_0)
		Y_i = Y[i, :]
		Y_is_zero = np.abs(Y_i) < 1e-6
		zero_indices = np.array([True for a in range(K)] + [np.abs(Y_i[j]) < 1e-6 for j in range(D)])
		dim = K + D#this is dimension of matrix
		
		#2. compute P(Z, X_0 | Y_+)
		mu_c, sigma_c, augmentedA_0, augmentedA_plus, augmented_D, sigma_22_inv = calcConditionalDistribution(A, mus, sigmas, np.array([np.abs(Y_i[j]) < 1e-6 for j in range(D)]), Y_i[~Y_is_zero])
		
		#3. compute P(Z, X_0 | Y_+, Y_0)
		dim = len(sigma_c)
		matrixToInvert = computeMatrixInLastStep(A, np.abs(Y[i, :]) < 1e-6, sigmas, K, sigma_c, decay_coef, sigma_22_inv)
		if (Y_is_zero).sum() < D:
			magical_matrix = 2 * decay_coef * (mult_diag(augmented_D, matrixToInvert) + augmentedA_0 * (np.eye(K)-augmentedA_plus.T * sigma_22_inv * augmentedA_plus) * (augmentedA_0.T * matrixToInvert))
		else:
			magical_matrix = 2 * decay_coef * (mult_diag(augmented_D, matrixToInvert) + augmentedA_0 * (augmentedA_0.T * matrixToInvert))
		magical_matrix[:, :K] = 0
		sigma_11 = augmentedA_0 * augmentedA_0.T + np.diag(augmented_D)
		if (Y_is_zero).sum() < D:
			sigma_xz = np.array(sigma_c - mult_diag(augmented_D, np.array(magical_matrix), left = False) -  (magical_matrix *augmentedA_0) * ((np.eye(K) - augmentedA_plus.T * sigma_22_inv * augmentedA_plus) * augmentedA_0.T))
		else:
			sigma_xz = np.array(sigma_c - mult_diag(augmented_D, np.array(magical_matrix), left = False) -  (magical_matrix *augmentedA_0) * augmentedA_0.T)
		mu_xz = np.array(np.matrix(np.eye(dim) - magical_matrix) * np.matrix(mu_c))
		EZ[i, :] = mu_xz[:K, 0]
		EX[i, Y_is_zero] = mu_xz[K:, 0]
		EX2[i, Y_is_zero] = mu_xz[K:, 0] ** 2 + np.diag(sigma_xz[K:, K:])
		EZZT[i, :, :] = np.dot(np.atleast_2d(mu_xz[:K, :]), np.atleast_2d(mu_xz[:K,:].transpose()))  + sigma_xz[:K, :K]
		EXZ[i, Y_is_zero, :] =  np.dot(mu_xz[K:], mu_xz[:K].transpose()) + sigma_xz[K:, :K]	
	return EZ, EZZT, EX, EXZ, EX2
estimates the requisite latent expectations in the E-step. 
dataset/ETHPy150Open epierson9/ZIFA/ZIFA/ZIFA.py/Estep