def METHOD_NAME(lam, k, voxelsize):
    lpk = lp_poisson_n(lam, k)
    pk = np.exp(lpk)
    #print 'number of NaNs = %d' % isnan(pk).sum()
    dx, dy, dz = np.gradient(lpk)
    dx *= 1./voxelsize[0]
    dy *= 1./voxelsize[1]
    dz *= 1./voxelsize[2]
    Exx = (dx*dx*pk).sum(1).sum(0)
    Eyy = (dy*dy*pk).sum(1).sum(0)
    Ezz = (dz*dz*pk).sum(1).sum(0)
    Exy = (dx*dy*pk).sum(1).sum(0)
    Exz = (dx*dz*pk).sum(1).sum(0)
    Eyz = (dy*dz*pk).sum(1).sum(0)
    return np.array([([[Exx[i], Exy[i], Exz[i]],[Exy[i],Eyy[i], Eyz[i]],[Exz[i],Eyz[i], Ezz[i]]]) for i in range(lam.shape[2])])