diff --git a/REQUIRE b/REQUIRE
index 7650995..7755e4a 100644
--- a/REQUIRE
+++ b/REQUIRE
@@ -8,3 +8,4 @@ LightGraphs 0.7.0
 PDMats 0.3.6
 Showoff 0.0.6
 StatsBase 0.7.4
+ForwardDiff 0.5.0
diff --git a/src/Mamba.jl b/src/Mamba.jl
index 8d2a342..478a9c4 100644
--- a/src/Mamba.jl
+++ b/src/Mamba.jl
@@ -1,4 +1,5 @@
 using Distributions
+using ForwardDiff
 
 module Mamba
 
@@ -55,7 +56,7 @@ module Mamba
   #################### Variate Types ####################
 
   abstract type ScalarVariate <: Real end
-  abstract type ArrayVariate{N} <: DenseArray{Float64, N} end
+  abstract type ArrayVariate{N} <: DenseArray{Real, N} end
 
   const AbstractVariate = Union{ScalarVariate, ArrayVariate}
   const VectorVariate = ArrayVariate{1}
@@ -72,7 +73,7 @@ module Mamba
   #################### Dependent Types ####################
 
   type ScalarLogical <: ScalarVariate
-    value::Float64
+    value::Real
     symbol::Symbol
     monitor::Vector{Int}
     eval::Function
@@ -81,7 +82,7 @@ module Mamba
   end
 
   type ArrayLogical{N} <: ArrayVariate{N}
-    value::Array{Float64, N}
+    value::Array{Real, N}
     symbol::Symbol
     monitor::Vector{Int}
     eval::Function
@@ -90,7 +91,7 @@ module Mamba
   end
 
   type ScalarStochastic <: ScalarVariate
-    value::Float64
+    value::Real
     symbol::Symbol
     monitor::Vector{Int}
     eval::Function
@@ -100,7 +101,7 @@ module Mamba
   end
 
   type ArrayStochastic{N} <: ArrayVariate{N}
-    value::Array{Float64, N}
+    value::Array{Real, N}
     symbol::Symbol
     monitor::Vector{Int}
     eval::Function
@@ -127,7 +128,7 @@ module Mamba
   abstract type SamplerTune end
 
   type SamplerVariate{T<:SamplerTune} <: VectorVariate
-    value::Vector{Float64}
+    value::Vector{Real}
     tune::T
 
     function SamplerVariate{T}(x::AbstractVector, tune::T) where T<:SamplerTune
@@ -136,7 +137,7 @@ module Mamba
     end
 
     function SamplerVariate{T}(x::AbstractVector, pargs...; kargs...) where T<:SamplerTune
-      value = convert(Vector{Float64}, x)
+      value = convert(Vector{Real}, x)
       SamplerVariate{T}(value, T(value, pargs...; kargs...))
     end
   end
@@ -150,7 +151,7 @@ module Mamba
   end
 
   type ModelState
-    value::Vector{Float64}
+    value::Vector{Real}
     tune::Vector{Any}
   end
 
@@ -170,14 +171,14 @@ module Mamba
   abstract type AbstractChains end
 
   immutable Chains <: AbstractChains
-    value::Array{Float64, 3}
+    value::Array{Real, 3}
     range::Range{Int}
     names::Vector{AbstractString}
     chains::Vector{Int}
   end
 
   immutable ModelChains <: AbstractChains
-    value::Array{Float64, 3}
+    value::Array{Real, 3}
     range::Range{Int}
     names::Vector{AbstractString}
     chains::Vector{Int}
diff --git a/src/model/dependent.jl b/src/model/dependent.jl
index 75e60f3..fa77198 100644
--- a/src/model/dependent.jl
+++ b/src/model/dependent.jl
@@ -73,7 +73,7 @@ logpdf(d::AbstractDependent, x, transform::Bool=false) = 0.0
 #################### Constructors ####################
 
 function Logical(f::Function, monitor::Union{Bool, Vector{Int}}=true)
-  value = Float64(NaN)
+  value = Real(NaN)
   fx, src = modelfxsrc(depfxargs, f)
   l = ScalarLogical(value, :nothing, Int[], fx, src, Symbol[])
   setmonitor!(l, monitor)
@@ -81,7 +81,7 @@ end
 
 function Logical(d::Integer, f::Function,
                  monitor::Union{Bool, Vector{Int}}=true)
-  value = Array{Float64}(fill(0, d)...)
+  value = Array{Real}(fill(0, d)...)
   fx, src = modelfxsrc(depfxargs, f)
   l = ArrayLogical(value, :nothing, Int[], fx, src, Symbol[])
   setmonitor!(l, monitor)
@@ -135,7 +135,7 @@ end
 #################### Constructors ####################
 
 function Stochastic(f::Function, monitor::Union{Bool, Vector{Int}}=true)
-  value = Float64(NaN)
+  value = Real(NaN)
   fx, src = modelfxsrc(depfxargs, f)
   s = ScalarStochastic(value, :nothing, Int[], fx, src, Symbol[],
                        NullUnivariateDistribution())
