diff --git a/examples/ou_integration_test.jl b/examples/ou_integration_test.jl
index 5025499..aa4c4aa 100644
--- a/examples/ou_integration_test.jl
+++ b/examples/ou_integration_test.jl
@@ -29,7 +29,7 @@ u0 = [copy(zeros(dimensions)) .+ randn(dimensions)*0.71 for i in 1:nensemble]
 k = 0.3
 ϵ = 0.002
 metric(y) = y[1]
-t_end = 100.0
+t_end = 50.0
 t_start = 0.0
 tspan = (t_start, t_end)
 # averaging window
@@ -37,7 +37,7 @@ T = 50
 steps_per_window = Int64.(round(T/dt))
 
 
-# Run GKLT -  bootstrap iterations - and compute necessary statistics per iteration
+# Run GKLT
 iters = 5
 em = zeros(nensemble*iters)
 em = reshape(em, (nensemble,iters))
@@ -60,7 +60,7 @@ for iter in 1:iters
     for i in 1:nensemble
         trajectory = [metric(sim.ensemble[i][j]) for j in 1:length(sim.ensemble[i])]
         lr_vector[i] = likelihood_ratio(sim, trajectory, λ, t_end-t_start)
-        # This just counts the number of (treated as independent) block means per bin, over all trajectories, over all iterations
+        # This just counts the number of (treated as independent) block means per bin, over all trajectories, over all iterations. Not doing anything with it, yet
         counts .+=  count_events_per_bin(block_mean(trajectory, steps_per_window), lr_vector[i],a_range)
         # This is for estimating return times via the Ragone paper method; via `return_curve`.
         a_m[i] = maximum(block_mean(trajectory, steps_per_window))
@@ -73,24 +73,31 @@ rmprocs([p for p in procs() if p != myid()])
 
 
 # Evolve the solution directly; compute statistics
-tspan_direct = (0.0, 5e5)
+tspan_direct = (0.0, 2.5e5)
 direct_solution = evolve_stochastic_system(model, zeros(dimensions), tspan_direct, dt, (maxiters = 1e8,));
 direct_cost = (tspan_direct[2]-tspan_direct[1])*dt
 u = [metric(direct_solution[k]) for k in 1:length(direct_solution)]
 A = block_mean(u, steps_per_window)
-m = 100
 em_direct, r_direct, r_paper_direct, σr_direct = return_curve(A, FT(T),  FT.(ones(length(A))))
 
 # Fit GEV
 m = 100 
 blocks = block_maxima(A, m)
 params, nll = fit_gev(blocks, [std(blocks), mean(blocks), 0.1])
-cdf_A= (gev_cdf.(a_range, params[1], params[2], params[3])).^(1.0/m)
-pdf_A= (cdf_A[2:end] .- cdf_A[1:end-1])/bin
+chain = evolve_mcmc(blocks, params, [0.01,0.01,0.05], log_likelihood_gev,1000000; save_at = 100)
+n_samples = 100
+chain_samples = chain[Int.(round.(rand(n_samples)*length(chain)))]
+gev_rtn = zeros(n_samples, length(a_range))
+for (i, sample) in enumerate(chain_samples)
+    cdf = min.(gev_cdf.(a_range, sample[1], sample[2], sample[3]).^(1.0/m), 1.0-eps(FT))
+    gev_rtn[i,:] .= 50.0 ./ (1.0 .- cdf)
+end
+median_curve = [median(gev_rtn[:,k]) for k in 1:length(a_range)]
+upper_curve = [sort(gev_rtn[:,k])[83] for k in 1:length(a_range)]
+lower_curve = [sort(gev_rtn[:,k])[17] for k in 1:length(a_range)]
 
 
-# Make plots
-plot1 = plot()
+# GKLT curve - we need to combine iterations into single result
 mean_a = zeros(length(a_range)-1)
 mean_rtn = zeros(length(a_range)-1)
 std_rtn = zeros(length(a_range)-1)
@@ -98,23 +105,24 @@ for i in 1:(length(a_range)-1)
     mask = (em[:] .< a_range[i+1]) .& ( em[:] .>= a_range[i])
     if sum(mask) >0
         mean_a[i] = mean(em[:][mask])
