diff --git a/cppbugs/distributions/mcmc.exponential.censored.hpp b/cppbugs/distributions/mcmc.exponential.censored.hpp
deleted file mode 100644
index 3e41c01..0000000
--- a/cppbugs/distributions/mcmc.exponential.censored.hpp
+++ /dev/null
@@ -1,51 +0,0 @@
-///////////////////////////////////////////////////////////////////////////
-// Copyright (C) 2011 Jacques-Henri Jourdan                              //
-//                                                                       //
-// This program is free software: you can redistribute it and/or modify  //
-// it under the terms of the GNU General Public License as published by  //
-// the Free Software Foundation, either version 3 of the License, or     //
-// (at your option) any later version.                                   //
-//                                                                       //
-// This program is distributed in the hope that it will be useful,       //
-// but WITHOUT ANY WARRANTY; without even the implied warranty of        //
-// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         //
-// GNU General Public License for more details.                          //
-//                                                                       //
-// You should have received a copy of the GNU General Public License     //
-// along with this program.  If not, see <http://www.gnu.org/licenses/>. //
-///////////////////////////////////////////////////////////////////////////
-
-#pragma once
-
-
-#include <cmath>
-#include <armadillo>
-#include <cppbugs/mcmc.stochastic.hpp>
-#include <cstdio>
-namespace cppbugs {
-
-  template <typename T,typename U, typename V>
-  class ExponentialCensoredLikelihiood : public Likelihiood {
-    const T& x_;
-    const U& lambda_;
-    const V& delta_;
-  public:
-    ExponentialCensoredLikelihiood(const T& x, const U& lambda, const V& delta): x_(x), lambda_(lambda), delta_(delta) { dimension_check(x_, lambda_, delta_); }
-    inline double calc() const {
-      return arma::accu(arma::schur(delta_, log_approx(lambda_)) - arma::schur(lambda_, x_));
-    }
-  };
-
-  template<typename T>
-  class ObservedExponentialCensored : public Observed<T> {
-  public:
-    ObservedExponentialCensored(const T& value): Observed<T>(value) {}
-
-    template<typename U, typename V>
-    ObservedExponentialCensored<T>& dexpcens(const U& lambda, const V& delta) {
-      Stochastic::likelihood_functor = new ExponentialCensoredLikelihiood<T,U,V>(Observed<T>::value,lambda,delta);
-      return *this;
-    }
-  };
-
-} // namespace cppbugs
diff --git a/cppbugs/mcmc.arma.extensions.hpp b/cppbugs/mcmc.arma.extensions.hpp
index 317ac6a..f265e66 100644
--- a/cppbugs/mcmc.arma.extensions.hpp
+++ b/cppbugs/mcmc.arma.extensions.hpp
@@ -121,94 +121,6 @@ namespace arma {
     return eOpCube<T1, eop_factln>(A.get_ref());
   }
 