@@ -144,7 +144,7 @@ end
 
 function Stochastic(d::Integer, f::Function,
                     monitor::Union{Bool, Vector{Int}}=true)
-  value = Array{Float64}(fill(0, d)...)
+  value = Array{Real}(fill(0, d)...)
   fx, src = modelfxsrc(depfxargs, f)
   s = ArrayStochastic(value, :nothing, Int[], fx, src, Symbol[],
                       NullUnivariateDistribution())
@@ -155,7 +155,7 @@ end
 #################### Updating ####################
 
 function setinits!(s::ScalarStochastic, m::Model, x::Real)
-  s.value = convert(Float64, x)
+  s.value = convert(Real, x)
   s.distr = s.eval(m)
   setmonitor!(s, s.monitor)
 end
diff --git a/src/model/mcmc.jl b/src/model/mcmc.jl
index 7dfd07f..3fccecb 100644
--- a/src/model/mcmc.jl
+++ b/src/model/mcmc.jl
@@ -45,6 +45,11 @@ function mcmc_master!(m::Model, window::UnitRange{Int}, burnin::Integer,
     "MCMC Simulation of $N Iterations x $K Chain" * "s"^(K > 1), verbose
   )
 
+  #lsts = Array(Array{Any,1}, 6)
+  #for k in chains
+  #  lsts[k] = [deepcopy(m), states[k],  window, burnin, thin, ChainProgress(frame, k, N)]
+  #end
+
   lsts = [
     Any[m, states[k], window, burnin, thin, ChainProgress(frame, k, N)]
     for k in chains
diff --git a/src/model/simulation.jl b/src/model/simulation.jl
index 6ac1e83..d4a8e0f 100644
--- a/src/model/simulation.jl
+++ b/src/model/simulation.jl
@@ -47,7 +47,8 @@ end
 function gradlogpdf!{T<:Real}(m::Model, x::AbstractVector{T}, block::Integer=0,
                               transform::Bool=false; dtype::Symbol=:forward)
   f = x -> logpdf!(m, x, block, transform)
-  gradient(f, convert(Vector{T}, x), dtype)
+  #gradient(f, convert(Vector{T}, x), dtype)
+  ForwardDiff.gradient(f, x)
 end
 
 
diff --git a/src/output/chains.jl b/src/output/chains.jl
index 164c969..e97fdc8 100644
--- a/src/output/chains.jl
+++ b/src/output/chains.jl
@@ -5,7 +5,7 @@
 function Chains{T<:AbstractString}(iters::Integer, params::Integer;
                start::Integer=1, thin::Integer=1, chains::Integer=1,
                names::Vector{T}=AbstractString[])
-  value = Array{Float64}(length(start:thin:iters), params, chains)
+  value = Array{Real}(length(start:thin:iters), params, chains)
   fill!(value, NaN)
   Chains(value, start=start, thin=thin, names=names)
 end
@@ -27,7 +27,7 @@ function Chains{T<:Real, U<:AbstractString, V<:Integer}(value::Array{T, 3};
     throw(DimensionMismatch("size(value, 3) and chains length differ"))
   end
 
-  v = convert(Array{Float64, 3}, value)
+  v = convert(Array{Real, 3}, value)
   Chains(v, range(start, thin, n), AbstractString[names...], Int[chains...])
 end
 
@@ -197,7 +197,7 @@ Base.last(c::AbstractChains) = last(c.range)
 
 function combine(c::AbstractChains)
   n, p, m = size(c.value)
-  value = Array{Float64}(n * m, p)
+  value = Array{Real}(n * m, p)
   for j in 1:p
     idx = 1
     for i in 1:n, k in 1:m
diff --git a/src/output/chainsummary.jl b/src/output/chainsummary.jl
index 711287f..91a4607 100644
--- a/src/output/chainsummary.jl
+++ b/src/output/chainsummary.jl
@@ -3,12 +3,12 @@
 #################### Types and Constructors ####################
 
 immutable ChainSummary
-  value::Array{Float64, 3}
+  value::Array{Real, 3}
   rownames::Vector{AbstractString}
   colnames::Vector{AbstractString}
   header::AbstractString
 
-  function ChainSummary(value::Array{Float64, 3},
+  function ChainSummary(value::Array{Real, 3},
                         rownames::Vector{AbstractString},
                         colnames::Vector{AbstractString},
                         header::AbstractString)
@@ -22,14 +22,14 @@ immutable ChainSummary
 end
 
 function ChainSummary{T<:AbstractString, U<:AbstractString}(
-                     value::Array{Float64, 3}, rownames::Vector{T},
+                     value::Array{Real, 3}, rownames::Vector{T},
                      colnames::Vector{U}, header::AbstractString)
   ChainSummary(copy(value), AbstractString[rownames...],
                AbstractString[colnames...], header)
 end
 
 function ChainSummary{T<:AbstractString, U<:AbstractString}(
-                     value::Matrix{Float64}, rownames::Vector{T},
+                     value::Matrix{Real}, rownames::Vector{T},
                      colnames::Vector{U}, header::AbstractString)
   dim = size(value)
   ChainSummary(reshape(value, dim[1], dim[2], 1), AbstractString[rownames...],
diff --git a/src/samplers/nuts.jl b/src/samplers/nuts.jl
index 8853727..ed79e6d 100644
--- a/src/samplers/nuts.jl
+++ b/src/samplers/nuts.jl
@@ -26,10 +26,10 @@ type NUTSTune <: SamplerTune
   end
 end
 
-NUTSTune(x::Vector{Float64}, logfgrad::Function, ::NullFunction; args...) =
+NUTSTune(x::Vector{Real}, logfgrad::Function, ::NullFunction; args...) =
   NUTSTune(x, nutsepsilon(x, logfgrad); args...)
 
-NUTSTune(x::Vector{Float64}, logfgrad::Function; args...) =
+NUTSTune(x::Vector{Real}, logfgrad::Function; args...) =
   NUTSTune(x, nutsepsilon(x, logfgrad), logfgrad; args...)
 
 NUTSTune(x::Vector, epsilon::Real; args...) =
@@ -126,7 +126,7 @@ function nuts_sub!(v::NUTSVariate, epsilon::Real, logfgrad::Function)
 end
 
 
-function leapfrog(x::Vector{Float64}, r::Vector{Float64}, grad::Vector{Float64},
+function leapfrog(x::Vector{Real}, r::Vector{Real}, grad::Vector{Real},
                   epsilon::Real, logfgrad::Function)
   r += (0.5 * epsilon) * grad
   x += epsilon * r
@@ -136,8 +136,8 @@ function leapfrog(x::Vector{Float64}, r::Vector{Float64}, grad::Vector{Float64},
 end
 
 
-function buildtree(x::Vector{Float64}, r::Vector{Float64},
-                   grad::Vector{Float64}, pm::Integer, j::Integer,
+function buildtree(x::Vector{Real}, r::Vector{Real},
+                   grad::Vector{Real}, pm::Integer, j::Integer,
                    epsilon::Real, logfgrad::Function, logp0::Real, logu0::Real)
   if j == 0
     xprime, rprime, logfprime, gradprime = leapfrog(x, r, grad, pm * epsilon,
@@ -180,8 +180,8 @@ function buildtree(x::Vector{Float64}, r::Vector{Float64},
 end
 
 
-function nouturn(xminus::Vector{Float64}, xplus::Vector{Float64},
-                 rminus::Vector{Float64}, rplus::Vector{Float64})
+function nouturn(xminus::Vector{Real}, xplus::Vector{Real},
+                 rminus::Vector{Real}, rplus::Vector{Real})
   xdiff = xplus - xminus
   dot(xdiff, rminus) >= 0 && dot(xdiff, rplus) >= 0
 end
@@ -189,7 +189,7 @@ end
 
 #################### Auxilliary Functions ####################
 
-function nutsepsilon(x::Vector{Float64}, logfgrad::Function)
+function nutsepsilon(x::Vector{Real}, logfgrad::Function)
   n = length(x)
   _, r0, logf0, grad0 = leapfrog(x, randn(n), zeros(n), 0.0, logfgrad)
   epsilon = 1.0
diff --git a/src/variate.jl b/src/variate.jl
index 8875ad0..a930041 100644
--- a/src/variate.jl
+++ b/src/variate.jl
@@ -6,6 +6,8 @@ Base.convert(::Type{Bool}, v::ScalarVariate) = convert(Bool, v.value)
 Base.convert{T<:Integer}(::Type{T}, v::ScalarVariate) = convert(T, v.value)
 Base.convert{T<:AbstractFloat}(::Type{T}, v::ScalarVariate) =
   convert(T, v.value)
+Base.convert{T<:Real}(::Type{T}, v::ScalarVariate) =
+  convert(T, v.value)
 
 Base.convert(::Type{Matrix}, v::MatrixVariate) = v.value
 Base.convert(::Type{Vector}, v::VectorVariate) = v.value
@@ -17,7 +19,7 @@ Base.unsafe_convert(::Type{Ptr{Float64}}, v::ArrayVariate) = pointer(v.value)
 
 macro promote_scalarvariate(V)
   quote
-    Base.promote_rule{T<:Real}(::Type{$(esc(V))}, ::Type{T}) = Float64
+    Base.promote_rule{T<:Real}(::Type{$(esc(V))}, ::Type{T}) = Real#Float64
   end
 end
 
@@ -34,7 +36,7 @@ Base.stride(v::ArrayVariate, k::Int) = stride(v.value, k)
 Base.getindex(v::ScalarVariate, ind::Int) = v.value[ind]
 
 Base.getindex(v::ScalarVariate, inds::Union{Range{Int}, Vector{Int}}) =
-  Float64[v[i] for i in inds]
+  Real[v[i] for i in inds]
 
 Base.getindex(v::ArrayVariate, inds::Int...) = getindex(v.value, inds...)