code: def _twofilter_smoothing_ON2(self, t, ti, info, phi, lwinfo):
          
          sp, sw = 0., 0.
          upb = lwinfo.max() + self.wgt[t].lw.max()
          if hasattr(self.model, 'upper_bound_trans'):
              upb += self.model.upper_bound_trans(t + 1)
          # Loop over n, to avoid having in memory a NxN matrix
          for n in range(self.N):
              omegan = np.exp(lwinfo + self.wgt[t].lw[n] - upb
                              + self.model.logpt(t + 1, self.X[t][n],
                                                 info.hist.X[ti]))
              sp += np.sum(omegan * phi(self.X[t][n], info.hist.X[ti]))
              sw += np.sum(omegan)
          return sp / sw
docs: O(N^2) version of two-filter smoothing. 
  
          This method should not be called directly, see twofilter_smoothing.