diff --git a/src/particle_filters/particle_smoother_RS-FFBSi.jl b/src/particle_filters/particle_smoother_RS-FFBSi.jl
index 2f725a885fed08708ddb68bf3b6dc91a9955fc23..4e5a46419b701f3ba8c73836420935985d988c32 100644
--- a/src/particle_filters/particle_smoother_RS-FFBSi.jl
+++ b/src/particle_filters/particle_smoother_RS-FFBSi.jl
@@ -1,4 +1,5 @@
 function pf_smoother_FFBSi( w,xp, M, f_density )
+    error("Not yet implemented")
     (N,T) = size(w)
     xb = Array(Float64,M,T)
     wexp = exp(w)
diff --git a/src/particle_filters/pf_PMMH.jl b/src/particle_filters/pf_PMMH.jl
new file mode 100644
index 0000000000000000000000000000000000000000..b87b5f03ce16e7b07f234f20173b017174f51bc5
--- /dev/null
+++ b/src/particle_filters/pf_PMMH.jl
@@ -0,0 +1,68 @@
+using Debug
+function pf_PMMH(y, N, R, pf, f_theta, draw_theta, f_density_theta, g_density_theta, sigmaW0, theta0 )
+
+    n_params = length(theta0)
+    T = length(y)
+    xp = Array(Float64,(N,T))
+    w = Array(Float64,(N,T))
+    wnn = Array(Float64,(N,T))
+    Zf(wnn) = sum(log(mean(exp(wnn),1)))
+    Z = Array(Float64,R)
+
+    theta = Array(Float64,n_params,R)
+    theta[:,1] = theta0
+
+    g_density(y_t, x_t,w)  = g_density_theta(y_t, x_t,w, theta0)
+    f(x,t) = f_theta(x,t,theta0)
+    pf(xp, w, wnn, y, N, g_density, f, sigmaW0)
+
+    Z[1] = Zf(wnn)
+
+    breakPoints = 10
+
+    #     noiseMultiplier = 0.1^(1/R)
+    Nstart = N
+    breakTime = floor(Int64, R / breakPoints)
+    particleIncrease = 10*N
+    thetaP = draw_theta(theta0)
+    f(x,t) = f_theta(x,t,thetaP)
+    g_density(y_t, x_t,w)  = g_density_theta(y_t, x_t,w, thetaP)
+
+    for r = 2:R
+        thetaP = draw_theta(theta[:,r-1])
+
+        #         f_density(x_t1, x_t, t) = f_density_theta(x_t1, x_t, t, thetaP)
+
+
+        pf(xp, w, wnn, y, N, g_density, f, sigmaW0)
+        Zp = Zf(wnn)
+#         @show exp(Zp-Z[r-1])
+        Iaccept = rand() < exp(Zp-Z[r-1])
+
+        if Iaccept
+            theta[:,r] = thetaP
+            Z[r] = Zp
+        else
+            theta[:,r] = theta[:,r-1]
+            Z[r] = Z[r-1]
+        end
+
+        if (r % breakTime == 0)
+            display(r/R)
+#             N = round(Int64, particleIncrease*(r/R)^3 + Nstart  )
+#             xp = Array(Float64,(N,T))
+#             w = Array(Float64,(N,T))
+#             wnn = Array(Float64,(N,T))
+        end
+
+
+    end
+
+    return theta
+end
+
+
+
+
+
+
diff --git a/src/particle_filters/resample.jl b/src/particle_filters/resample.jl
index a72fa33ab45bb6e9b2d91048e84f2304f1cd7e1c..293b90ce95d8eebfe67499c72d5d403ac0a57854 100644
--- a/src/particle_filters/resample.jl
+++ b/src/particle_filters/resample.jl
@@ -22,7 +22,7 @@ function resample_systematic!(j,w,M)
     s = collect((rand()/M+0):1/M:bins[end])
     bo = 1
     for i = 1:M
-        @inbounds for b = bo:N
+        for b = bo:N
             if s[i] < bins[b]
                 j[i] = b
                 bo = b
@@ -51,6 +51,24 @@ function resample_systematic_exp(w,M)
     return j
 end
 
+function resample_systematic_exp!(j,w)
+    N = size(w,1)
+    M = size(j,1)
+    bins = cumsum(w) # devec
+    s = collect((rand()/M+0):1/M:bins[end])
+    bo = 1
+    for i = 1:M
+        for b = bo:N
+            if s[i] <= bins[b]
+                j[i] = b
+                bo = b
+                break
+            end
+        end
+    end
+    return j
+end
+
 # """
 # There is probably lots of room for improvement here. All bins need not be formed in the beginning.
 # One only has to keep 1 values, the current upper limit, no array needed.
diff --git a/src/sobol.jl b/src/sobol.jl
index 64b455dfd7a8429f9310652cbf14b0f65b07328e..5ab9290a27ad9d14665fd4f366fbd2c6aab02182 100644
--- a/src/sobol.jl
+++ b/src/sobol.jl
@@ -44,9 +44,6 @@ function sobol(X, nextseed, MeM)
     return X, nextseed, MeM
 end
 
-## Functions below are inner functions, I have not determined if this is needed, the fastest method or anything else
-#  ------------------------------------------------------------------------
-
 function getNewSobolVector(dim, seed, MeM)
 
     seed::Int64 = max(floor(seed),0)
@@ -95,7 +92,6 @@ function getNewSobolVector(dim, seed, MeM)
 
     return qrv, nextseed, MeM
 end
-# -------------------------------------------------------------------------
 
 function smallest0bit(b)
     i = 0
@@ -111,7 +107,7 @@ function smallest0bit(b)
     return i
 end
 
-#   ----------------------------------------------------------------------
+
 function InitSobol(dim::Int64)
     MeM = Mem(dim)
     MeM.seed_save = -1