diff --git a/src/fastcpd_class_nll.cc b/src/fastcpd_class_nll.cc
index b7444a3e..f3b6b96c 100644
--- a/src/fastcpd_class_nll.cc
+++ b/src/fastcpd_class_nll.cc
@@ -2,6 +2,327 @@
 
 namespace fastcpd::classes {
 
+colvec Fastcpd::get_gradient_arma(
+  const unsigned int segment_start,
+  const unsigned int segment_end,
+  const colvec& theta
+) {
+  const mat data_segment = data.rows(segment_start, segment_end);
+  const unsigned int segment_length = segment_end - segment_start + 1;
+  mat reversed_data = reverse(data_segment, 0);
+  colvec reversed_theta = reverse(theta);
+  if (segment_length < max(order) + 1) {
+    return ones(theta.n_elem);
+  }
+  colvec variance_term = zeros(segment_length);
+  for (unsigned int i = max(order); i < segment_length; i++) {
+    variance_term(i) = data_segment(i, 0) - dot(
+        reversed_theta.rows(order(1) + 1, sum(order)),
+        data_segment.rows(i - order(0), i - 1)
+      ) - dot(
+        reversed_theta.rows(1, order(1)),
+        variance_term.rows(i - order(1), i - 1)
+      );
+  }
+  colvec reversed_variance_term = reverse(variance_term);
+  mat phi_coefficient = zeros(segment_length, order(0)),
+      psi_coefficient = zeros(segment_length, order(1));
+  for (unsigned int i = max(order); i < segment_length; i++) {
+    phi_coefficient.row(i) = -reversed_data.rows(
+      segment_length - i, segment_length - i + order(0) - 1
+    ).t() - reversed_theta.rows(1, order(1)).t() *
+    phi_coefficient.rows(i - order(1), i - 1);
+  }
+  for (unsigned int i = order(1); i < segment_length; i++) {
+    psi_coefficient.row(i) = -reversed_variance_term.rows(
+        segment_length - i, segment_length - i + order(1) - 1
+      ).t() - reversed_theta.rows(1, order(1)).t() *
+      psi_coefficient.rows(i - order(1), i - 1);
+  }
+  colvec gradient = zeros(sum(order) + 1);
+  gradient.rows(0, order(0) - 1) =
+    phi_coefficient.row(segment_length - 1).t() *
+    variance_term(segment_length - 1) / theta(sum(order));
+  gradient.rows(order(0), sum(order) - 1) =
+    psi_coefficient.row(segment_length - 1).t() *
+    variance_term(segment_length - 1) / theta(sum(order));
+  gradient(sum(order)) = 1.0 / 2.0 / theta(sum(order)) -
+    std::pow(variance_term(segment_length - 1), 2) / 2.0 /
+    std::pow(theta(sum(order)), 2);
+  return gradient;
+}
+
+colvec Fastcpd::get_gradient_binomial(
+  const unsigned int segment_start,
+  const unsigned int segment_end,
+  const colvec& theta
+) {
+  const mat data_segment = data.rows(segment_start, segment_end);
+  const unsigned int segment_length = segment_end - segment_start + 1;
+  rowvec new_data = data_segment.row(segment_length - 1);
+  rowvec x = new_data.tail(new_data.n_elem - 1);
+  double y = new_data(0);
+  return - (y - 1 / (1 + exp(-as_scalar(x * theta)))) * x.t();
+}
+
+colvec Fastcpd::get_gradient_lm(
+  const unsigned int segment_start,
+  const unsigned int segment_end,
+  const colvec& theta
+) {
+  const mat data_segment = data.rows(segment_start, segment_end);
+  const unsigned int segment_length = segment_end - segment_start + 1;
+  rowvec new_data = data_segment.row(segment_length - 1);
+  rowvec x = new_data.tail(new_data.n_elem - 1);
+  double y = new_data(0);
+  return - (y - as_scalar(x * theta)) * x.t();
+}
+
+colvec Fastcpd::get_gradient_ma(
+  const unsigned int segment_start,
+  const unsigned int segment_end,
+  const colvec& theta
+) {
+  const mat data_segment = data.rows(segment_start, segment_end);
+  const unsigned int segment_length = segment_end - segment_start + 1;
+  const unsigned int q = order(1);
+  mat reversed_data = reverse(data_segment, 0);
+  colvec reversed_theta = reverse(theta);
+  if (segment_length < q + 1) {
+    return ones(theta.n_elem);
+  }
+  colvec variance_term = zeros(segment_length);
+  for (unsigned int i = q; i < segment_length; i++) {
+    variance_term(i) = data_segment(i, 0) -
+      dot(reversed_theta.rows(1, q), variance_term.rows(i - q, i - 1));
+  }
+  colvec reversed_variance_term = reverse(variance_term);
+  mat psi_coefficient = zeros(segment_length, q);
+  for (unsigned int i = q; i < segment_length; i++) {
+    psi_coefficient.row(i) = -reversed_variance_term.rows(
+        segment_length - i, segment_length - i + q - 1
+      ).t() - reversed_theta.rows(1, q).t() *
+      psi_coefficient.rows(i - q, i - 1);
+  }
+  colvec gradient = zeros(q + 1);
+  gradient.rows(0, q - 1) =
+    psi_coefficient.row(segment_length - 1).t() *
+    variance_term(segment_length - 1) / theta(q);
+  gradient(q) = 1.0 / 2.0 / theta(q) -
+    std::pow(variance_term(segment_length - 1), 2) / 2.0 /
+    std::pow(theta(q), 2);
+  return gradient;
+}
+
+colvec Fastcpd::get_gradient_poisson(
+  const unsigned int segment_start,
+  const unsigned int segment_end,
+  const colvec& theta
+) {
+  const mat data_segment = data.rows(segment_start, segment_end);
+  const unsigned int segment_length = segment_end - segment_start + 1;
+  rowvec new_data = data_segment.row(segment_length - 1);
+  rowvec x = new_data.tail(new_data.n_elem - 1);
+  double y = new_data(0);
+  return - (y - exp(as_scalar(x * theta))) * x.t();
+}
+
+mat Fastcpd::get_hessian_arma(
+  const unsigned int segment_start,
+  const unsigned int segment_end,
+  const colvec& theta
+) {
+  const mat data_segment = data.rows(segment_start, segment_end);
+  const unsigned int segment_length = segment_end - segment_start + 1;
+  // TODO(doccstat): Maybe we can store all these computations
+  mat reversed_data = reverse(data_segment, 0);
+  colvec reversed_theta = reverse(theta);
+  if (segment_length < max(order) + 1) {
+    return eye(theta.n_elem, theta.n_elem);
+  }
+  colvec variance_term = zeros(segment_length);
+  for (unsigned int i = max(order); i < segment_length; i++) {
+    variance_term(i) = data_segment(i, 0) - dot(
+        reversed_theta.rows(order(1) + 1, sum(order)),
+        data_segment.rows(i - order(0), i - 1)
+      ) - dot(
+        reversed_theta.rows(1, order(1)),
+        variance_term.rows(i - order(1), i - 1)
+      );
+  }
+  colvec reversed_variance_term = reverse(variance_term);
+  mat phi_coefficient = zeros(segment_length, order(0)),
+      psi_coefficient = zeros(segment_length, order(1));
+  for (unsigned int i = max(order); i < segment_length; i++) {
+    phi_coefficient.row(i) = -reversed_data.rows(
+      segment_length - i, segment_length - i + order(0) - 1
+    ).t() - reversed_theta.rows(1, order(1)).t() *
+    phi_coefficient.rows(i - order(1), i - 1);
+  }
+  for (unsigned int i = order(1); i < segment_length; i++) {
+    psi_coefficient.row(i) = -reversed_variance_term.rows(
+        segment_length - i, segment_length - i + order(1) - 1
+      ).t() - reversed_theta.rows(1, order(1)).t() *
+      psi_coefficient.rows(i - order(1), i - 1);
+  }
+  mat reversed_coef_phi = reverse(phi_coefficient, 0),
+      reversed_coef_psi = reverse(psi_coefficient, 0);
+  cube phi_psi_coefficient = zeros(order(1), order(0), segment_length),
+        psi_psi_coefficient = zeros(order(1), order(1), segment_length);
+  for (unsigned int i = order(1); i < segment_length; i++) {
+    mat phi_psi_coefficient_part = zeros(order(1), order(0)),
+        psi_psi_coefficient_part = zeros(order(1), order(1));
+    for (unsigned int j = 1; j <= order(1); j++) {
+      phi_psi_coefficient_part +=
+        phi_psi_coefficient.slice(i - j) * theta(order(0) - 1 + j);
+    }
+    phi_psi_coefficient.slice(i) = -reversed_coef_phi.rows(
+      segment_length - i, segment_length - i + order(1) - 1
+    ) - phi_psi_coefficient_part;
+    for (unsigned int j = 1; j <= order(1); j++) {
+      psi_psi_coefficient_part +=
+        psi_psi_coefficient.slice(i - j) * theta(order(0) - 1 + j);
+    }
+    psi_psi_coefficient.slice(i) = -reversed_coef_psi.rows(
+        segment_length - i, segment_length - i + order(1) - 1
+      ) - reversed_coef_psi.rows(
+        segment_length - i, segment_length - i + order(1) - 1
+      ).t() - psi_psi_coefficient_part;
+  }
+  mat hessian = zeros(sum(order) + 1, sum(order) + 1);
+  hessian.submat(0, 0, order(0) - 1, order(0) - 1) =
+    phi_coefficient.row(segment_length - 1).t() *
+    phi_coefficient.row(segment_length - 1) / theta(sum(order));
+  hessian.submat(0, order(0), order(0) - 1, sum(order) - 1) = (
+    phi_psi_coefficient.slice(segment_length - 1).t() *
+      variance_term(segment_length - 1) +
+      phi_coefficient.row(segment_length - 1).t() *
+      psi_coefficient.row(segment_length - 1)
+  ) / theta(sum(order));
+  hessian.submat(order(0), 0, sum(order) - 1, order(0) - 1) =
+    hessian.submat(0, order(0), order(0) - 1, sum(order) - 1).t();
+  hessian.submat(0, sum(order), order(0) - 1, sum(order)) =
+    -phi_coefficient.row(segment_length - 1).t() *
+    variance_term(segment_length - 1) / theta(sum(order)) / theta(sum(order));
+  hessian.submat(sum(order), 0, sum(order), order(0) - 1) =
+    hessian.submat(0, sum(order), order(0) - 1, sum(order)).t();
+  hessian.submat(order(0), order(0), sum(order) - 1, sum(order) - 1) = (
+    psi_coefficient.row(segment_length - 1).t() *
+    psi_coefficient.row(segment_length - 1) +
+    psi_psi_coefficient.slice(segment_length - 1) *
+    variance_term(segment_length - 1)
+  ) / theta(sum(order));
+  hessian.submat(order(0), sum(order), sum(order) - 1, sum(order)) =
+    -psi_coefficient.row(segment_length - 1).t() *
+    variance_term(segment_length - 1) / theta(sum(order)) / theta(sum(order));
+  hessian.submat(sum(order), order(0), sum(order), sum(order) - 1) =
+    hessian.submat(order(0), sum(order), sum(order) - 1, sum(order)).t();
+  hessian(sum(order), sum(order)) =
+    std::pow(variance_term(segment_length - 1), 2) /
+    std::pow(theta(sum(order)), 3) -
+    1.0 / 2.0 / std::pow(theta(sum(order)), 2);
+  return hessian;
+}
+
+mat Fastcpd::get_hessian_binomial(
+  const unsigned int segment_start,
+  const unsigned int segment_end,
+  const colvec& theta
+) {
+  const mat data_segment = data.rows(segment_start, segment_end);
+  const unsigned int segment_length = segment_end - segment_start + 1;
+  rowvec new_data = data_segment.row(segment_length - 1);
+  rowvec x = new_data.tail(new_data.n_elem - 1);
+  double prob = 1 / (1 + exp(-dot(x, theta)));
+  return (x.t() * x) * ((1 - prob) * prob);
+}
+
+mat Fastcpd::get_hessian_lm(
+  const unsigned int segment_start,
+  const unsigned int segment_end,
+  const colvec& theta
+) {
+  const mat data_segment = data.rows(segment_start, segment_end);
+  const unsigned int segment_length = segment_end - segment_start + 1;
+  rowvec new_data = data_segment.row(segment_length - 1);
+  rowvec x = new_data.tail(new_data.n_elem - 1);
+  return x.t() * x;
+}
+
+mat Fastcpd::get_hessian_ma(
+  const unsigned int segment_start,
+  const unsigned int segment_end,
+  const colvec& theta
+) {
+  const mat data_segment = data.rows(segment_start, segment_end);
+  const unsigned int segment_length = segment_end - segment_start + 1;
+  const unsigned int q = order(1);
+  // TODO(doccstat): Maybe we can store all these computations
+  mat reversed_data = reverse(data_segment, 0);
+  colvec reversed_theta = reverse(theta);
+  if (segment_length < q + 1) {
+    return eye(theta.n_elem, theta.n_elem);
+  }
+  colvec variance_term = zeros(segment_length);
+  for (unsigned int i = q; i < segment_length; i++) {
+    variance_term(i) = data_segment(i, 0) - dot(
+      reversed_theta.rows(1, q), variance_term.rows(i - q, i - 1)
+    );
+  }
+  colvec reversed_variance_term = reverse(variance_term);
+  mat psi_coefficient = zeros(segment_length, q);
+  for (unsigned int i = q; i < segment_length; i++) {
+    psi_coefficient.row(i) = -reversed_variance_term.rows(
+        segment_length - i, segment_length - i + q - 1
+      ).t() - reversed_theta.rows(1, q).t() *
+      psi_coefficient.rows(i - q, i - 1);
+  }
+  mat reversed_coef_psi = reverse(psi_coefficient, 0);
+  cube psi_psi_coefficient = zeros(q, q, segment_length);
+  for (unsigned int i = q; i < segment_length; i++) {
+    mat psi_psi_coefficient_part = zeros(q, q);
+    for (unsigned int j = 1; j <= q; j++) {
+      psi_psi_coefficient_part +=
+        psi_psi_coefficient.slice(i - j) * theta(j - 1);
+    }
+    psi_psi_coefficient.slice(i) = -reversed_coef_psi.rows(
+        segment_length - i, segment_length - i + q - 1
+      ) - reversed_coef_psi.rows(
+        segment_length - i, segment_length - i + q - 1
+      ).t() - psi_psi_coefficient_part;
+  }
+  mat hessian = zeros(q + 1, q + 1);
+  hessian.submat(0, 0, q - 1, q - 1) = (
+    psi_coefficient.row(segment_length - 1).t() *
+    psi_coefficient.row(segment_length - 1) +
+    psi_psi_coefficient.slice(segment_length - 1) *
+    variance_term(segment_length - 1)
+  ) / theta(q);
+  hessian.submat(0, q, q - 1, q) =
+    -psi_coefficient.row(segment_length - 1).t() *
+    variance_term(segment_length - 1) / theta(q) / theta(q);
+  hessian.submat(q, 0, q, q - 1) = hessian.submat(0, q, q - 1, q).t();
+  hessian(q, q) =
+    std::pow(variance_term(segment_length - 1), 2) /
+    std::pow(theta(q), 3) -
+    1.0 / 2.0 / std::pow(theta(q), 2);
+  return hessian;
+}
+
+mat Fastcpd::get_hessian_poisson(
+  const unsigned int segment_start,
+  const unsigned int segment_end,
+  const colvec& theta
+) {
+  const mat data_segment = data.rows(segment_start, segment_end);
+  const unsigned int segment_length = segment_end - segment_start + 1;
+  rowvec new_data = data_segment.row(segment_length - 1);
+  rowvec x = new_data.tail(new_data.n_elem - 1);
+  double prob = exp(as_scalar(x * theta));
+  // Prevent numerical issues if `prob` is too large.
+  return (x.t() * x) * std::min(as_scalar(prob), 1e10);
+}
+
 CostResult Fastcpd::get_nll_arma(
   const unsigned int segment_start,
   const unsigned int segment_end
diff --git a/src/fastcpd_class_private.cc b/src/fastcpd_class_private.cc
index bc06b89f..a5e90f4a 100644
--- a/src/fastcpd_class_private.cc
+++ b/src/fastcpd_class_private.cc
@@ -32,7 +32,7 @@ void Fastcpd::create_cost_function_wrapper(Nullable<Function> cost) {
 void Fastcpd::create_cost_gradient_wrapper(Nullable<Function> cost_gradient) {
   if (contain(FASTCPD_FAMILIES, family)) {
     cost_gradient_wrapper = std::bind(  // # nocov start
-      &Fastcpd::cost_update_gradient,  // # nocov end
+      get_gradient,  // # nocov end
       this,
       std::placeholders::_1,
       std::placeholders::_2,
@@ -53,7 +53,7 @@ void Fastcpd::create_cost_gradient_wrapper(Nullable<Function> cost_gradient) {
 void Fastcpd::create_cost_hessian_wrapper(Nullable<Function> cost_hessian) {
   if (contain(FASTCPD_FAMILIES, family)) {
     cost_hessian_wrapper = std::bind(  // # nocov start
-      &Fastcpd::cost_update_hessian,  // # nocov end
+      get_hessian,  // # nocov end
       this,
       std::placeholders::_1,
       std::placeholders::_2,
@@ -335,103 +335,6 @@ double Fastcpd::get_cval_sen(
   return cval;
 }
 
-mat Fastcpd::get_hessian_arma(
-  const unsigned int segment_start,
-  const unsigned int segment_end,
-  const colvec& theta
-) {
-  const mat data_segment = data.rows(segment_start, segment_end);
-  const unsigned int segment_length = segment_end - segment_start + 1;
-  // TODO(doccstat): Maybe we can store all these computations
-  mat reversed_data = reverse(data_segment, 0);
-  colvec reversed_theta = reverse(theta);
-  if (segment_length < max(order) + 1) {
-    return eye(theta.n_elem, theta.n_elem);
-  }
-  colvec variance_term = zeros(segment_length);
-  for (unsigned int i = max(order); i < segment_length; i++) {
-    variance_term(i) = data_segment(i, 0) - dot(
-        reversed_theta.rows(order(1) + 1, sum(order)),
-        data_segment.rows(i - order(0), i - 1)
-      ) - dot(
-        reversed_theta.rows(1, order(1)),
-        variance_term.rows(i - order(1), i - 1)
-      );
-  }
-  colvec reversed_variance_term = reverse(variance_term);
-  mat phi_coefficient = zeros(segment_length, order(0)),
-      psi_coefficient = zeros(segment_length, order(1));
-  for (unsigned int i = max(order); i < segment_length; i++) {
-    phi_coefficient.row(i) = -reversed_data.rows(
-      segment_length - i, segment_length - i + order(0) - 1
-    ).t() - reversed_theta.rows(1, order(1)).t() *
-    phi_coefficient.rows(i - order(1), i - 1);
-  }
-  for (unsigned int i = order(1); i < segment_length; i++) {
-    psi_coefficient.row(i) = -reversed_variance_term.rows(
-        segment_length - i, segment_length - i + order(1) - 1
-      ).t() - reversed_theta.rows(1, order(1)).t() *
-      psi_coefficient.rows(i - order(1), i - 1);
-  }
-  mat reversed_coef_phi = reverse(phi_coefficient, 0),
-      reversed_coef_psi = reverse(psi_coefficient, 0);
-  cube phi_psi_coefficient = zeros(order(1), order(0), segment_length),
-        psi_psi_coefficient = zeros(order(1), order(1), segment_length);
-  for (unsigned int i = order(1); i < segment_length; i++) {
-    mat phi_psi_coefficient_part = zeros(order(1), order(0)),
-        psi_psi_coefficient_part = zeros(order(1), order(1));
-    for (unsigned int j = 1; j <= order(1); j++) {
-      phi_psi_coefficient_part +=
-        phi_psi_coefficient.slice(i - j) * theta(order(0) - 1 + j);
-    }
-    phi_psi_coefficient.slice(i) = -reversed_coef_phi.rows(
-      segment_length - i, segment_length - i + order(1) - 1
-    ) - phi_psi_coefficient_part;
-    for (unsigned int j = 1; j <= order(1); j++) {
-      psi_psi_coefficient_part +=
-        psi_psi_coefficient.slice(i - j) * theta(order(0) - 1 + j);
-    }
-    psi_psi_coefficient.slice(i) = -reversed_coef_psi.rows(
-        segment_length - i, segment_length - i + order(1) - 1
-      ) - reversed_coef_psi.rows(
-        segment_length - i, segment_length - i + order(1) - 1
-      ).t() - psi_psi_coefficient_part;
-  }
-  mat hessian = zeros(sum(order) + 1, sum(order) + 1);
-  hessian.submat(0, 0, order(0) - 1, order(0) - 1) =
-    phi_coefficient.row(segment_length - 1).t() *
-    phi_coefficient.row(segment_length - 1) / theta(sum(order));
-  hessian.submat(0, order(0), order(0) - 1, sum(order) - 1) = (
-    phi_psi_coefficient.slice(segment_length - 1).t() *
-      variance_term(segment_length - 1) +
-      phi_coefficient.row(segment_length - 1).t() *
-      psi_coefficient.row(segment_length - 1)
-  ) / theta(sum(order));
-  hessian.submat(order(0), 0, sum(order) - 1, order(0) - 1) =
-    hessian.submat(0, order(0), order(0) - 1, sum(order) - 1).t();
-  hessian.submat(0, sum(order), order(0) - 1, sum(order)) =
-    -phi_coefficient.row(segment_length - 1).t() *
-    variance_term(segment_length - 1) / theta(sum(order)) / theta(sum(order));
-  hessian.submat(sum(order), 0, sum(order), order(0) - 1) =
-    hessian.submat(0, sum(order), order(0) - 1, sum(order)).t();
-  hessian.submat(order(0), order(0), sum(order) - 1, sum(order) - 1) = (
-    psi_coefficient.row(segment_length - 1).t() *
-    psi_coefficient.row(segment_length - 1) +
-    psi_psi_coefficient.slice(segment_length - 1) *
-    variance_term(segment_length - 1)
-  ) / theta(sum(order));
-  hessian.submat(order(0), sum(order), sum(order) - 1, sum(order)) =
-    -psi_coefficient.row(segment_length - 1).t() *
-    variance_term(segment_length - 1) / theta(sum(order)) / theta(sum(order));
-  hessian.submat(sum(order), order(0), sum(order), sum(order) - 1) =
-    hessian.submat(order(0), sum(order), sum(order) - 1, sum(order)).t();
-  hessian(sum(order), sum(order)) =
-    std::pow(variance_term(segment_length - 1), 2) /
-    std::pow(theta(sum(order)), 3) -
-    1.0 / 2.0 / std::pow(theta(sum(order)), 2);
-  return hessian;
-}
-
 CostResult Fastcpd::get_optimized_cost(
   const unsigned int segment_start,
   const unsigned int segment_end
@@ -525,10 +428,10 @@ void Fastcpd::update_cost_parameters_step(
     );
     gradient = cost_gradient_result;
   } else {
-    hessian_i += cost_update_hessian(
+    hessian_i += (this->*get_hessian)(
       segment_start + data_start, segment_start + data_end, theta_hat.col(i)
     );
-    gradient = cost_update_gradient(
+    gradient = (this->*get_gradient)(
       segment_start + data_start, segment_start + data_end, theta_hat.col(i)
     );
   }
@@ -691,7 +594,7 @@ void Fastcpd::update_fastcpd_parameters(const unsigned int t) {
     const rowvec x = data.row(t - 1).tail(p);
     hessian_new = x.t() * x + epsilon * eye<mat>(p, p);
   } else if (family == "arma") {
-    hessian_new = cost_update_hessian(0, t - 1, coef_add.t());
+    hessian_new = (this->*get_hessian)(0, t - 1, coef_add.t());
   } else if (!contain(FASTCPD_FAMILIES, family)) {
     hessian_new = zeros<mat>(p, p);
   }
diff --git a/src/fastcpd_class_public.cc b/src/fastcpd_class_public.cc
index 198316f0..3fef51bf 100644
--- a/src/fastcpd_class_public.cc
+++ b/src/fastcpd_class_public.cc
@@ -71,180 +71,39 @@ Fastcpd::Fastcpd(
     "[:bar] :current/:total in :elapsed", data_n_rows
   );
 
-  create_cost_function_wrapper(cost);
-  create_cost_gradient_wrapper(cost_gradient);
-  create_cost_hessian_wrapper(cost_hessian);
-
-  // TODO(doccstat): Store environment functions from R.
-}
-
-colvec Fastcpd::cost_update_gradient(
-  const unsigned int segment_start,
-  const unsigned int segment_end,
-  const colvec& theta
-) {
-  const unsigned int segment_length = segment_end - segment_start + 1;
-  const mat data_segment = data.rows(segment_start, segment_end);
-  rowvec new_data = data_segment.row(segment_length - 1);
-  rowvec x = new_data.tail(new_data.n_elem - 1);
-  double y = new_data(0);
-  colvec gradient;
   if (family.compare("binomial") == 0) {
-    gradient = - (y - 1 / (1 + exp(-as_scalar(x * theta)))) * x.t();
+    get_gradient = &Fastcpd::get_gradient_binomial;
   } else if (family.compare("poisson") == 0) {
-    gradient = - (y - exp(as_scalar(x * theta))) * x.t();
+    get_gradient = &Fastcpd::get_gradient_poisson;
   } else if (family == "lasso" || family == "gaussian") {
-    gradient = - (y - as_scalar(x * theta)) * x.t();
+    get_gradient = &Fastcpd::get_gradient_lm;
   } else if (family == "arma" && order(0) == 0) {
-    const unsigned int q = order(1);
-    mat reversed_data = reverse(data_segment, 0);
-    colvec reversed_theta = reverse(theta);
-    if (segment_length < q + 1) {
-      return ones(theta.n_elem);
-    }
-    colvec variance_term = zeros(segment_length);
-    for (unsigned int i = q; i < segment_length; i++) {
-      variance_term(i) = data_segment(i, 0) -
-        dot(reversed_theta.rows(1, q), variance_term.rows(i - q, i - 1));
-    }
-    colvec reversed_variance_term = reverse(variance_term);
-    mat psi_coefficient = zeros(segment_length, q);
-    for (unsigned int i = q; i < segment_length; i++) {
-      psi_coefficient.row(i) = -reversed_variance_term.rows(
-          segment_length - i, segment_length - i + q - 1
-        ).t() - reversed_theta.rows(1, q).t() *
-        psi_coefficient.rows(i - q, i - 1);
-    }
-    gradient = zeros(q + 1);
-    gradient.rows(0, q - 1) =
-      psi_coefficient.row(segment_length - 1).t() *
-      variance_term(segment_length - 1) / theta(q);
-    gradient(q) = 1.0 / 2.0 / theta(q) -
-      std::pow(variance_term(segment_length - 1), 2) / 2.0 /
-      std::pow(theta(q), 2);
+    get_gradient = &Fastcpd::get_gradient_ma;
   } else if (family == "arma" && order(0) > 0) {
-    mat reversed_data = reverse(data_segment, 0);
-    colvec reversed_theta = reverse(theta);
-    if (segment_length < max(order) + 1) {
-      return ones(theta.n_elem);
-    }
-    colvec variance_term = zeros(segment_length);
-    for (unsigned int i = max(order); i < segment_length; i++) {
-      variance_term(i) = data_segment(i, 0) - dot(
-          reversed_theta.rows(order(1) + 1, sum(order)),
-          data_segment.rows(i - order(0), i - 1)
-        ) - dot(
-          reversed_theta.rows(1, order(1)),
-          variance_term.rows(i - order(1), i - 1)
-        );
-    }
-    colvec reversed_variance_term = reverse(variance_term);
-    mat phi_coefficient = zeros(segment_length, order(0)),
-        psi_coefficient = zeros(segment_length, order(1));
-    for (unsigned int i = max(order); i < segment_length; i++) {
-      phi_coefficient.row(i) = -reversed_data.rows(
-        segment_length - i, segment_length - i + order(0) - 1
-      ).t() - reversed_theta.rows(1, order(1)).t() *
-      phi_coefficient.rows(i - order(1), i - 1);
-    }
-    for (unsigned int i = order(1); i < segment_length; i++) {
-      psi_coefficient.row(i) = -reversed_variance_term.rows(
-          segment_length - i, segment_length - i + order(1) - 1
-        ).t() - reversed_theta.rows(1, order(1)).t() *
-        psi_coefficient.rows(i - order(1), i - 1);
-    }
-    gradient = zeros(sum(order) + 1);
-    gradient.rows(0, order(0) - 1) =
-      phi_coefficient.row(segment_length - 1).t() *
-      variance_term(segment_length - 1) / theta(sum(order));
-    gradient.rows(order(0), sum(order) - 1) =
-      psi_coefficient.row(segment_length - 1).t() *
-      variance_term(segment_length - 1) / theta(sum(order));
-    gradient(sum(order)) = 1.0 / 2.0 / theta(sum(order)) -
-      std::pow(variance_term(segment_length - 1), 2) / 2.0 /
-      std::pow(theta(sum(order)), 2);
+    get_gradient = &Fastcpd::get_gradient_arma;
   } else {
-    // # nocov start
-    stop("This branch should not be reached at functions.cc: 188.");
-    // # nocov end
+    ERROR("Invalid family.");
   }
-  return gradient;
-}
 
-mat Fastcpd::cost_update_hessian(
-  const unsigned int segment_start,
-  const unsigned int segment_end,
-  const colvec& theta
-) {
-  const mat data_segment = data.rows(segment_start, segment_end);
-  const unsigned int segment_length = segment_end - segment_start + 1;
-  rowvec new_data = data_segment.row(segment_length - 1);
-  rowvec x = new_data.tail(new_data.n_elem - 1);
-  mat hessian;
   if (family.compare("binomial") == 0) {
-    double prob = 1 / (1 + exp(-as_scalar(x * theta)));
-    hessian = (x.t() * x) * as_scalar((1 - prob) * prob);
+    get_hessian = &Fastcpd::get_hessian_binomial;
   } else if (family.compare("poisson") == 0) {
-    double prob = exp(as_scalar(x * theta));
-    // Prevent numerical issues if `prob` is too large.
-    hessian = (x.t() * x) * std::min(as_scalar(prob), 1e10);
+    get_hessian = &Fastcpd::get_hessian_poisson;
   } else if (family == "lasso" || family == "gaussian") {
-    hessian = x.t() * x;
+    get_hessian = &Fastcpd::get_hessian_lm;
   } else if (family == "arma" && order(0) == 0) {
-    const unsigned int q = order(1);
-    // TODO(doccstat): Maybe we can store all these computations
-    mat reversed_data = reverse(data_segment, 0);
-    colvec reversed_theta = reverse(theta);
-    if (segment_length < q + 1) {
-      return eye(theta.n_elem, theta.n_elem);
-    }
-    colvec variance_term = zeros(segment_length);
-    for (unsigned int i = q; i < segment_length; i++) {
-      variance_term(i) = data_segment(i, 0) - dot(
-        reversed_theta.rows(1, q), variance_term.rows(i - q, i - 1)
-      );
-    }
-    colvec reversed_variance_term = reverse(variance_term);
-    mat psi_coefficient = zeros(segment_length, q);
-    for (unsigned int i = q; i < segment_length; i++) {
-      psi_coefficient.row(i) = -reversed_variance_term.rows(
-          segment_length - i, segment_length - i + q - 1
-        ).t() - reversed_theta.rows(1, q).t() *
-        psi_coefficient.rows(i - q, i - 1);
-    }
-    mat reversed_coef_psi = reverse(psi_coefficient, 0);
-    cube psi_psi_coefficient = zeros(q, q, segment_length);
-    for (unsigned int i = q; i < segment_length; i++) {
-      mat psi_psi_coefficient_part = zeros(q, q);
-      for (unsigned int j = 1; j <= q; j++) {
-        psi_psi_coefficient_part +=
-          psi_psi_coefficient.slice(i - j) * theta(j - 1);
-      }
-      psi_psi_coefficient.slice(i) = -reversed_coef_psi.rows(
-          segment_length - i, segment_length - i + q - 1
-        ) - reversed_coef_psi.rows(
-          segment_length - i, segment_length - i + q - 1
-        ).t() - psi_psi_coefficient_part;
-    }
-    hessian = zeros(q + 1, q + 1);
-    hessian.submat(0, 0, q - 1, q - 1) = (
-      psi_coefficient.row(segment_length - 1).t() *
-      psi_coefficient.row(segment_length - 1) +
-      psi_psi_coefficient.slice(segment_length - 1) *
-      variance_term(segment_length - 1)
-    ) / theta(q);
-    hessian.submat(0, q, q - 1, q) =
-      -psi_coefficient.row(segment_length - 1).t() *
-      variance_term(segment_length - 1) / theta(q) / theta(q);
-    hessian.submat(q, 0, q, q - 1) = hessian.submat(0, q, q - 1, q).t();
-    hessian(q, q) =
-      std::pow(variance_term(segment_length - 1), 2) /
-      std::pow(theta(q), 3) -
-      1.0 / 2.0 / std::pow(theta(q), 2);
+    get_hessian = &Fastcpd::get_hessian_ma;
   } else if (family == "arma" && order(0) > 0) {
-    return get_hessian_arma(segment_start, segment_end, theta);
+    get_hessian = &Fastcpd::get_hessian_arma;
+  } else {
+    ERROR("Invalid family.");
   }
-  return hessian;
+
+  create_cost_function_wrapper(cost);
+  create_cost_gradient_wrapper(cost_gradient);
+  create_cost_hessian_wrapper(cost_hessian);
+
+  // TODO(doccstat): Store environment functions from R.
 }
 
 void Fastcpd::create_theta_sum(
diff --git a/src/fastcpd_classes.h b/src/fastcpd_classes.h
index 1ef18e62..6cf414a2 100644
--- a/src/fastcpd_classes.h
+++ b/src/fastcpd_classes.h
@@ -93,32 +93,6 @@ class Fastcpd {
     const bool warm_start
   );
 
-  // Function to calculate the gradient at the current data.
-  //
-  // @param data A data frame containing the data to be segmented.
-  // @param theta Estimated theta from the previous iteration.
-  // @param family Family of the model.
-  // @param order Order of the time series models.
-  //
-  // @return Gradient at the current data.
-  colvec cost_update_gradient(
-    const unsigned int segment_start,
-    const unsigned int segment_end,
-    const colvec& theta
-  );
-
-  // Function to calculate the Hessian matrix at the current data.
-  //
-  // @param data A data frame containing the data to be segmented.
-  // @param theta Estimated theta from the previous iteration.
-  //
-  // @return Hessian at the current data.
-  mat cost_update_hessian(
-    const unsigned int segment_start,
-    const unsigned int segment_end,
-    const colvec& theta
-  );
-
   // Set \code{theta_sum} for a specific column.
   void create_theta_sum(const unsigned int col, colvec new_theta_sum);
 
@@ -145,6 +119,32 @@ class Fastcpd {
   // Update \code{theta_sum} for a specific column by adding to that column.
   void update_theta_sum(const unsigned int col, colvec new_theta_sum);
 
+  // Function to calculate the gradient at the current data.
+  //
+  // @param data A data frame containing the data to be segmented.
+  // @param theta Estimated theta from the previous iteration.
+  // @param family Family of the model.
+  // @param order Order of the time series models.
+  //
+  // @return Gradient at the current data.
+  colvec (Fastcpd::*get_gradient)(
+    const unsigned int segment_start,
+    const unsigned int segment_end,
+    const colvec& theta
+  );
+
+  // Function to calculate the Hessian matrix at the current data.
+  //
+  // @param data A data frame containing the data to be segmented.
+  // @param theta Estimated theta from the previous iteration.
+  //
+  // @return Hessian at the current data.
+  mat (Fastcpd::*get_hessian)(
+    const unsigned int segment_start,
+    const unsigned int segment_end,
+    const colvec& theta
+  );
+
  private:
   // `act_num` is used in Lasso and Gaussian families only.
   colvec act_num;
@@ -339,12 +339,66 @@ class Fastcpd {
     const double lambda
   );
 
+  colvec get_gradient_arma(
+    const unsigned int segment_start,
+    const unsigned int segment_end,
+    const colvec& theta
+  );
+
+  colvec get_gradient_binomial(
+    const unsigned int segment_start,
+    const unsigned int segment_end,
+    const colvec& theta
+  );
+
+  colvec get_gradient_lm(
+    const unsigned int segment_start,
+    const unsigned int segment_end,
+    const colvec& theta
+  );
+
+  colvec get_gradient_ma(
+    const unsigned int segment_start,
+    const unsigned int segment_end,
+    const colvec& theta
+  );
+
+  colvec get_gradient_poisson(
+    const unsigned int segment_start,
+    const unsigned int segment_end,
+    const colvec& theta
+  );
+
   mat get_hessian_arma(
     const unsigned int segment_start,
     const unsigned int segment_end,
     const colvec& theta
   );
 
+  mat get_hessian_binomial(
+    const unsigned int segment_start,
+    const unsigned int segment_end,
+    const colvec& theta
+  );
+
+  mat get_hessian_lm(
+    const unsigned int segment_start,
+    const unsigned int segment_end,
+    const colvec& theta
+  );
+
+  mat get_hessian_ma(
+    const unsigned int segment_start,
+    const unsigned int segment_end,
+    const colvec& theta
+  );
+
+  mat get_hessian_poisson(
+    const unsigned int segment_start,
+    const unsigned int segment_end,
+    const colvec& theta
+  );
+
   CostResult get_nll_arma(
     const unsigned int segment_start,
     const unsigned int segment_end
diff --git a/src/fastcpd_types.h b/src/fastcpd_types.h
index dc737dda..22b61b11 100644
--- a/src/fastcpd_types.h
+++ b/src/fastcpd_types.h
@@ -5,10 +5,15 @@
 
 #include "RcppArmadillo.h"
 
+#define ERROR(msg) \
+  Rcout << "error: " << __FILE__ << ":" << __LINE__ << ": " << msg << std::endl
 #ifdef DEBUG
-#  define DEBUG_RCOUT(x) Rcout << #x << ": " << x << std::endl
+#define DEBUG_RCOUT(x) Rcout << #x << ": " << x << std::endl
+#define INFO(msg) \
+  Rcout << "info: " << __FILE__ << ":" << __LINE__ << ": " << msg << std::endl
 #else
-#  define DEBUG_RCOUT(x) do {} while (0)
+#define DEBUG_RCOUT(x) do {} while (0)
+#define INFO(msg) do {} while (0)
 #endif
 
 using ::arma::abs;
diff --git a/src/test-functions.cc b/src/test-functions.cc
index 27162e17..6f395f19 100644
--- a/src/test-functions.cc
+++ b/src/test-functions.cc
@@ -95,7 +95,7 @@ context("get_nll_wo_cv Unit Test") {
   }
 }
 
-context("cost_update_gradient Unit Test") {
+context("get_gradient Unit Test") {
   test_that("arma(3, 2) is correct for 200 data points") {
     fastcpd::classes::Fastcpd fastcpd_class(
       /* beta */ 0,
@@ -127,14 +127,15 @@ context("cost_update_gradient Unit Test") {
     );
 
     const colvec theta = 0.1 * ones<colvec>(6);
-    const colvec gradient = fastcpd_class.cost_update_gradient(0, 199, theta);
+    const colvec gradient =
+      (fastcpd_class.*fastcpd_class.get_gradient)(0, 199, theta);
     const colvec expected_gradient =
       {4.401258, 6.600128, -7.591818, 4.151778, 7.503752, -2.806806};
     expect_true(norm(gradient - expected_gradient, "fro") < 1e-6);
   }
 }
 
-context("cost_update_hessian Unit Test") {
+context("get_hessian Unit Test") {
   test_that("binomal is correct for a two dimensional data") {
     fastcpd::classes::Fastcpd fastcpd_class(
       /* beta */ 0,
@@ -166,7 +167,7 @@ context("cost_update_hessian Unit Test") {
     );
 
     const colvec theta = {-0.5, 0.3};
-    const mat hessian = fastcpd_class.cost_update_hessian(0, 0, theta);
+    const mat hessian = (fastcpd_class.*fastcpd_class.get_hessian)(0, 0, theta);
     const mat expected_hessian =
         {{0.238'28, 0.047'656}, {0.047'656, 0.009'531'2}};
     expect_true(norm(hessian - expected_hessian, "fro") < 0.000'001);
@@ -203,7 +204,7 @@ context("cost_update_hessian Unit Test") {
     );
 
     const colvec theta = {-0.5, 0.3};
-    const mat hessian = fastcpd_class.cost_update_hessian(0, 0, theta);
+    const mat hessian = (fastcpd_class.*fastcpd_class.get_hessian)(0, 0, theta);
     const mat expected_hessian =
         {{0.644'036'4, 0.128'807}, {0.128'807, 0.025'761'6}};
     expect_true(norm(hessian - expected_hessian, "fro") < 0.000'001);
@@ -240,7 +241,8 @@ context("cost_update_hessian Unit Test") {
     );
 
     const colvec theta = 0.1 * ones<colvec>(6);
-    const mat hessian = fastcpd_class.cost_update_hessian(0, 199, theta);
+    const mat hessian =
+      (fastcpd_class.*fastcpd_class.get_hessian)(0, 199, theta);
     const mat expected_hessian = {
       { 12.406525,  18.60483, -21.40027,   4.743794,  28.98263, -44.01258},
       { 18.604831,  27.89981, -32.09185,  25.380851,  27.48253, -66.00128},