diff --git a/mvbernoulli/src/ising_model.cpp b/mvbernoulli/src/ising_model.cpp
index 38f83e5..329d22d 100644
--- a/mvbernoulli/src/ising_model.cpp
+++ b/mvbernoulli/src/ising_model.cpp
@@ -310,7 +310,7 @@ Rcpp::NumericVector ising_marginal_probs(const VechView& theta) {
     return score;
 }
 
-//' Natural parameters from the sufficient conditional probability statistis `pi`
+//' Natural parameters from the sufficient conditional probability statistic `pi`
 //'
 //' Computes the natural parameters `theta` of the Ising model from zero
 //' conditioned probabilites for single and two way effects.
@@ -348,7 +348,7 @@ Rcpp::NumericVector ising_theta_from_cond_prob(const VechView& pi) {
     return theta;
 }
 
-//' Computes the log-lokelihood at natural parameters `theta` of the Ising model
+//' Computes the log-likelihood at natural parameters `theta` of the Ising model
 //' given a set of observations `Y`
 //'
 //'      l(theta) = log(p_0(theta)) + n^-1 sum_i vech(y_i y_i')' theta
@@ -583,6 +583,8 @@ Rcpp::NumericVector ising_conditional_score(const Rcpp::NumericMatrix& alpha,
 //       as the garbage collector stack gets out of sync.
 // SEE: https://stackoverflow.com/questions/42119609/how-to-handle-mismatching-in-calling-convention-in-r-and-c-using-rcpp
 //
+// NOTE: The below algorithm implements the _old_ `theta` to `eta` relation
+//
 // // [[Rcpp::export(rng = false)]]
 // Rcpp::NumericVector ising_conditional_score_mt(const Rcpp::NumericMatrix& alpha,
 //     const Rcpp::NumericMatrix& X, const MVBinary& Y
diff --git a/tensorPredictors/R/make_gmlm_family.R b/tensorPredictors/R/make_gmlm_family.R
index 8d06a93..4027680 100644
--- a/tensorPredictors/R/make_gmlm_family.R
+++ b/tensorPredictors/R/make_gmlm_family.R
@@ -132,17 +132,58 @@ make.gmlm.family <- function(name) {
                 q <- head(dim(Fy), -1)      # response dimensions
                 r <- length(p)              # single predictor/response tensor order
 
-                # Hij = Cov(ti(X) %x% tj(X) | Y = y), i, j = 1, 2
-                H11 <- Reduce(`%x%`, rev(Map(solve, Omegas)))   # covariance
-                # 3rd central moment is zero
-                H12 <- H21 <- 0
-                # 4th moment by "Isserlis' theorem" a.k.a. "Wick's theorem"
-                H22 <- kronecker(H11, H11)
-                dim(H22) <- rep(prod(p), 4)
-                H22 <- H22 + aperm(H22, c(1, 3, 2, 4)) + aperm(H22, c(1, 3, 2, 4))
+                # Note: independent of Y
+                Deltas <- Map(solve, Omegas)
 
-                dim(H11) <- c(p, p)
-                dim(H22) <- c(p, p, p, p)
+
+                # Conditional Mean E1 = mu_y = E[X | Y = y]
+
+
+
+                # Get the i#th observation from Fy
+                FyGet <- function(i) {
+                    obs <- eval(str2lang(paste(
+                        "Fy[",
+                        paste(rep(",", r), collapse = ""),
+                        "i, drop = FALSE]"
+                    , collapse = "")))
+                    dim(obs) <- head(dim(obs), -1)
+                    obs
+                }
+
+                E1 <- mlm(mlm(Fy, alphas) + c(eta1), Deltas)
+
+                ## Hij = Cov(t_i(X), t_j(X) | Y = y) for i, j = 1, 2
+                #      = partial^2 b(eta_y) / partial eta_yi partial eta_yj
+                #
+                # Cov(vec(X) | Y = y)
+                # Note: Independent of Y
+                H11 <- Reduce(`%x%`, rev(Deltas))
+                # Cov(vec(X), vec(X) %x% vec(X) | Y = y)
+                #   = (mu_y %o% Sigma)_(3, (1, 2)) + (mu_y %o% Sigma)_(3, (2, 1))
+                #   = (Sigma %o% mu_y)_(1, (2, 3)) + (Sigma %o% mu_y)_(1, (3, 2))
+                H12 <- outer(H11, E1)
+                H12 <- H12 + aperm(H12, c(1, 3, 2))
+                H21 <- aperm(H12, c(2, 3, 1))
+
+                # H22 = Cov(vec X %x% vec X | Y = y)
+                H22 <- local({
+                    # mu_i %o% mu_i %o% Sigma + Sigma %o% Sigma + Sigma %o% mu_i %o% mu_i
+                    h22 <- outer(outer(E1, E1), H11) + outer(H11, H11) + outer(H11, outer(E1, E1))
+                    aperm(h22, c(1, 3, 2, 4)) + aperm(h22, c(1, 3, 4, 2))
+                })
+
+                # # Hij = Cov(ti(X) %x% tj(X) | Y = y), i, j = 1, 2
+                # H11 <- Reduce(`%x%`, rev(Map(solve, Omegas)))   # covariance
+                # # 3rd central moment is zero
+                # H12 <- H21 <- 0
+                # # 4th moment by "Isserlis' theorem" a.k.a. "Wick's theorem"
+                # H22 <- kronecker(H11, H11)
+                # dim(H22) <- rep(prod(p), 4)
+                # H22 <- H22 + aperm(H22, c(1, 3, 2, 4)) + aperm(H22, c(1, 3, 2, 4))
+
+                # dim(H11) <- c(p, p)
+                # dim(H22) <- c(p, p, p, p)
 
                 ## Fisher Information: Tensor Normal Specific
                 # known exponential family constants