diff --git a/src/utility.cpp b/src/utility.cpp
index becd918..76f6b6e 100644
--- a/src/utility.cpp
+++ b/src/utility.cpp
@@ -111,6 +111,24 @@ vector <double> calcLLKs(const vector <double> &refCount,
     return tmpLLKs;
 }
 
+vector <log_double_t> calcSiteLikelihoods(const vector <double> &refCount,
+                                          const vector <double> &altCount,
+                                          const vector <double> &expectedWsaf,
+                                          size_t firstIndex, size_t length,
+                                          double fac, double err) {
+    assert(expectedWsaf.size() == length);
+    vector <log_double_t> siteLikelihoods(length);
+    size_t index = firstIndex;
+    for (size_t i = 0; i < length; i++) {
+        assert(expectedWsaf[i] >= 0);
+        // assert (expectedWsaf[i] <= 1);
+        siteLikelihoods[i] = calcSiteLikelihood(refCount[index], altCount[index],
+                                                expectedWsaf[i], err, fac);
+        index++;
+    }
+    return siteLikelihoods;
+}
+
 
 log_double_t calcSiteLikelihood(double ref, double alt, double unadjustedWsaf, double err,
     double fac) {
diff --git a/src/utility.hpp b/src/utility.hpp
index 85d69f3..b5b6303 100644
--- a/src/utility.hpp
+++ b/src/utility.hpp
@@ -58,17 +58,34 @@ vector <T> vecFromTo(const vector<T> &vec, size_t start, size_t end) {
 template <typename T>
 vector <T> vecSum(const vector<T> &vecA, const vector<T> &vecB) {
     assert(vecA.size() == vecB.size());
-    vector <T> tmpSum(vecA.size(), (T)0);
+    vector <T> tmpSum(vecA.size());
     for (size_t i = 0; i < vecA.size(); i++) {
         tmpSum[i] = vecA[i] + vecB[i];
     }
     return tmpSum;
 }
 
+template <typename T>
+vector <T> vecSum(const T &A, const vector<T> &vecB) {
+    vector <T> tmpSum(vecB.size());
+    for (size_t i = 0; i < vecB.size(); i++) {
+        tmpSum[i] = A + vecB[i];
+    }
+    return tmpSum;
+}
+
+template <typename T>
+vector <T> vecSum(const vector<T> &vecA, const T &B) {
+    vector <T> tmpSum(vecA.size());
+    for (size_t i = 0; i < vecA.size(); i++) {
+        tmpSum[i] = vecA[i] + B;
+    }
+    return tmpSum;
+}
+
 
 template <typename T>
-vector <T> vecProd(const vector<T> &vecA,
-                   const vector<T> &vecB) {
+vector <T> vecProd(const vector<T> &vecA, const vector<T> &vecB) {
     assert(vecA.size() == vecB.size());
     vector <T> tmpProd(vecA.size());
     for (size_t i = 0; i < vecA.size(); i++) {
@@ -77,6 +94,24 @@ vector <T> vecProd(const vector<T> &vecA,
     return tmpProd;
 }
 
+template <typename T>
+vector <T> vecProd(const T &A, const vector<T> &vecB) {
+    vector <T> tmpProd(vecB.size());
+    for (size_t i = 0; i < vecB.size(); i++) {
+        tmpProd[i] = A * vecB[i];
+    }
+    return tmpProd;
+}
+
+template <typename T>
+vector <T> vecProd(const vector<T> &vecA, const T& B) {
+    vector <T> tmpProd(vecA.size());
+    for (size_t i = 0; i < vecA.size(); i++) {
+        tmpProd[i] = vecA[i] * B;
+    }
+    return tmpProd;
+}
+
 
 template <typename T>
 T sumOfVec(const vector <T>& array ) {
@@ -87,6 +122,15 @@ T sumOfVec(const vector <T>& array ) {
     return tmp;
 }
 
+template <typename T>
+T sum(const std::vector <T>& array ) {
+    T tmp = 0;
+    for (auto const& value : array) {
+        tmp += value;
+    }
+    return tmp;
+}
+
 template <typename T>
 T product(const std::vector<T>& xs)
 {
@@ -143,6 +187,10 @@ vector <double> calcLLKs(const vector <double> &refCount,
     const vector <double> &altCount,
     const vector <double> &expectedWsaf, size_t firstIndex, size_t length,
     double fac, double err = 0.01);
+vector <log_double_t> calcSiteLikelihoods(const vector <double> &refCount,
+                                          const vector <double> &altCount,
+                                          const vector <double> &expectedWsaf, size_t firstIndex, size_t length,
+                                          double fac, double err = 0.01);
 log_double_t calcSiteLikelihood(double ref, double alt,
                                 double unadjustedWsaf, double err, double fac);
 size_t sampleIndexGivenProp(RandomGenerator* rg, vector <double> proportion);