diff --git a/src/initialization.jl b/src/initialization.jl
index d2f6ddf..c1f15db 100644
--- a/src/initialization.jl
+++ b/src/initialization.jl
@@ -1,7 +1,8 @@
 """
 ```
-function one_draw(loglikelihood::Function, parameters::ParameterVector{U},
-                  data::Matrix{Float64}; regime_switching::Bool = false) where {U<:Number}
+one_draw(loglikelihood::Function, parameters::ParameterVector{U},
+         data::Matrix{Float64}; regime_switching::Bool = false,
+         toggle::Bool = true) where {U<:Number}
 ```
 
 Finds and returns one valid draw from parameter distribution, along with its
@@ -10,11 +11,20 @@ log likelihood and log posterior.
 Set `regime_switching` to true if
 there are regime-switching parameters. Otherwise, not all the values of the
 regimes will be used or saved.
+
+Set `toggle` to false if, after calculating the loglikelihood,
+the values in the fields of every parameter in `parameters`
+are set to their regime 1 values. The regime-switching version of `rand`
+requires that the fields of all parameters take their regime 1 values,
+or else sampling may be wrong. The default is `true` as a safety, but
+if speed is a paramount concern, setting `toggle = true` will avoid
+unnecessary computations.
 """
 function one_draw(loglikelihood::Function, parameters::ParameterVector{U},
-                  data::Matrix{Float64}; regime_switching::Bool = false) where {U<:Number}
+                  data::Matrix{Float64}; regime_switching::Bool = false,
+                  toggle::Bool = true) where {U<:Number}
     success    = false
-    draw       = vec(rand(parameters, 1, regime_switching = regime_switching))
+    draw       = vec(rand(parameters, 1, regime_switching = regime_switching, toggle = toggle))
 
     draw_loglh = draw_logpost = 0.0
 
