@@ -0,0 +1,1540 @@
+import numpy as np%0Aimport scipy as sp%0Afrom scipy import special%0Afrom __future__ import division%0Aimport numpy.matlib as M%0A%0Adef gnmf_vb_poisson_mult_fast(x,%0A                            a_tm,%0A                            b_tm,%0A                            a_ve,%0A                            b_ve,%0A                            EPOCH =1000,%0A                            METHOD = 'vb',%0A                            UPDATE = np.inf,%0A                            t_init = np.random.gamma(a_tm, b_tm/a_tm),%0A                            v_init = np.random.gamma(a_ve, b_ve/a_ve),%0A                            tie_a_ve = 'clamp',%0A                            tie_b_ve = 'clamp',%0A                            tie_a_tm = 'clamp',%0A                            tie_b_tm = 'clamp',%0A                            print_period = 500%0A                            ):%0A    W = x.shape%5B0%5D%0A    K = x.shape%5B1%5D%0A    I = b_tm.shape%5B1%5D%0A%0A    M = ~np.isnan(x)%0A    X = np.zeros(x.shape)%0A    X%5BM%5D = x%5BM%5D%0A%0A    L_t = t_init%0A    L_v = v_init%0A    E_t = t_init%0A    Sig_t = t_init%0A    Sig_v = v_init%0A%0A    B = np.zeros(1,EPOCH)%0A    gammalnX = special.gammaln(X+1)%0A%0A    for e in range(EPOCH):%0A        %0A        LtLv = L_t.dot(L_v)%0A        tmp = X / (LtLv)%0A        #check Tranpose%0A        Sig_t = L_t.dot(tmp.dot(L_v.T))%0A        Sig_v = L_v.dot(L_t.T.dot(tmp))%0A%0A        alpha_tm = a_tm + Sig_t%0A        beta_tm = 1/((a_tm/b_tm) + M.dot(E_v.T))%0A        E_t = alpha_tm.dot(beta_tm)%0A%0A        alpha_ve = a_ve + Sig_v%0A        beta_ve = 1/((a_ve/b_ve) + E_t.T.dot(M))%0A%0A        E_v = alpha_ve.dot(beta_ve)%0A

gnmf_vb_poisson_multi_fast mark 1