-        mean_rtn[i] = log10.(mean(r_paper[:][mask]))
-        std_rtn[i] = std(r_paper[:][mask])/mean(r_paper[:][mask])./log(10.0)
+        mean_rtn[i] = mean(r_paper[:][mask])
+        std_rtn[i] = std(r_paper[:][mask])
     end
 end
 nonzero = (mean_a .!= 0.0) .& (isnan.(std_rtn) .== 0)
 final_a = mean_a[nonzero]
 final_r = mean_rtn[nonzero]
 final_σr = std_rtn[nonzero]
-plot!(plot1,final_a, final_r, ribbon = final_σr, label = "GKLT Algorithm; using maxes")
-denominator = 1.0 .- cumsum(counts)/sum(counts)
-safe_denominator = denominator .> 10.0*eps(FT)
-plot!(plot1,a_range[safe_denominator], log10.(50.0 ./denominator[safe_denominator]), label = "GKLT Algorithm; using counts")
-plot!(plot1,em_direct, log10.(r_paper_direct), ribbon = σr_direct ./ r_paper_direct ./ log(10.0), label = "Direct Integration")
-plot!(plot1, a_range, log10.(50.0 ./ (1.0 .- cdf_A)), label = "GEV")
-plot!(plot1, a_range, log10.(50.0 ./ (1.0 .- 0.5 .*(1 .+erf.(a_range ./2 .^0.5 ./std(A))))), label = "Gaussian")
+
+plot1 = plot()
+plot!(plot1,final_a, final_r, ribbon = final_σr, label = "GKLT Algorithm")
+plot!(plot1,em_direct[2:end], r_paper_direct[2:end], ribbon = σr_direct[2:end], label = "Direct Integration")
+plot!(plot1, a_range, median_curve, ribbon = (median_curve .- lower_curve, upper_curve .- median_curve), label = "GEV", color = :purple)
 plot!(plot1,xlabel = "Event magnitude")
-plot!(plot1,ylabel = "Log10(Return Time)")
+plot!(plot1,ylabel = "Return Time")
 plot!(plot1, legend = :topleft)
-plot!(plot1, xrange = [0.0,1.0])
+plot!(plot1, xrange = [0.0,0.8])
+plot!(plot1, yaxis = :log)
+plot!(plot1, yrange = (1,1e10), yticks = [1, 1e2, 1e4, 1e6, 1e8, 1e10])
 savefig(plot1, "return_curve_ou.png")
+@assert direct_cost ≈ algorithm_cost
diff --git a/src/RareEvents.jl b/src/RareEvents.jl
index 3d9cb1a..63fa66b 100644
--- a/src/RareEvents.jl
+++ b/src/RareEvents.jl
@@ -254,6 +254,7 @@ end
 
 include("utils.jl")
 include("gev_utils.jl")
+include("mcmc_utils.jl")
 
 # I'm keeping this around for now; we use it in tests.
 function orig_sample_and_rewrite_history!(ensemble::Vector, frequencies::Array, idx_current::Int,rng)
diff --git a/src/mcmc_utils.jl b/src/mcmc_utils.jl
new file mode 100644
index 0000000..954ae0f
--- /dev/null
+++ b/src/mcmc_utils.jl
@@ -0,0 +1,33 @@
+
+function accept_or_reject(x::Vector, x_proposed::Vector, data::Vector, log_likelihood::Function)
+    u = rand()
+    if log(u) <= log_likelihood(data, x_proposed) - log_likelihood(data, x)
+        return x_proposed
+    else
+        return x
+    end
+end
+
+function mcmc_sample(data, x::Vector, k::Int, σ_k::Float64, log_likelihood::Function)
+    x_proposed = copy(x)
+    x_proposed[k] += randn()*σ_k
+    x_next = accept_or_reject(x, x_proposed, data, log_likelihood)
+    return vcat(x_next, log_likelihood(data, x_next))
+end
+
+function evolve_mcmc(data::Vector, x0::Vector, σ::Vector, log_likelihood::Function, nsteps::Int; save_at::Int)
+    nparams = length(x0)
+    chain = [zeros(nparams+1) for _ in 1:(div(nsteps,save_at)+1)]
+    chain[1] .= vcat(x0, log_likelihood(data, x0))
+    for i in 2:nsteps
+        k = Int64(ceil(rand()*3))
+        if i % save_at == 0
+            curr = div(i, save_at)+1
+            prev = curr -1
+            chain[curr] .= mcmc_sample(data, chain[prev][1:nparams],k, σ[k], log_likelihood)
+        end
+    end
+    return chain
+end
+
+export evolve_mcmc