-  // cube
-  //! element-wise multiplication of BaseCube objects with same element type
-  template<typename T1, typename T2>
-  arma_inline
-  const eGlueCube<T1, T2, eglue_schur>
-  schur
-  (
-   const BaseCube<typename T1::elem_type,T1>& X,
-   const BaseCube<typename T1::elem_type,T2>& Y
-   )
-  {
-    arma_extra_debug_sigprint();
-    return eGlueCube<T1, T2, eglue_schur>(X.get_ref(), Y.get_ref());
-  }
-
-  //! element-wise multiplication of BaseCube objects with different element types
-  template<typename T1, typename T2>
-  inline
-  const mtGlueCube<typename promote_type<typename T1::elem_type, typename T2::elem_type>::result, T1, T2, glue_mixed_schur>
-  schur
-  (
-   const BaseCube< typename force_different_type<typename T1::elem_type, typename T2::elem_type>::T1_result, T1>& X,
-   const BaseCube< typename force_different_type<typename T1::elem_type, typename T2::elem_type>::T2_result, T2>& Y
-   )
-  {
-    arma_extra_debug_sigprint();
-    typedef typename T1::elem_type eT1;
-    typedef typename T2::elem_type eT2;
-    typedef typename promote_type<eT1,eT2>::result out_eT;
-    promote_type<eT1,eT2>::check();
-    return mtGlueCube<out_eT, T1, T2, glue_mixed_schur>( X.get_ref(), Y.get_ref() );
-  }
-
-  // matrix
-  template<typename T1, typename T2>
-  arma_inline
-  const eGlue<T1, T2, eglue_schur>
-  schur(const Base<typename T1::elem_type,T1>& X, const Base<typename T1::elem_type,T2>& Y) {
-    arma_extra_debug_sigprint();
-    return eGlue<T1, T2, eglue_schur>(X.get_ref(), Y.get_ref());
-  }
-
-  //! element-wise multiplication of Base objects with different element types
-  template<typename T1, typename T2>
-  inline
-  const mtGlue<typename promote_type<typename T1::elem_type, typename T2::elem_type>::result, T1, T2, glue_mixed_schur>
-  schur
-  (
-   const Base< typename force_different_type<typename T1::elem_type, typename T2::elem_type>::T1_result, T1>& X,
-   const Base< typename force_different_type<typename T1::elem_type, typename T2::elem_type>::T2_result, T2>& Y
-   )
-  {
-    arma_extra_debug_sigprint();
-    typedef typename T1::elem_type eT1;
-    typedef typename T2::elem_type eT2;
-    typedef typename promote_type<eT1,eT2>::result out_eT;
-    promote_type<eT1,eT2>::check();
-    return mtGlue<out_eT, T1, T2, glue_mixed_schur>( X.get_ref(), Y.get_ref() );
-  }
-
-
-  //! Base * scalar
-  template<typename T1>
-  arma_inline
-  const eOp<T1, eop_scalar_times>
-  schur
-  (const Base<typename T1::elem_type,T1>& X, const typename T1::elem_type k)
-  {
-    arma_extra_debug_sigprint();
-    return eOp<T1, eop_scalar_times>(X.get_ref(),k);
-  }
-
-  //! scalar * Base
-  template<typename T1>
-  arma_inline
-  const eOp<T1, eop_scalar_times>
-  schur
-  (const typename T1::elem_type k, const Base<typename T1::elem_type,T1>& X)
-  {
-    arma_extra_debug_sigprint();
-    return eOp<T1, eop_scalar_times>(X.get_ref(),k);  // NOTE: order is swapped
-  }
-
-  double schur(const int x, const double y) { return x * y; }
-  double schur(const double x, const int y) { return x * y; }
-  double schur(const double& x, const double& y) { return x * y; }
-  double schur(const int& x, const int& y) { return x * y; }
-
   // insert an 'any' function for bools into the arma namespace
   bool any(const bool x) {
     return x;
diff --git a/cppbugs/mcmc.math.hpp b/cppbugs/mcmc.math.hpp
index 75fae22..b595440 100644
--- a/cppbugs/mcmc.math.hpp
+++ b/cppbugs/mcmc.math.hpp
@@ -25,6 +25,37 @@
 // Stochastic/Math related functions
 namespace cppbugs {
 
+  template<typename T, typename U>
+  inline
+  auto schur_product(T&& x, U&& y) ->
+    decltype (arma::operator%(std::forward<T>(x), std::forward<U>(y))){
+    return arma::operator%(std::forward<T>(x), std::forward<U>(y));
+  }
+
+  template<typename T>
+  inline
+  auto schur_product(const typename T::elem_type& x, T&& y) ->
+    decltype (arma::operator*(x, std::forward<T>(y))){
+    return arma::operator*(x, std::forward<T>(y));
+  }
+
+  template<typename T>
+  inline
+  auto schur_product(T&& x, const typename T::elem_type& y) ->
+    decltype (arma::operator*(std::forward<T>(x), y)){
+    return arma::operator*(std::forward<T>(x), y);
+  }
+
+  double schur_product(const int x, const double y) { return x * y; }
+  double schur_product(const double x, const int y) { return x * y; }
+  double schur_product(const double x, const double y) { return x * y; }
+  double schur_product(const float x, const double y) { return x * y; }
+  double schur_product(const double x, const float y) { return x * y; }
+  float schur_product(const int x, const float y) { return x * y; }
+  float schur_product(const float x, const int y) { return x * y; }
+  float schur_product(const float x, const float y) { return x * y; }
+  int schur_product(const int x, const int y) { return x * y; }
+
   static inline double square(double x) {
     return x*x;
   }
@@ -55,7 +86,7 @@ namespace cppbugs {
 
   template<typename T, typename U, typename V>
   double normal_logp(const T& x, const U& mu, const V& tau) {
-    return arma::accu(0.5*log_approx(0.5*tau/arma::math::pi()) - 0.5 * arma::schur(tau, square(x - mu)));
+    return arma::accu(0.5*log_approx(0.5*tau/arma::math::pi()) - 0.5 * schur_product(tau, square(x - mu)));
   }
 
   template<typename T, typename U, typename V>
@@ -67,7 +98,7 @@ namespace cppbugs {
   double gamma_logp(const T& x, const U& alpha, const V& beta) {
     return arma::any(arma::vectorise(x < 0)) ?
       -std::numeric_limits<double>::infinity() :
-      arma::accu(arma::schur((alpha - 1.0),log_approx(x)) - arma::schur(beta,x) - lgamma(alpha) + arma::schur(alpha,log_approx(beta)));
+      arma::accu(schur_product((alpha - 1.0),log_approx(x)) - schur_product(beta,x) - lgamma(alpha) + schur_product(alpha,log_approx(beta)));
   }
 
   template<typename T, typename U, typename V>
@@ -75,7 +106,7 @@ namespace cppbugs {
     const double one = 1.0;
     return arma::any(arma::vectorise(x <= 0)) || arma::any(arma::vectorise(x >= 1)) || arma::any(arma::vectorise(alpha <= 0)) || arma::any(arma::vectorise(beta <= 0)) ?
       -std::numeric_limits<double>::infinity() :
-      arma::accu(lgamma(alpha+beta) - lgamma(alpha) - lgamma(beta) + arma::schur((alpha-one),log_approx(x)) + arma::schur((beta-one),log_approx(one-x)));
+      arma::accu(lgamma(alpha+beta) - lgamma(alpha) - lgamma(beta) + schur_product((alpha-one),log_approx(x)) + schur_product((beta-one),log_approx(one-x)));
   }
 
   double categorical_logp(const arma::ivec& x, const arma::mat& p) {
@@ -111,7 +142,7 @@ namespace cppbugs {
     if(arma::any(arma::vectorise(p <= 0)) || arma::any(arma::vectorise(p >= 1)) || arma::any(arma::vectorise(x < 0))  || arma::any(arma::vectorise(x > n))) {
       return -std::numeric_limits<double>::infinity();
     }
-    return arma::accu(arma::schur(x,log_approx(p)) + arma::schur((n-x),log_approx(1-p)) + arma::factln(n) - arma::factln(x) - arma::factln(n-x));
+    return arma::accu(schur_product(x,log_approx(p)) + schur_product((n-x),log_approx(1-p)) + arma::factln(n) - arma::factln(x) - arma::factln(n-x));
   }
 
   template<typename T, typename U>
@@ -119,7 +150,7 @@ namespace cppbugs {
     if( arma::any(arma::vectorise(p <= 0)) || arma::any(arma::vectorise(p >= 1)) || arma::any(arma::vectorise(x < 0))  || arma::any(arma::vectorise(x > 1)) ) {
       return -std::numeric_limits<double>::infinity();
     } else {
-      return arma::accu(arma::schur(x,log_approx(p)) + arma::schur((1-x), log_approx(1-p)));
+      return arma::accu(schur_product(x,log_approx(p)) + schur_product((1-x), log_approx(1-p)));
     }
   }
 
@@ -134,7 +165,7 @@ namespace cppbugs {
 
   template<typename T, typename U>
   double exponential_logp(const T& x, const U& lambda) {
-    return arma::accu(log_approx(lambda) - arma::schur(lambda, x));
+    return arma::accu(log_approx(lambda) - schur_product(lambda, x));
   }
 
   template<typename T, typename U>