code: def cv_gradient(self,z):
          
          gradient = np.zeros(np.sum(self.approx_param_no))
          z_t = z.T      
          log_q = self.normal_log_q(z.T)
          log_p = self.log_p(z.T)
          grad_log_q = self.grad_log_q(z)
          gradient = grad_log_q*(log_p-log_q)
  
          alpha0 = alpha_recursion(np.zeros(np.sum(self.approx_param_no)), grad_log_q, gradient, np.sum(self.approx_param_no))           
  
          vectorized = gradient - ((alpha0/np.var(grad_log_q,axis=1))*grad_log_q.T).T
  
          return np.mean(vectorized,axis=1)
docs: The control variate augmented Monte Carlo gradient estimate