@@ -22,7 +32,12 @@ function one_draw(loglikelihood::Function, parameters::ParameterVector{U},
         try
             update!(parameters, draw)
 
-            draw_loglh   = loglikelihood(parameters, data)
+            draw_loglh = loglikelihood(parameters, data)
+
+            if toggle
+                toggle_regime!(parameters, 1)
+            end
+
             draw_logpost = prior(parameters)
 
             if (draw_loglh == -Inf) || (draw_loglh === NaN)
@@ -39,7 +54,7 @@ function one_draw(loglikelihood::Function, parameters::ParameterVector{U},
         end
 
         if any(isinf.(draw_loglh))
-            draw = vec(rand(parameters, 1, regime_switching = regime_switching))
+            draw = vec(rand(parameters, 1, regime_switching = regime_switching, toggle = false))
         else
             success = true
         end
@@ -51,7 +66,7 @@ end
 ```
 function initial_draw!(loglikelihood::Function, parameters::ParameterVector{U},
                        data::Matrix{Float64}, c::Cloud; parallel::Bool = false,
-                       regime_switching::Bool = false) where {U<:Number}
+                       regime_switching::Bool = false, toggle::Bool = true) where {U<:Number}
 ```
 
 Draw from a general starting distribution (set by default to be from the prior) to
@@ -61,10 +76,18 @@ particle objects in the particle cloud in place.
 Set `regime_switching` to true if
 there are regime-switching parameters. Otherwise, not all the values of the
 regimes will be used or saved.
+
+Set `toggle` to false if, after calculating the loglikelihood,
+the values in the fields of every parameter in `parameters`
+are set to their regime 1 values. The regime-switching version of `rand`
+requires that the fields of all parameters take their regime 1 values,
+or else sampling may be wrong. The default is `true` as a safety, but
+if speed is a paramount concern, setting `toggle = true` will avoid
+unnecessary computations.
 """
 function initial_draw!(loglikelihood::Function, parameters::ParameterVector{U},
                        data::Matrix{Float64}, c::Cloud; parallel::Bool = false,
-                       regime_switching::Bool = false) where {U<:Number}
+                       regime_switching::Bool = false, toggle::Bool = true) where {U<:Number}
     n_parts = length(c)
 
     # ================== Define closure on one_draw function ==================
@@ -72,8 +95,8 @@ function initial_draw!(loglikelihood::Function, parameters::ParameterVector{U},
     sendto(workers(), parameters = parameters)
     sendto(workers(), data       = data)
 
-    one_draw_closure() = one_draw(loglikelihood, parameters, data, regime_switching = regime_switching)
-    @everywhere one_draw_closure() = one_draw(loglikelihood, parameters, data, regime_switching = regime_switching)
+    one_draw_closure() = one_draw(loglikelihood, parameters, data, regime_switching = regime_switching, toggle = toggle)
+    @everywhere one_draw_closure() = one_draw(loglikelihood, parameters, data, regime_switching = regime_switching, toggle = toggle)
     # =========================================================================
 
     # For each particle, finds valid parameter draw and returns loglikelihood & posterior
@@ -97,24 +120,30 @@ end
 
 """
 ```
-function draw_likelihood(loglikelihood::Function, parameters::ParameterVector{U},
-                         data::Matrix{Float64}, draw::Vector{Float64}) where {U<:Number}
+draw_likelihood(loglikelihood::Function, parameters::ParameterVector{U},
+                data::Matrix{Float64}, draw::Vector{Float64};
+                toggle::Bool = true) where {U<:Number}
 ```
 Computes likelihood of a particular parameter draw; returns loglh and logpost.
 """
 function draw_likelihood(loglikelihood::Function, parameters::ParameterVector{U},
-                         data::Matrix{Float64}, draw::Vector{Float64}) where {U<:Number}
+                         data::Matrix{Float64}, draw::Vector{Float64};
+                         toggle::Bool = true) where {U<:Number}
     update!(parameters, draw)
     loglh   = loglikelihood(parameters, data)
+    if toggle
+        ModelConstructors.toggle_regime!(parameters, 1)
+    end
     logpost = prior(parameters)
     return scalar_reshape(loglh, logpost)
 end
 
 """
 ```
-function initialize_likelihoods!(loglikelihood::Function, parameters::ParameterVector{U},
-                                 data::Matrix{Float64}, c::Cloud;
-                                 parallel::Bool = false) where {U<:Number}
+initialize_likelihoods!(loglikelihood::Function, parameters::ParameterVector{U},
+                        data::Matrix{Float64}, c::Cloud;
+                        parallel::Bool = false,
+                        toggle::Bool = true) where {U<:Number}
 ```
 This function is made for transfering the log-likelihood values saved in the
 Cloud from a previous estimation to each particle's respective old_loglh
@@ -123,7 +152,8 @@ here is just the argument, data.
 """
 function initialize_likelihoods!(loglikelihood::Function, parameters::ParameterVector{U},
                                  data::Matrix{Float64}, c::Cloud;
-                                 parallel::Bool = false) where {U<:Number}
+                                 parallel::Bool = false,
+                                 toggle::Bool = true) where {U<:Number}
     n_parts = length(c)
     draws   = get_vals(c; transpose = false)
 
@@ -136,10 +166,10 @@ function initialize_likelihoods!(loglikelihood::Function, parameters::ParameterV
     sendto(workers(), data = data)
 
     draw_likelihood_closure(draw::Vector{Float64}) = draw_likelihood(loglikelihood, parameters,
-                                                                     data, draw)
+                                                                     data, draw, toggle = toggle)
     @everywhere draw_likelihood_closure(draw::Vector{Float64}) = draw_likelihood(loglikelihood,
                                                                                  parameters,
-                                                                                 data, draw)
+                                                                                 data, draw, toggle = toggle)
     # =========================================================================
 
     # TODO: handle when the likelihood with new data cannot be evaluated (returns -Inf),
diff --git a/src/mutation.jl b/src/mutation.jl
index 5c3b771..98caa8d 100755
--- a/src/mutation.jl
+++ b/src/mutation.jl
@@ -36,6 +36,10 @@ Execute one proposed move of the Metropolis-Hastings algorithm for a given param
 - `regime_switching::Bool`: Set to true if
     there are regime-switching parameters. Otherwise, not all the values of the
     regimes will be used or saved.
+- `toggle::Bool = true`: Flag for resetting the fields of parameter values to regime 1 anytime
+    the loglikelihood is computed. The regime-switching version of SMC assumes at various points
+    that this resetting occurs. If speed is important, then ensure that the fields of parameters
+    take their regime 1 values at the end of the loglikelihood computation and set `toggle = false`.
 
 ### Outputs:
 - `p::Vector{Float64}`: An updated particle containing updated parameter values,
@@ -48,7 +52,8 @@ function mutation(loglikelihood::Function, parameters::ParameterVector{U},
                   blocks_free::Vector{Vector{Int}}, blocks_all::Vector{Vector{Int}},
                   ϕ_n::S, ϕ_n1::S; c::S = 1., α::S = 1., n_mh_steps::Int = 1,
                   old_data::T = T(undef, size(data, 1), 0),
-                  regime_switching::Bool = false) where {S<:AbstractFloat,T<:AbstractMatrix, U<:Number}
+                  regime_switching::Bool = false,
+                  toggle::Bool = true) where {S<:AbstractFloat,T<:AbstractMatrix, U<:Number}
 
     step_prob   = rand() # Draw initial step probability
 
@@ -94,12 +99,20 @@ function mutation(loglikelihood::Function, parameters::ParameterVector{U},
                 prior_new = prior(parameters)
                 like_new  = loglikelihood(parameters, data) #, regime_switching = regime_switching)
 
+                if toggle
+                    toggle_regime!(parameters, 1)
+                end
+
                 if like_new == -Inf
                     prior_new = like_old_data = -Inf
                 end
 
                 like_old_data = isempty(old_data) ? 0. : loglikelihood(parameters, old_data)
 
+                if toggle && isempty(old_data)
+                    toggle_regime!(parameters, 1)
+                end
+
             catch err
                 if isa(err, ParamBoundsError) || isa(err, LinearAlgebra.LAPACKException) ||
                    isa(err, PosDefException)  || isa(err, SingularException)             ||
diff --git a/src/smc_main.jl b/src/smc_main.jl
index 4daf526..9bc0f96 100644
--- a/src/smc_main.jl
+++ b/src/smc_main.jl
@@ -74,6 +74,10 @@ function smc(loglikelihood::Function, parameters::ParameterVector{U}, data::Matr
     the Hessian is negative definite).
 - `regime_switching::Bool = false`: Flag if there are regime-switching parameters. Otherwise, not all the values of the
     regimes will be used or saved.
+- `toggle::Bool = true`: Flag for resetting the fields of parameter values to regime 1 anytime
+    the loglikelihood is computed. The regime-switching version of SMC assumes at various points
+    that this resetting occurs. If speed is important, then ensure that the fields of parameters
+    take their regime 1 values at the end of the loglikelihood computation and set `toggle = false`.
 
 ### Outputs
 
@@ -141,7 +145,8 @@ function smc(loglikelihood::Function, parameters::ParameterVector{U}, data::Matr
              continue_intermediate::Bool = false,
              intermediate_stage_start::Int = 0,
              tempered_update_prior_weight::S = 0.0,
-             regime_switching::Bool = false) where {S<:AbstractFloat, U<:Number}
+             regime_switching::Bool = false,
+             toggle::Bool = true) where {S<:AbstractFloat, U<:Number}
 
     ########################################################################################
     ### Settings
@@ -157,16 +162,19 @@ function smc(loglikelihood::Function, parameters::ParameterVector{U}, data::Matr
                               blocks_free::Vector{Vector{Int64}}, blocks_all::Vector{Vector{Int64}},
                               ϕ_n::S, ϕ_n1::S; c::S = 1.0, α::S = 1.0, n_mh_steps::Int = 1,
                               old_data::T = Matrix{S}(undef, size(data, 1), 0),
-                              regime_switching::Bool = false) where {S<:Float64, T<:Matrix}
+                              regime_switching::Bool = false, toggle::Bool = true) where {S<:Float64, T<:Matrix}
         return mutation(loglikelihood, parameters, data, p, d_μ, d_Σ, n_free_para, blocks_free, blocks_all,
-                        ϕ_n, ϕ_n1; c = c, α = α, n_mh_steps = n_mh_steps, old_data = old_data, regime_switching = regime_switching)
+                        ϕ_n, ϕ_n1; c = c, α = α, n_mh_steps = n_mh_steps, old_data = old_data, regime_switching = regime_switching,
+                        toggle = toggle)
     end
     @everywhere function mutation_closure(p::Vector{S}, d_μ::Vector{S}, d_Σ::Matrix{S},
-             blocks_free::Vector{Vector{Int64}}, blocks_all::Vector{Vector{Int64}}, n_free_para::Int,
-             ϕ_n::S, ϕ_n1::S; c::S = 1.0, α::S = 1.0, n_mh_steps::Int = 1,
-             old_data::T = Matrix{S}(undef, size(data, 1), 0), regime_switching::Bool = false) where {S<:Float64, T<:Matrix}
+                                          blocks_free::Vector{Vector{Int64}}, blocks_all::Vector{Vector{Int64}}, n_free_para::Int,
+                                          ϕ_n::S, ϕ_n1::S; c::S = 1.0, α::S = 1.0, n_mh_steps::Int = 1,
+                                          old_data::T = Matrix{S}(undef, size(data, 1), 0), regime_switching::Bool = false,
+                                          toggle::Bool = true) where {S<:Float64, T<:Matrix}
         return mutation(loglikelihood, parameters, data, p, d_μ, d_Σ, blocks_free, blocks_all, n_free_para,
-                        ϕ_n, ϕ_n1; c = c, α = α, n_mh_steps = n_mh_steps, old_data = old_data, regime_switching = regime_switching)
+                        ϕ_n, ϕ_n1; c = c, α = α, n_mh_steps = n_mh_steps, old_data = old_data,
+                        regime_switching = regime_switching, toggle = toggle)
     end
 
     # Check that if there's a tempered update, old and current vintages are different
@@ -189,9 +197,7 @@ function smc(loglikelihood::Function, parameters::ParameterVector{U}, data::Matr
     for para in parameters
         if !isempty(para.regimes)
             for (ind, val) in para.regimes[:value]
-                if ind == 1
-                    nothing
-                else # ind != 1
+                if ind != 1
                     n_para_rs += 1
                 end
             end
@@ -224,7 +230,8 @@ function smc(loglikelihood::Function, parameters::ParameterVector{U}, data::Matr
             initialize_cloud_settings!(cloud; tempered_update = tempered_update,
                                        n_parts = n_parts, n_Φ = n_Φ, c = c, accept = target)
 
-            initialize_likelihoods!(loglikelihood, parameters, data, cloud; parallel = parallel)
+            initialize_likelihoods!(loglikelihood, parameters, data, cloud; parallel = parallel,
+                                    toggle = toggle)
 
         elseif (tempered_update_prior_weight > 0.0) || (old_n_parts != n_parts)
             # Resample from bridge distribution
@@ -242,7 +249,8 @@ function smc(loglikelihood::Function, parameters::ParameterVector{U}, data::Matr
             # Make a cloud by drawing from the prior
             if n_from_prior > 0
                 prior_cloud = Cloud(n_para, n_from_prior)
-                initial_draw!(loglikelihood, parameters, old_data, prior_cloud, parallel = parallel, regime_switching = regime_switching)
+                initial_draw!(loglikelihood, parameters, old_data, prior_cloud, parallel = parallel,
+                              regime_switching = regime_switching, toggle = toggle)
 
                 cloud = Cloud(n_para, n_to_resample + n_from_prior)
                 update_cloud!(cloud, vcat(bridge_cloud.particles, prior_cloud.particles))
@@ -253,7 +261,8 @@ function smc(loglikelihood::Function, parameters::ParameterVector{U}, data::Matr
                 cloud = bridge_cloud
             end
 
-            initialize_likelihoods!(loglikelihood, parameters, data, cloud; parallel = parallel)
+            initialize_likelihoods!(loglikelihood, parameters, data, cloud; parallel = parallel,
+                                    toggle = toggle)
 
             reset_weights!(cloud)
             cloud.resamples      += 1
@@ -266,7 +275,8 @@ function smc(loglikelihood::Function, parameters::ParameterVector{U}, data::Matr
         cloud = load(loadpath, "cloud")
     else
         # Instantiating Cloud object, update draws, loglh, & logpost
-        initial_draw!(loglikelihood, parameters, data, cloud; parallel = parallel, regime_switching = regime_switching)
+        initial_draw!(loglikelihood, parameters, data, cloud; parallel = parallel, regime_switching = regime_switching,
+                      toggle = toggle)
         initialize_cloud_settings!(cloud; tempered_update = tempered_update,
                                    n_parts = n_parts, n_Φ = n_Φ, c = c, accept = target)
     end
@@ -389,14 +399,16 @@ function smc(loglikelihood::Function, parameters::ParameterVector{U}, data::Matr
                 mutation_closure(cloud.particles[k, :], θ_bar_fr, R_fr, n_free_para,
                                  blocks_free, blocks_all, ϕ_n, ϕ_n1; c = c, α = α,
                                  n_mh_steps = n_mh_steps, old_data = old_data,
-                                 regime_switching = regime_switching)
+                                 regime_switching = regime_switching,
+                                 toggle = toggle)
             end
         else
             hcat([mutation_closure(cloud.particles[k, :], θ_bar_fr, R_fr, n_free_para,
                                    blocks_free, blocks_all, ϕ_n, ϕ_n1; c = c,
                                    α = α, n_mh_steps = n_mh_steps,
                                    old_data = old_data,
-                                   regime_switching = regime_switching) for k=1:n_parts]...)
+                                   regime_switching = regime_switching,
+                                   toggle = toggle) for k=1:n_parts]...)
         end
         update_cloud!(cloud, new_particles)
         update_acceptance_rate!(cloud)
diff --git a/test/modelsetup.jl b/test/modelsetup.jl
index c55d269..3e4f0bb 100644
--- a/test/modelsetup.jl
+++ b/test/modelsetup.jl
@@ -1,17 +1,21 @@
+using ModelConstructors, HDF5
+
+regenerate_data = false
+
 ###################################################################
 # Set Up Linear Model
 ###################################################################
 
-function setup_linear_model()
+function setup_linear_model(; regime_switching::Bool = false)
     m = GenericModel()
 
     # Set up linear parameters
     m <= parameter(:α1, 0., (-1e5, 1e5), (-1e5, 1e5), Untransformed(), Normal(0, 1e3),
                   fixed = false)
     m <= parameter(:β1, 0., (-1e5, 1e5), (-1e5, 1e5), Untransformed(), Normal(0, 1e3),
-                  fixed = false)  
+                  fixed = false)
     m <= parameter(:σ1, 1., (1e-5, 1e5), (1e-5, 1e5), SquareRoot(), Uniform(0, 1e3),
-                  fixed = false) 
+                  fixed = false)
     m <= parameter(:α2, 0., (-1e5, 1e5), (-1e5, 1e5), Untransformed(), Normal(0, 1e3),
                   fixed = false)
     m <= parameter(:β2, 0., (-1e5, 1e5), (-1e5, 1e5), Untransformed(), Normal(0, 1e3),
@@ -39,57 +43,117 @@ function setup_linear_model()
     m <= Setting(:resampling_threshold, .5)
     m <= Setting(:use_fixed_schedule, true)
 
+    if regime_switching
+        for i in 1:3
+            ModelConstructors.set_regime_val!(m[Symbol("α$(i)")], 1, 0.)
+            ModelConstructors.set_regime_fixed!(m[Symbol("α$(i)")], 1, false)
+            ModelConstructors.set_regime_val!(m[Symbol("α$(i)")], 2, 0.)
+            ModelConstructors.set_regime_fixed!(m[Symbol("α$(i)")], 2, false)
+            ModelConstructors.set_regime_val!(m[Symbol("α$(i)")], 3, float(i))
+            ModelConstructors.set_regime_fixed!(m[Symbol("α$(i)")], 3, true) # Do not estimate this parameter, just to check this functionality
+        end
+        for i in 1:3
+            ModelConstructors.set_regime_val!(m[Symbol("β$(i)")], 1, 0.)
+            ModelConstructors.set_regime_prior!(m[Symbol("β$(i)")], 1, Normal(0, 1e3)) # regime-switching prior, just to check functionality
+            ModelConstructors.set_regime_val!(m[Symbol("β$(i)")], 2, 0.)
+            ModelConstructors.set_regime_prior!(m[Symbol("β$(i)")], 2, Normal(0, 1e3))
+            ModelConstructors.set_regime_val!(m[Symbol("β$(i)")], 3, 0.)
+            ModelConstructors.set_regime_prior!(m[Symbol("β$(i)")], 3, Normal(0, 1e2))
+        end
+    end
+
     return m
 end
 
 
 # Generate Data
-#
-#  X = rand(Float64, (3, 100))
-#  err = rand(Float64, (3, 100))
-#  β = fill(2, 3) #(1,2,3)
-#  β[1]=1
-#  β[2]=2
-#  β[3]=3
-#
-#  α = fill(2,3)
-#  α[1]=1
-#  α[2]=2
-#  α[3]=3
-#
-#  data = β .* X .+ α .+ err
-
-# Save Data
-
-#  h5open("reference/test_data.h5", "w") do file
-#      write(file, "data", data)
-#      write(file, "X", X)
-#  end
+if regenerate_data
+    Random.seed!(1793)
+
+    X = randn(Float64, (3, 100))
+    err = randn(Float64, (3, 100))
+    β = collect(1:3)
+    α = collect(1:3)
+
+    data = β .* X .+ α .+ err
+
+    # Regime-switching version
+    β₁ = collect(1:3)
+    β₂ = collect(2:4)
+    β₃ = collect(3:5)
+
+    α₁ = collect(1:3)
+    α₂ = collect(1:3)
+
+    rsdata = similar(data)
+    rsdata[:, 1:50] = β₁ .* X[:, 1:50] .+ α₁ .+ err[:, 1:50]
+    rsdata[:, 51:70] = β₂ .* X[:, 51:70] .+ α₂ .+ err[:, 51:70]
+    rsdata[:, 71:end] = β₃ .* X[:, 71:end] .+ α₂ .+ err[:, 71:end]
+
+    # Save Data
+
+    h5open("reference/test_data.h5", "w") do file
+        write(file, "data", data)
+        write(file, "rsdata", rsdata)
+        write(file, "X", X)
+    end
+end
 
 # Read Data
 data = h5read("reference/test_data.h5", "data")
+rsdata = h5read("reference/test_data.h5", "rsdata")
 X = h5read("reference/test_data.h5", "X")
 
 # Log Likelihood Function
 N = 3
 function loglik_fn(p, d)
     # we assume the ordering of (α_i, β_i, σ_i)
-    Σ = zeros(N,N) 
+    Σ = zeros(N,N)
     α = Vector{Float64}(undef,N)
     β = Vector{Float64}(undef,N)
     for i in 1:N
-        α[i]   = p[i * 3 - 2]
-        β[i]   = p[i * 3 - 1]
-        Σ[i,i] = p[i * 3]^2
+        α[i]   = p[i * 3 - 2] # note p is a ParameterVector
+        β[i]   = p[i * 3 - 1] # but since LHS are Vector{Float64},
+        Σ[i,i] = p[i * 3]^2   # the value of p is autoamtically stored
     end
     det_Σ = det(Σ)
     inv_Σ = inv(Σ)
-    term1 = -size(d,2) / 2 * log(2 * π) - 1 /2 * log(det_Σ)
+    term1 = -N / 2 * log(2 * π) - 1 /2 * log(det_Σ) # N should be number of equations
     logprob = 0.
-    errors = d .- α .- β .* X[:, 1:size(sparse(d),2)] 
-    for t in 1:size(d,2) #see above
-        logprob += term1 - 1/2 * dot(errors, inv_Σ * errors)
+    errors = d .- α .- β .* X[:, 1:size(d, 2)] # bridging test uses first 1/2 of sample so need to be robust there
+    for t in 1:size(d, 2) # Number of time periods
+        logprob += term1 - 1/2 * dot(errors[:, t], inv_Σ * errors[:, t]) # Should be error in each time period
     end
     return logprob
 end
 
+function rs_loglik_fn(p, d)
+    # we assume the ordering of (α_i, β_i, σ_i)
+    Σ = zeros(N,N)
+    α = Vector{Float64}(undef, 3 * N)
+    β = Vector{Float64}(undef, 3 * N)
+    for i in 1:N
+        α[i]   = regime_val(p[N * (i - 1) + 1], 1)
+        β[i]   = regime_val(p[N * (i - 1) + 2], 1)
+        Σ[i,i] = p[N * (i - 1) + 2]
+    end
+    for i in 1:N # Looping over each entry of the constant vector and betas in a given regime
+        α[N + i] = regime_val(p[N * (i - 1) + 1], 2) # α1, α2, α3 each has 3 regimes. This line is for regime 2.
+        α[(2 * N) + i] = regime_val(p[N * (i - 1) + 1], 3) # This line is for regime 3.
+        β[N + i] = regime_val(p[N * (i - 1) + 2], 2) # β1, β2, β3 each has 3 regimes.
+        β[(2 * N) + i] = regime_val(p[N * (i - 1) + 2], 3)
+    end
+
+    det_Σ = det(Σ)
+    inv_Σ = inv(Σ)
+    term1 = -N / 2 * log(2 * π) - 1 /2 * log(det_Σ)
+    logprob = 0.
+    errors = similar(d)
+    errors[:, 1:50] = d[:, 1:50] .- α[1:N] .- β[1:N] .* X[:, 1:50]
+    errors[:, 51:70] = d[:, 51:70] .- α[(N + 1):(2 * N)] .- β[(N + 1):(2 * N)] .* X[:, 51:70]
+    errors[:, 71:end] = d[:, 71:end] .- α[(2 * N + 1):(3 * N)] .- β[(2 * N + 1):(3 * N)] .* X[:, 71:end]
+    for t in 1:size(d, 2) # see above
+        logprob += term1 - 1/2 * dot(errors[:, t], inv_Σ * errors[:, t])
+    end
+    return logprob
+end
diff --git a/test/reference/smc_cloud_fix=true_version=111.jld2 b/test/reference/smc_cloud_fix=true_version=111.jld2
index ef60004..7240e9f 100644
Binary files a/test/reference/smc_cloud_fix=true_version=111.jld2 and b/test/reference/smc_cloud_fix=true_version=111.jld2 differ
diff --git a/test/reference/smc_cloud_fix=true_version=150.jld2 b/test/reference/smc_cloud_fix=true_version=150.jld2
index fc908f0..a7544a9 100644
Binary files a/test/reference/smc_cloud_fix=true_version=150.jld2 and b/test/reference/smc_cloud_fix=true_version=150.jld2 differ
diff --git a/test/reference/test_data.h5 b/test/reference/test_data.h5
index 89b4a08..7bb0b38 100644
Binary files a/test/reference/test_data.h5 and b/test/reference/test_data.h5 differ
diff --git a/test/runtests.jl b/test/runtests.jl
index f99a8eb..2170c3a 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -7,13 +7,13 @@ import SparseArrays.SparseMatrixCSC
 import ModelConstructors.Setting
 
 my_tests = [
-            "smc",
             "helpers",
             "initialization",
             "resample",
             "util",
             "mutation",
-            "particle"
+            "particle",
+            "smc"
             ]
 
 for test in my_tests
@@ -21,4 +21,3 @@ for test in my_tests
     @printf " * %s\n" test_file
     include(test_file)
 end
-
diff --git a/test/smc.jl b/test/smc.jl
index 1fd44cf..4f40f53 100644
--- a/test/smc.jl
+++ b/test/smc.jl
@@ -2,7 +2,7 @@ using ModelConstructors, HDF5, Random, JLD2, FileIO, SMC, Test
 include("modelsetup.jl")
 
 path = dirname(@__FILE__)
-writing_output = false
+writing_output = true
 
 if VERSION < v"1.5"
     ver = "111"
@@ -21,14 +21,15 @@ data = h5read("reference/test_data.h5", "data")
 
 @everywhere Random.seed!(42)
 
+println("Estimating Linear Model... (approx. 3 minutes)")
 
-println("Estimating Linear Model... (approx. 2 minutes)")
-
-SMC.smc(loglik_fn, m.parameters, data, verbose = :none, use_fixed_schedule = true, parallel = true, n_Φ = 100, n_mh_steps = 1, resampling_method = :polyalgo, data_vintage = "200707", target = 0.25, savepath = savepath, particle_store_path = particle_store_path, α = .9, threshold_ratio = .5, smc_iteration = 0)
+SMC.smc(loglik_fn, m.parameters, data, verbose = :none, use_fixed_schedule = true,
+        parallel = true, n_Φ = 120, n_mh_steps = 1, resampling_method = :polyalgo,
+        data_vintage = "200707", target = 0.25, savepath = savepath,
+        particle_store_path = particle_store_path, α = .9, threshold_ratio = .5, smc_iteration = 0)
 
 println("Estimation done!")
 
-
 test_file = load(rawpath(m, "estimate", "smc_cloud.jld2"))
 test_cloud  = test_file["cloud"]
 test_w      = test_file["w"]
@@ -49,6 +50,12 @@ saved_W     = saved_file["W"]
 
 ####################################################################
 cloud_fields = fieldnames(typeof(test_cloud))
+@testset "Linear Regression Parameter Estimates Are Close" begin
+    # Posterior mean should be close to true parameters
+    @test maximum(abs.(mean(SMC.get_vals(test_cloud), dims = 2) -
+                       [1., 1., 1., 2., 2., 1., 3., 3., 1.])) < .5
+end
+
 @testset "ParticleCloud Fields: Linear" begin
     @test @test_matrix_approx_eq SMC.get_vals(test_cloud) SMC.get_vals(saved_cloud)
     @test @test_matrix_approx_eq SMC.get_loglh(test_cloud) SMC.get_loglh(saved_cloud)
@@ -80,7 +87,7 @@ end
     @test @test_matrix_approx_eq test_W saved_W
 end
 
-
+# TODO: add a check here that the correct parameters are estimated
 
 
 ####################################################################
@@ -99,7 +106,12 @@ data = h5read("reference/test_data.h5", "data")
 # Estimate with 1st half of sample
 m_old = deepcopy(m)
 
-SMC.smc(loglik_fn, m_old.parameters, data[:, 1:Int(floor(end/2))], verbose = :none, use_fixed_schedule = true, parallel = true,  n_Φ = 100, n_mh_steps = 1, resampling_method = :polyalgo, data_vintage = "000000", target = 0.25, savepath = savepath, particle_store_path = particle_store_path, α = .9, threshold_ratio = .5, smc_iteration = 0, n_parts = 1000)
+println("Estimating Linear Model on 1st half of sample... (approx. 2 minutes)")
+SMC.smc(loglik_fn, m_old.parameters, data[:, 1:Int(floor(end/2))], verbose = :none,
+        use_fixed_schedule = true, parallel = true,  n_Φ = 100, n_mh_steps = 1,
+        resampling_method = :polyalgo, data_vintage = "000000", target = 0.25,
+        savepath = savepath, particle_store_path = particle_store_path, α = .9,
+        threshold_ratio = .5, smc_iteration = 0, n_parts = 1000)
 
 m_new = deepcopy(m)
 
@@ -113,12 +125,21 @@ loadpath = rawpath(m_old, "estimate", "smc_cloud.jld2")
 
 old_cloud = load(loadpath, "cloud")
 
-SMC.smc(loglik_fn, m_new.parameters, data, old_data=data[:,1:Int(floor(end/2))], old_cloud=old_cloud, verbose = :none, use_fixed_schedule = true, parallel = true,  n_Φ = 100, n_mh_steps = 1, resampling_method = :polyalgo, data_vintage = "200708", target = 0.25, savepath = savepath, particle_store_path = particle_store_path, α = .9, threshold_ratio = .5, smc_iteration = 0)
+println("Estimating Linear Model using a bridge distribution... (approx. 2 minutes)")
+SMC.smc(loglik_fn, m_new.parameters, data, old_data=data[:,1:Int(floor(end/2))], old_cloud=old_cloud,
+        verbose = :none, use_fixed_schedule = true, parallel = true,
+        n_Φ = 100, n_mh_steps = 1, resampling_method = :polyalgo, data_vintage = "200708",
+        target = 0.25, savepath = savepath, particle_store_path = particle_store_path,
+        α = .9, threshold_ratio = .5, smc_iteration = 0)
 
 loadpath = rawpath(m_new, "estimate", "smc_cloud.jld2")
 new_cloud = load(loadpath, "cloud")
- 
+
+@testset "Linear Regression Parameter Estimates Are Close" begin
+    # Posterior mean should be close to true parameters
+    @test maximum(abs.(mean(SMC.get_vals(test_cloud), dims = 2) - [1., 1., 1., 2., 2., 1., 3., 3., 1.])) < .25
+end
+
 # Clean output files up
 rm(rawpath(m_new, "estimate", "smc_cloud.jld2"))
 rm(rawpath(m_new, "estimate", "smcsave.h5"))
-