@staticmethod<EOL><INDENT>def _EM_gradient(A,Q,H,R,m_init,P_init,measurements, M, P, G, dA, dQ, dH, dR, dm_init, dP_init):<DEDENT>
import pdb; pdb.set_trace();<EOL>param_number = dA.shape[-<NUM_LIT:1>]<EOL>d_log_likelihood_update = np.empty((param_number,<NUM_LIT:1>))<EOL>sample_no = measurements.shape[<NUM_LIT:0>]<EOL>P_1 = P[<NUM_LIT:1>:,:,:] <EOL>P_2 = P[<NUM_LIT:0>:-<NUM_LIT:1>,:,:] <EOL>M_1 = M[<NUM_LIT:1>:,:] <EOL>M_2 = M[<NUM_LIT:0>:-<NUM_LIT:1>,:] <EOL>Sigma = np.mean(P_1,axis=<NUM_LIT:0>) + np.dot(M_1.T, M_1) / sample_no <EOL>Phi =   np.mean(P_2,axis=<NUM_LIT:0>) + np.dot(M_2.T, M_2) / sample_no <EOL>B = np.dot( measurements.T, M_1 )/ sample_no<EOL>C =   (sp.einsum( '<STR_LIT>', P_1, G) + np.dot(M_1.T, M_2)) / sample_no <EOL><INDENT>C1 = np.zeros( (P_1.shape[<NUM_LIT:1>],P_1.shape[<NUM_LIT:1>]) )<EOL>for k in range(P_1.shape[<NUM_LIT:0>]):<EOL><INDENT>C1 += np.dot(P_1[k,:,:],G[k,:,:]) + sp.outer( M_1[k,:], M_2[k,:] )<EOL><DEDENT>C1 = C1 / sample_no<EOL><DEDENT>D = np.dot( measurements.T, measurements ) / sample_no<EOL>try:<EOL><INDENT>P_init_inv = sp.linalg.inv(P_init)<EOL>if np.max( np.abs(P_init_inv)) > <NUM_LIT>:<EOL><INDENT>compute_P_init_terms = False<EOL><DEDENT>else:<EOL><INDENT>compute_P_init_terms = True<EOL><DEDENT><DEDENT>except np.linalg.LinAlgError:<EOL><INDENT>compute_P_init_terms = False<EOL><DEDENT>try:<EOL><INDENT>Q_inv = sp.linalg.inv(Q)<EOL>if np.max( np.abs(Q_inv)) > <NUM_LIT>:<EOL><INDENT>compute_Q_terms = False<EOL><DEDENT>else:<EOL><INDENT>compute_Q_terms = True<EOL><DEDENT><DEDENT>except np.linalg.LinAlgError:<EOL><INDENT>compute_Q_terms = False<EOL><DEDENT>try:<EOL><INDENT>R_inv = sp.linalg.inv(R)<EOL>if np.max( np.abs(R_inv)) > <NUM_LIT>:<EOL><INDENT>compute_R_terms = False<EOL><DEDENT>else:<EOL><INDENT>compute_R_terms = True<EOL><DEDENT><DEDENT>except np.linalg.LinAlgError:<EOL><INDENT>compute_R_terms = False<EOL><DEDENT>d_log_likelihood_update = np.zeros((param_number,<NUM_LIT:1>))<EOL>for j in range(param_number):<EOL><INDENT>if compute_P_init_terms:<EOL><INDENT>d_log_likelihood_update[j,:] -= <NUM_LIT:0.5> * np.sum(P_init_inv* dP_init[:,:,j].T ) <EOL>M0_smoothed = M[<NUM_LIT:0>]; M0_smoothed.shape = (M0_smoothed.shape[<NUM_LIT:0>],<NUM_LIT:1>)<EOL>tmp1 = np.dot( dP_init[:,:,j], np.dot( P_init_inv, (P[<NUM_LIT:0>,:,:] + sp.outer( (M0_smoothed - m_init), (M0_smoothed - m_init) )) )  ) <EOL>d_log_likelihood_update[j,:] += <NUM_LIT:0.5> * np.sum(P_init_inv* tmp1.T )<EOL>tmp2 = sp.outer( dm_init[:,j], M0_smoothed )<EOL>tmp2 += tmp2.T<EOL>d_log_likelihood_update[j,:] += <NUM_LIT:0.5> * np.sum(P_init_inv* tmp2.T )<EOL><DEDENT>if compute_Q_terms:<EOL><INDENT>d_log_likelihood_update[j,:] -=  sample_no/<NUM_LIT> * np.sum(Q_inv* dQ[:,:,j].T ) <EOL>tmp1 = np.dot(C,A.T); tmp1 += tmp1.T; tmp1 = Sigma - tmp1 + np.dot(A, np.dot(Phi,A.T)) <EOL>tmp1 = np.dot( dQ[:,:,j], np.dot( Q_inv, tmp1) )<EOL>d_log_likelihood_update[j,:] += sample_no/<NUM_LIT> * np.sum(Q_inv * tmp1.T)<EOL>tmp2 = np.dot( dA[:,:,j], C.T); tmp2 += tmp2.T;<EOL>tmp3 = np.dot(dA[:,:,j], np.dot(Phi,A.T)); tmp3 += tmp3.T<EOL>d_log_likelihood_update[j,:] -= sample_no/<NUM_LIT> * np.sum(Q_inv.T * (tmp3 - tmp2) )<EOL><DEDENT>if compute_R_terms:<EOL><INDENT>d_log_likelihood_update[j,:] -=  sample_no/<NUM_LIT> * np.sum(R_inv* dR[:,:,j].T )<EOL>tmp1 = np.dot(B,H.T); tmp1 += tmp1.T; tmp1 = D - tmp1 + np.dot(H, np.dot(Sigma,H.T))<EOL>tmp1 = np.dot( dR[:,:,j], np.dot( R_inv, tmp1) )<EOL>d_log_likelihood_update[j,:] += sample_no/<NUM_LIT> * np.sum(R_inv * tmp1.T)<EOL>tmp2 = np.dot( dH[:,:,j], B.T); tmp2 += tmp2.T;<EOL>tmp3 = np.dot(dH[:,:,j], np.dot(Sigma,H.T)); tmp3 += tmp3.T<EOL>d_log_likelihood_update[j,:] -= sample_no/<NUM_LIT> * np.sum(R_inv.T * (tmp3 - tmp2) )<EOL><DEDENT><DEDENT>return d_log_likelihood_update<EOL>
Gradient computation with the EM algorithm.

Input:
-----------------

M: Means from the smoother
P: Variances from the smoother
G: Gains? from the smoother