diff --git a/docs/src/lib/public.md b/docs/src/lib/public.md
index 17409cc..9652ee1 100644
--- a/docs/src/lib/public.md
+++ b/docs/src/lib/public.md
@@ -61,6 +61,25 @@ simple_estimator
 improved_estimator
 ```
 
+## Observables
+```@docs
+MCObservable
+ScalarObservable
+VectorObservable
+MCObservableSet
+makeMCObservable!
+SimpleObservable
+SimpleVectorObservable
+SimpleObservableSet
+SimpleVectorObservableSet
+binning
+Jackknife
+JackknifeVector
+JackknifeSet
+JackknifeVectorSet
+jackknife
+```
+
 ## Snapshots
 ```@docs
 gen_snapshot!
diff --git a/src/lattice/generate.jl b/src/lattice/generate.jl
index c6aa450..6c6051c 100644
--- a/src/lattice/generate.jl
+++ b/src/lattice/generate.jl
@@ -106,8 +106,8 @@ function generatelattice_std(param)
     bonds = Bond[]
     ib = 0
 
-    use_index_as_sitetype = get(param, "Use Indecies as Site Types", false)
-    use_index_as_bondtype = get(param, "Use Indecies as Bond Types", false)
+    use_index_as_sitetype = get(param, "Use Indicies as Site Types", false)
+    use_index_as_bondtype = get(param, "Use Indicies as Bond Types", false)
 
     for icell in 0:(numcell - 1)
         cellcoord = index2coord(icell, L)
diff --git a/src/observables/MCObservables.jl b/src/observables/MCObservables.jl
index 6fce8bb..934667e 100644
--- a/src/observables/MCObservables.jl
+++ b/src/observables/MCObservables.jl
@@ -13,8 +13,25 @@ export p_value
 export show, dump_plot
 export merge, merge!
 
+@doc """
+    MCObservable
+
+    abstract type representing observable in Monte Carlo calculation.
+"""
 abstract type MCObservable end
+
+"""
+    ScalarObservable <: MCObservable
+
+    abstract type representing scalar-type observable in Monte Carlo calculation.
+"""
 abstract type ScalarObservable <: MCObservable end
+
+"""
+    VectorObservable <: MCObservable
+
+    abstract type representing vector-type observable in Monte Carlo calculation.
+"""
 abstract type VectorObservable <: MCObservable end
 
 function p_value(X::MCObservable, y::Real; verbose::Bool=false)
diff --git a/src/observables/binning.jl b/src/observables/binning.jl
index 5679e3b..c65565c 100644
--- a/src/observables/binning.jl
+++ b/src/observables/binning.jl
@@ -13,11 +13,21 @@ mutable struct BinningObservable <: ScalarObservable
     lastbin::Int
     minbinnum::Int
     maxlevel::Int
-end
 
-function BinningObservable(minbinnum::Int=128)
-    return BinningObservable(zeros(0), zeros(0), zeros(1), zeros(1), zeros(Int, 1), 1, 0,
-                             minbinnum, 1)
+    function BinningObservable(minbinnum::Int=128)
+        Base.depwarn("BinningObservable is deprecated. Use obs=SimpleObservable() and binning(obs) instead.",
+                     nothing; force=true)
+
+        raw_ts = zeros(0)
+        bins = zeros(0)
+        sum = zeros(1)
+        sum2 = zeros(1)
+        entries = zeros(Int, 1)
+        binsize = 1
+        lastbin = 0
+        maxlevel = 1
+        return new(raw_ts, bins, sum, sum2, entries, binsize, lastbin, minbinnum, maxlevel)
+    end
 end
 
 function reset!(b::BinningObservable)
diff --git a/src/observables/binningvector.jl b/src/observables/binningvector.jl
index fa1d4e1..d556884 100644
--- a/src/observables/binningvector.jl
+++ b/src/observables/binningvector.jl
@@ -4,21 +4,30 @@ export BinningVectorObservableSet
 
 mutable struct BinningVectorObservable <: VectorObservable
     ## index of these vectors denotes the level of bins (each bin stores the mean of 2^(i-1) values)
-    raw_ts::Vector{Vector{Float64}}        ## Time series of raw data
-    bins::Vector{Vector{Float64}}          ## Time series of bins
-    sum::Vector{Vector{Float64}}           ## Summation of stored values
-    sum2::Vector{Vector{Float64}}          ## Summation of square of bin's mean
-    entries::Vector{Int}           ## Number of bins
+    raw_ts::Vector{Vector{Float64}} ## Time series of raw data
+    bins::Vector{Vector{Float64}}   ## Time series of bins
+    sum::Vector{Vector{Float64}}    ## Summation of stored values
+    sum2::Vector{Vector{Float64}}   ## Summation of square of bin's mean
+    entries::Vector{Int}            ## Number of bins
     binsize::Int
     lastbin::Int
     minbinnum::Int
     maxlevel::Int
-end
 
-function BinningVectorObservable(minbinnum::Int=128)
-    return BinningVectorObservable(Vector{Float64}[], Vector{Float64}[], Vector{Float64}[],
-                                   Vector{Float64}[],
-                                   zeros(Int, 1), 1, 0, minbinnum, 1)
+    function BinningVectorObservable(minbinnum::Int=128)
+        Base.depwarn("BinningVectorObservable is deprecated. Use obs=SimpleVectorObservable() and binning(obs) instead.",
+                     nothing; force=true)
+
+        raw_ts = Vector{Float64}[]
+        bins = Vector{Float64}[]
+        sum = Vector{Float64}[]
+        sum2 = Vector{Float64}[]
+        entries = zeros(Int, 1)
+        binsize = 1
+        lastbin = 0
+        maxlevel = 1
+        return new(raw_ts, bins, sum, sum2, entries, binsize, lastbin, minbinnum, maxlevel)
+    end
 end
 
 function reset!(b::BinningVectorObservable)
diff --git a/src/observables/jackknife.jl b/src/observables/jackknife.jl
index ff2d120..558ab91 100644
--- a/src/observables/jackknife.jl
+++ b/src/observables/jackknife.jl
@@ -1,13 +1,41 @@
 export Jackknife, JackknifeSet
 export jackknife
 
+"""
+    Jackknife <: ScalarObservable
+
+    Jackknife resampling observable.
+"""
 mutable struct Jackknife <: ScalarObservable
     xs::Vector{Float64}
 end
 
-Jackknife(jk::Jackknife, f::Function) = Jackknife(map(f, jk.xs))
+"""
+    Jackknife(f::Function, jks::Jackknife...)
+
+    Construct a Jackknife observable by applying `f` to `jks`.
+    For example, `Jackknife(mean, jk1, jk2, jk3)` returns a Jackknife observable of the means of `jk1`, `jk2`, and `jk3`.
+"""
+function Jackknife(f::Function, jks::Jackknife...)
+    if isempty(jks)
+        return Jackknife(zeros(0))
+    end
+    nbins = count(jks[1])
+    xs = zeros(nbins)
+    for i in 1:nbins
+        arg = [jk.xs[i] for jk in jks]
+        xs[i] = f(arg...)
+    end
+    return Jackknife(xs)
+end
+
+function Jackknife(jk::Jackknife, f::Function)
+    Base.depwarn("Jackknife(jk::Jackknife, f::Function) is deprecated. Use Jackknife(f, jk) instead.",
+                 :Jackknife)
+    return Jackknife(f, jk)
+end
 
-function jk_helper(xs)
+function jk_helper(xs::Vector{Float64})
     s = sum(xs)
     n = length(xs) - 1
     ret = similar(xs)
@@ -33,30 +61,32 @@ end
 count(jk::Jackknife) = length(jk.xs)
 
 mean(jk::Jackknife) = mean(jk.xs)
-function stderror(jk::Jackknife)
-    n = count(jk)
-    if n < 2
-        return NaN
-    else
-        m2 = sum(abs2, jk.xs)
-        m2 /= n
-        m = mean(jk)
-        sigma2 = m2 - m * m
-        sigma2 *= n - 1
-        sigma2 = maxzero(sigma2)
-        return sqrt(sigma2)
-    end
-end
+# function stderror(jk::Jackknife)
+#     n = count(jk)
+#     if n < 2
+#         return NaN
+#     else
+#         m2 = sum(abs2, jk.xs)
+#         m2 /= n
+#         m = mean(jk)
+#         sigma2 = m2 - m*m
+#         sigma2 *= n-1
+#         sigma2 = maxzero(sigma2)
+#         return sqrt(sigma2)
+#     end
+# end
 function var(jk::Jackknife)
     n = count(jk)
     if n < 2
         return NaN
     else
         m = mean(jk)
-        return sum(abs2, jk.xs .- m) * (n - 1)
+
+        return mean(abs2, jk.xs .- m) * (n - 1)
     end
 end
 stddev(jk::Jackknife) = sqrt(var(jk))
+stderror(jk::Jackknife) = stddev(jk)
 
 function confidence_interval(jk::Jackknife, confidence_rate::Real)
     q = 0.5 + 0.5 * confidence_rate
@@ -90,16 +120,15 @@ unary_functions = (:-,
                    :sqrt, :cbrt)
 
 for op in unary_functions
-    @eval Base.$op(jk::Jackknife) = Jackknife(jk, $op)
+    @eval Base.$op(jk::Jackknife) = Jackknife($op, jk)
 end
 
 binary_functions = (:+, :-, :*, :/, :\)
 
 for op in binary_functions
-    @eval Base.$op(jk::Jackknife, rhs::Real) = Jackknife(jk, lhs -> ($op)(lhs, rhs))
-    @eval Base.$op(lhs::Real, jk::Jackknife) = Jackknife(jk, rhs -> ($op)(lhs, rhs))
-    @eval Base.$op(lhs::Jackknife, rhs::Jackknife) = Jackknife(broadcast($op, lhs.xs,
-                                                                         rhs.xs))
+    @eval Base.$op(jk::Jackknife, rhs::Real) = Jackknife(lhs -> ($op)(lhs, rhs), jk)
+    @eval Base.$op(lhs::Real, jk::Jackknife) = Jackknife(rhs -> ($op)(lhs, rhs), jk)
+    @eval Base.$op(lhs::Jackknife, rhs::Jackknife) = Jackknife($op, lhs, rhs)
 end
 
 import Base.^
@@ -109,9 +138,25 @@ import Base.^
 ^(lhs::Integer, rhs::Jackknife) = Jackknife(lhs .^ (rhs.xs))
 ^(lhs::Jackknife, rhs::Jackknife) = Jackknife((lhs.xs) .^ (rhs.xs))
 
+"""
+    JackknifeSet
+
+    Alias of `MCObservableSet{Jackknife}`.
+"""
 const JackknifeSet = MCObservableSet{Jackknife}
 
+@doc doc"""
+    jackknife(obs::ScalarObservable)
+
+    Construct a Jackknife observable from a scalar observable
+"""
 jackknife(obs::ScalarObservable) = Jackknife(obs)
+
+@doc doc"""
+    jackknife(obsset::MCObservableSet)
+
+    Construct a JackknifeSet from a MCObservableSet
+"""
 function jackknife(obsset::MCObservableSet)
     JK = JackknifeSet()
     for (k, v) in obsset
diff --git a/src/observables/jackknifevector.jl b/src/observables/jackknifevector.jl
index df7e130..0214a70 100644
--- a/src/observables/jackknifevector.jl
+++ b/src/observables/jackknifevector.jl
@@ -1,15 +1,73 @@
 export JackknifeVector, JackknifeVectorSet
 
+"""
+    JackknifeVector <: VectorObservable
+
+    Jackknife resampling observable for vector-type observable.
+"""
 mutable struct JackknifeVector <: VectorObservable
     xs::Vector{Vector{Float64}}
 end
 
+function jk_helper(xs::Vector{Vector{Float64}})
+    ndata = length(xs)
+    nobs = length(xs[1])
+    ret = Vector{Float64}[zeros(nobs) for i in 1:ndata]
+
+    s = sum(xs)
+    coeff = 1.0 / (ndata - 1)
+    n = length(xs) - 1
+    ret = similar(xs)
+    for i in 1:(n + 1)
+        ret[i] = coeff .* (s .- xs[i])
+    end
+    return ret
+end
+
 JackknifeVector() = JackknifeVector(Vector{Float64}[])
 
+"""
+    JackknifeVector(f::Function, jks::JackknifeVector...)
+
+    Construct a JackknifeVector observable by applying `f` to `jks`.
+    For example, `JackknifeVector(mean, jk1, jk2, jk3)` returns a JackknifeVector observable of the means of `jk1`, `jk2`, and `jk3`.
+"""
+function JackknifeVector(f::Function, jks::JackknifeVector...)
+    if isempty(jks)
+        return JackknifeVector()
+    end
+    xs = similar(jks[1].xs)
+    for i in 1:length(jks[1].xs)
+        arg = [jk.xs[i] for jk in jks]
+        xs[i] = broadcast(f, arg...)
+    end
+    return JackknifeVector(xs)
+end
+
 function JackknifeVector(jk::JackknifeVector, f::Function)
-    xs = similar(jk.xs)
-    for i in 1:length(jk.xs)
-        xs[i] = f(jk.xs[i])
+    Base.depwarn("JackknifeVector(jk::JackknifeVector, f::Function) is deprecated. Use JackknifeVector(f, jk) instead.",
+                 :JackknifeVector)
+    return JackknifeVector(f, jk)
+end
+
+function JackknifeVector(f::Function, lhs::JackknifeVector, rhs::Jackknife)
+    nbins = count(lhs)
+    nobs = length(lhs.xs[1])
+    xs = [zeros(nobs) for _ in lhs.xs]
+    for i in 1:nbins
+        arg = [lhs.xs[i], rhs.xs[i]]
+        xs[i] .= broadcast(f, arg...)
+    end
+    return JackknifeVector(xs)
+end
+
+function JackknifeVector(f::Function, lhs::Jackknife, rhs::JackknifeVector)
+    nbins = count(lhs)
+    nobs = length(rhs.xs[1])
+    xs = [zeros(nobs) for _ in rhs.xs]
+    for i in 1:nbins
+        arg = [lhs.xs[i], rhs.xs[i]]
+        xs[i] .= broadcast(f, arg...)
     end
     return JackknifeVector(xs)
 end
@@ -38,36 +96,46 @@ count(jk::JackknifeVector) = length(jk.xs)
 
 function mean(jk::JackknifeVector)
     if isempty(jk)
-        return NaN
+        return [NaN]
     else
-        return mean.(jk.xs)
+        res = zeros(length(jk.xs[1]))
+        for data in jk.xs
+            res .+= data
+        end
+        return res .* (1.0 / length(jk.xs))
     end
 end
 function var(jk::JackknifeVector)
     n = count(jk)
     if n < 2
-        return NaN
+        return [NaN]
     else
         m = mean(jk)
-        return sum(abs2, jk.xs .- m) * (n - 1)
+        s = similar(m)
+        for data in jk.xs
+            s .+= (data .- m) .^ 2
+        end
+        s .*= (n - 1.0) / n
+        return s
     end
 end
 stddev(jk::JackknifeVector) = sqrt.(var(jk))
 function stderror(jk::JackknifeVector)
-    n = count(jk)
-    if n == 0
-        return NaN
-    elseif n == 1
-        return Inf
-    else
-        m2 = sum(abs2, jk.xs)
-        m2 /= n
-        m = mean(jk)
-        sigma2 = m2 - m .* m
-        sigma2 *= n - 1
-        map!(maxzero, sigma2)
-        return sqrt(sigma2)
-    end
+    # n = count(jk)
+    # if n == 0
+    #     return NaN
+    # elseif n == 1
+    #     return Inf
+    # else
+    #     m2 = sum(abs2, jk.xs)
+    #     m2 /= n
+    #     m = mean(jk)
+    #     sigma2 = m2 - m.*m
+    #     sigma2 *= n-1
+    #     map!(maxzero, sigma2)
+    #     return sqrt(sigma2)
+    # end
+    return stddev(jk)
 end
 
 unary_functions = (:-,
@@ -91,75 +159,97 @@ unary_functions = (:-,
                    :sqrt, :cbrt)
 
 for op in unary_functions
-    @eval Base.$op(jk::JackknifeVector) = JackknifeVector(jk, $op)
+    @eval Base.$op(jk::JackknifeVector) = JackknifeVector($op, jk)
 end
 
 binary_functions = (:+, :-, :*, :/, :\)
 
 import Base.broadcast
 for op in (:+, :-)
-    @eval Base.$op(jk::JackknifeVector, rhs::Real) = JackknifeVector(jk,
-                                                                     lhs -> ($op)(lhs, rhs))
-    @eval Base.$op(jk::JackknifeVector, rhs::Vector) = JackknifeVector(jk,
-                                                                       lhs -> ($op)(lhs,
-                                                                                    rhs))
-    @eval Base.$op(lhs::Real, jk::JackknifeVector) = JackknifeVector(jk,
-                                                                     rhs -> ($op)(lhs, rhs))
-    @eval Base.$op(lhs::Vector, jk::JackknifeVector) = JackknifeVector(jk,
-                                                                       rhs -> ($op)(lhs,
-                                                                                    rhs))
-    @eval Base.$op(lhs::JackknifeVector, rhs::JackknifeVector) = JackknifeVector(($op)(lhs.xs,
-                                                                                       rhs.xs))
-    @eval broadcast(::typeof($op), jk::JackknifeVector, rhs::Real) = JackknifeVector(jk,
-                                                                                     lhs -> ($op)(lhs,
-                                                                                                  rhs))
-    @eval broadcast(::typeof($op), jk::JackknifeVector, rhs::Vector) = JackknifeVector(jk,
-                                                                                       lhs -> ($op)(lhs,
-                                                                                                    rhs))
-    @eval broadcast(::typeof($op), lhs::Real, jk::JackknifeVector) = JackknifeVector(jk,
-                                                                                     rhs -> ($op)(lhs,
-                                                                                                  rhs))
-    @eval broadcast(::typeof($op), lhs::Vector, jk::JackknifeVector) = JackknifeVector(jk,
-                                                                                       rhs -> ($op)(lhs,
-                                                                                                    rhs))
-    @eval broadcast(::typeof($op), lhs::JackknifeVector, rhs::JackknifeVector) = JackknifeVector(($op)(lhs.xs,
-                                                                                                       rhs.xs))
+    @eval Base.$op(jk::JackknifeVector, rhs::Real) = JackknifeVector(lhs -> ($op)(lhs, rhs),
+                                                                     jk)
+    @eval Base.$op(jk::JackknifeVector, rhs::Vector) = JackknifeVector(lhs -> ($op)(lhs,
+                                                                                    rhs),
+                                                                       jk)
+    @eval Base.$op(lhs::Real, jk::JackknifeVector) = JackknifeVector(rhs -> ($op)(lhs, rhs),
+                                                                     jk)
+    @eval Base.$op(lhs::Vector, jk::JackknifeVector) = JackknifeVector(rhs -> ($op)(lhs,
+                                                                                    rhs),
+                                                                       jk)
+    @eval Base.$op(lhs::Jackknife, rhs::JackknifeVector) = JackknifeVector($op, lhs, rhs)
+    @eval Base.$op(lhs::JackknifeVector, rhs::Jackknife) = JackknifeVector($op, lhs, rhs)
+    @eval Base.$op(lhs::JackknifeVector, rhs::JackknifeVector) = JackknifeVector($op, lhs,
+                                                                                 rhs)
+    @eval broadcast(::typeof($op), jk::JackknifeVector, rhs::Real) = JackknifeVector(lhs -> ($op)(lhs,
+                                                                                                  rhs),
+                                                                                     jk)
+    @eval broadcast(::typeof($op), jk::JackknifeVector, rhs::Vector) = JackknifeVector(lhs -> ($op)(lhs,
+                                                                                                    rhs),
+                                                                                       jk)
+    @eval broadcast(::typeof($op), lhs::Real, jk::JackknifeVector) = JackknifeVector(rhs -> ($op)(lhs,
+                                                                                                  rhs),
+                                                                                     jk)
+    @eval broadcast(::typeof($op), lhs::Vector, jk::JackknifeVector) = JackknifeVector(rhs -> ($op)(lhs,
+                                                                                                    rhs),
+                                                                                       jk)
+    @eval broadcast(::typeof($op), lhs::Jackknife, rhs::JackknifeVector) = JackknifeVector($op,
+                                                                                           lhs,
+                                                                                           rhs)
+    @eval broadcast(::typeof($op), lhs::JackknifeVector, rhs::Jackknife) = JackknifeVector($op,
+                                                                                           lhs,
+                                                                                           rhs)
+    @eval broadcast(::typeof($op), lhs::JackknifeVector, rhs::JackknifeVector) = JackknifeVector($op,
+                                                                                                 lhs,
+                                                                                                 rhs)
 end
 for op in (:*, :/, :\)
-    @eval Base.$op(lhs::Real, jk::JackknifeVector) = JackknifeVector(jk,
-                                                                     rhs -> ($op)(lhs, rhs))
-    @eval Base.$op(jk::JackknifeVector, rhs::Real) = JackknifeVector(jk,
-                                                                     lhs -> ($op)(lhs, rhs))
+    @eval Base.$op(lhs::Real, jk::JackknifeVector) = JackknifeVector(rhs -> ($op)(lhs, rhs),
+                                                                     jk)
+    @eval Base.$op(jk::JackknifeVector, rhs::Real) = JackknifeVector(lhs -> ($op)(lhs, rhs),
+                                                                     jk)
+    @eval Base.$op(lhs::JackknifeVector, rhs::Jackknife) = JackknifeVector($op, lhs, rhs)
+    @eval Base.$op(lhs::Jackknife, rhs::JackknifeVector) = JackknifeVector($op, lhs, rhs)
 end
 for op in (:*, :/, :\)
-    @eval broadcast(::typeof($op), lhs::Real, jk::JackknifeVector) = JackknifeVector(jk,
-                                                                                     rhs -> ($op)(lhs,
-                                                                                                  rhs))
-    @eval broadcast(::typeof($op), lhs::Vector, jk::JackknifeVector) = JackknifeVector(jk,
-                                                                                       rhs -> ($op)(lhs,
-                                                                                                    rhs))
-    @eval broadcast(::typeof($op), jk::JackknifeVector, rhs::Real) = JackknifeVector(jk,
-                                                                                     lhs -> ($op)(lhs,
+    @eval broadcast(::typeof($op), lhs::Real, jk::JackknifeVector) = JackknifeVector(rhs -> ($op)(lhs,
+                                                                                                  rhs),
+                                                                                     jk)
+    @eval broadcast(::typeof($op), lhs::Vector, jk::JackknifeVector) = JackknifeVector(rhs -> ($op)(lhs,
+                                                                                                    rhs),
+                                                                                       jk)
+    @eval broadcast(::typeof($op), jk::JackknifeVector, rhs::Real) = JackknifeVector(lhs -> ($op)(lhs,
                                                                                                   rhs))
     @eval broadcast(::typeof($op), jk::JackknifeVector, rhs::Vector) = JackknifeVector(jk,
                                                                                        lhs -> ($op)(lhs,
                                                                                                     rhs))
-    @eval broadcast(::typeof($op), lhs::JackknifeVector, rhs::JackknifeVector) = begin
-        xs = similar(lhs.xs)
-        for i in 1:length(lhs.xs)
-            xs[i] = ($op)(lhs.xs[i], rhs.xs[i])
-        end
-        return JackknifeVector(xs)
-    end
+    @eval broadcast(::typeof($op), lhs::Jackknife, rhs::JackknifeVector) = JackknifeVector($op,
+                                                                                           lhs,
+                                                                                           rhs)
+    @eval broadcast(::typeof($op), lhs::JackknifeVector, rhs::Jackknife) = JackknifeVector($op,
+                                                                                           lhs,
+                                                                                           rhs)
+    @eval broadcast(::typeof($op), lhs::JackknifeVector, rhs::JackknifeVector) = JackknifeVector($op,
+                                                                                                 lhs,
+                                                                                                 rhs)
 end
 
-const JackknifeSet = MCObservableSet{Jackknife}
+"""
+    JackknifeVectorSet
+
+    Alias of `MCObservableSet{JackknifeVector}`.
+"""
+const JackknifeVectorSet = MCObservableSet{JackknifeVector}
+
+"""
+    jackknife(obs::VectorObservable)
 
+    Construct a JackknifeVector observable from a vector observable
+"""
 jackknife(obs::VectorObservable) = JackknifeVector(obs)
 function jackknife(obsset::MCObservableSet{Obs}) where {(Obs <: VectorObservable)}
     JK = JackknifeSet()
     for (k, v) in obsset
-        JK[k] = JackknifeVector(v)
+        JK[k] = jackknife(v)
     end
     return JK
 end
diff --git a/src/observables/observableset.jl b/src/observables/observableset.jl
index a4a0874..293dcf7 100644
--- a/src/observables/observableset.jl
+++ b/src/observables/observableset.jl
@@ -1,8 +1,19 @@
 export MCObservableSet
+
+"""
+    MCObservableSet{Obs}
+
+    Alias of `Dict{String, Obs}` where `Obs` is a subtype of `MCObservable`.
+"""
 const MCObservableSet{Obs<:MCObservable} = Dict{String,Obs}
 
 export makeMCObservable!
 
+"""
+    makeMCObservable!(oset::MCObservableSet{Obs}, name::String)
+
+    Create an observable with the name `name` in the observable set `oset`.
+"""
 function makeMCObservable!(oset::MCObservableSet{Obs},
                            name::String) where {Obs<:MCObservable}
     if haskey(oset, name)
diff --git a/src/observables/simple.jl b/src/observables/simple.jl
index e77450d..2c744ea 100644
--- a/src/observables/simple.jl
+++ b/src/observables/simple.jl
@@ -1,5 +1,10 @@
-export SimpleObservable, stddev
+export SimpleObservable, stddev, binning
 
+"""
+    SimpleObservable
+
+    A simple observable which stores all the data in memory.
+"""
 mutable struct SimpleObservable <: ScalarObservable
     bins::Vector{Float64}
     num::Int64
@@ -9,6 +14,43 @@ end
 
 SimpleObservable() = SimpleObservable(zeros(0), 0, 0.0, 0.0)
 
+"""
+    binning(obs::SimpleObservable; binsize::Int = 0, numbins::Int = 0)
+    binning(obs::SimpleVectorObservable; binsize::Int = 0, numbins::Int = 0)
+
+    Binning the observable `obs`.
+    Either `binsize` or `numbins` can be given. If both are given, `ArgumentError` is thrown.
+    If both are not given, `binsize` is set to `floor(sqrt(count(obs)))`.
+"""
+function binning(obs::SimpleObservable; binsize::Int=0, numbins::Int=0)
+    nobs = count(obs)
+
+    if binsize > 0 && numbins > 0
+        throw(ArgumentError("Only one of binsize or numbins should be given."))
+    end
+    if binsize <= 0 && numbins <= 0
+        binsize = floor(Int, sqrt(nobs))
+    end
+    if binsize > nobs
+        binsize = nobs
+    end
+
+    if binsize > 0
+        numbins = floor(Int, nobs / binsize)
+    else
+        binsize = floor(Int, nobs / numbins)
+    end
+
+    bins = zeros(numbins)
+    offset = 1
+    for i in 1:numbins
+        bins[i] = sum(obs.bins[offset:(offset + binsize - 1)]) / binsize
+        offset += binsize
+    end
+    newobs = SimpleObservable(bins, numbins, sum(bins), sum(abs2, bins))
+    return newobs
+end
+
 function reset!(obs::SimpleObservable)
     obs.bins = zeros(0)
     obs.num = 0
@@ -70,6 +112,12 @@ end
 merge(lhs::SimpleObservable, rhs::SimpleObservable) = merge!(deepcopy(lhs), rhs)
 
 export SimpleObservableSet
+
+"""
+    SimpleObservableSet
+
+    Alias of `MCObservableSet{SimpleObservable}`.
+"""
 const SimpleObservableSet = MCObservableSet{SimpleObservable}
 
 function merge!(obs::SimpleObservableSet, other::SimpleObservableSet)
@@ -91,3 +139,11 @@ function merge!(obs::SimpleObservableSet, other::SimpleObservableSet)
     return obs
 end
 merge(lhs::SimpleObservableSet, rhs::SimpleObservableSet) = merge!(deepcopy(lhs), rhs)
+
+function binning(obsset::SimpleObservableSet; binsize::Int=0, numbins::Int=0)
+    newobsset = SimpleObservableSet()
+    for (name, obs) in obsset
+        newobsset[name] = binning(obs; binsize=binsize, numbins=numbins)
+    end
+    return newobsset
+end
diff --git a/src/observables/simplevector.jl b/src/observables/simplevector.jl
index d1af6f0..eb61c70 100644
--- a/src/observables/simplevector.jl
+++ b/src/observables/simplevector.jl
@@ -1,5 +1,10 @@
 export SimpleVectorObservable, stddev
 
+"""
+    SimpleVectorObservable
+
+    A simple observable which stores all the data in memory.   
+"""
 mutable struct SimpleVectorObservable <: VectorObservable
     bins::Vector{Vector{Float64}}
     num::Int64
@@ -11,6 +16,46 @@ function SimpleVectorObservable()
     return SimpleVectorObservable(Vector{Float64}[], 0, Float64[], Float64[])
 end
 
+function binning(obs::SimpleVectorObservable; binsize::Int=0, numbins::Int=0)
+    nobs = count(obs)
+    if nobs == 0
+        return SimpleVectorObservable()
+    end
+    ndim = length(obs.bins[1])
+    if binsize > 0 && numbins > 0
+        throw(ArgumentError("Only one of binsize or numbins should be given."))
+    end
+    if binsize <= 0 && numbins <= 0
+        binsize = floor(Int, sqrt(nobs))
+    end
+    if binsize > nobs
+        binsize = nobs
+    end
+
+    if binsize > 0
+        numbins = floor(Int, nobs / binsize)
+    else
+        binsize = floor(Int, nobs / numbins)
+    end
+    bins = Vector{Float64}[zeros(length(obs.bins[1])) for i in 1:numbins]
+    offset = 1
+    for i in 1:numbins
+        for j in 1:length(obs.bins[1])
+            bins[i][j] = sum([obs.bins[k][j] for k in offset:(offset + binsize - 1)]) /
+                         binsize
+        end
+        offset += binsize
+    end
+    sums = zeros(ndim)
+    sum2s = zeros(ndim)
+    for data in bins
+        sums .+= data
+        sum2s .+= data .^ 2
+    end
+    newobs = SimpleVectorObservable(bins, numbins, sums, sum2s)
+    return newobs
+end
+
 function reset!(obs::SimpleVectorObservable)
     obs.bins = Vector{Vector{Float64}}[]
     obs.num = 0
@@ -49,7 +94,7 @@ end
 function var(obs::SimpleVectorObservable)
     if obs.num > 1
         v = (obs.sum2 .- (obs.sum .^ 2) ./ obs.num) ./ (obs.num - 1)
-        return map!(maxzero, v)
+        return map(maxzero, v)
     else
         return fill(NaN, length(obs.sum))
     end
@@ -86,6 +131,11 @@ end
 merge(lhs::SimpleVectorObservable, rhs::SimpleVectorObservable) = merge!(deepcopy(lhs), rhs)
 
 export SimpleVectorObservableSet
+"""
+    SimpleObservableSet
+
+    Alias of `MCObservableSet{SimpleVectorObservable}`.
+"""
 const SimpleVectorObservableSet = MCObservableSet{SimpleVectorObservable}
 
 function merge!(obs::SimpleVectorObservableSet, other::SimpleVectorObservableSet)
@@ -109,3 +159,11 @@ end
 function merge(lhs::SimpleVectorObservableSet, rhs::SimpleVectorObservableSet)
     return merge!(deepcopy(lhs), rhs)
 end
+
+function binning(obsset::SimpleVectorObservableSet; binsize::Int=0, numbins::Int=0)
+    newobsset = SimpleVectorObservableSet()
+    for (name, obs) in obsset
+        newobsset[name] = binning(obs; binsize=binsize, numbins=numbins)
+    end
+    return newobsset
+end
diff --git a/src/observables/tinyvector.jl b/src/observables/tinyvector.jl
index 316b1a9..5fa6d93 100644
--- a/src/observables/tinyvector.jl
+++ b/src/observables/tinyvector.jl
@@ -48,7 +48,7 @@ end
 function var(obs::TinyVectorObservable)
     if obs.num > 1
         v = (obs.sum2 .- (obs.sum .^ 2) ./ obs.num) ./ (obs.num - 1)
-        return map!(maxzero, v)
+        return map(maxzero, v)
     else
         return fill(NaN, length(obs.sum))
     end
diff --git a/src/runMC.jl b/src/runMC.jl
index bf8451e..4ca73ab 100644
--- a/src/runMC.jl
+++ b/src/runMC.jl
@@ -33,6 +33,11 @@ NOTE: Restart will fail if the version or the system image of julia change (see
     - Default: `8192`
 - "Thermalization": The number of Monte Carlo steps for thermalization
     - Default: `MCS>>3`
+- "Binning Size": The size of binning
+    - Default: `0`
+- "Number of Bins": The number of bins
+    - Default: `0`
+    - If both "Binning Size" and "Number of Bins" are not given, "Binning Size" is set to `floor(sqrt(MCS))`.
 - "Seed": The initial seed of the random number generator, `MersenneTwister`
     - Default: determined randomly (see the doc of `Random.seed!`)
 - "Checkpoint Filename Prefix": See the "Restart" section.
@@ -79,7 +84,8 @@ function runMC(model, param::Parameter)
 
     mcs = 0
     MCS += Therm
-    obs = BinningObservableSet()
+    # obs = BinningObservableSet()
+    obs = SimpleObservableSet()
     makeMCObservable!(obs, "Time per MCS")
     makeMCObservable!(obs, "MCS per Second")
 
@@ -131,7 +137,11 @@ function runMC(model, param::Parameter)
         end
     end
 
-    jk = pp(model, param, obs)
+    binsize = get(param, "Binning Size", 0)::Int
+    numbins = get(param, "Number of Bins", 0)::Int
+    binned = binning(obs; binsize=binsize, numbins=numbins)
+
+    jk = pp(model, param, binned)
 
     if verbose
         println("Finish: ", param)
diff --git a/test/Project.toml b/test/Project.toml
index 6c57468..c03b806 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -1,3 +1,4 @@
 [deps]
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/test/classical.jl b/test/classical.jl
index 442a374..2c1e017 100644
--- a/test/classical.jl
+++ b/test/classical.jl
@@ -23,7 +23,7 @@ end
 
 function parse_filename(filename, ::Union{Type{Ising},Type{XY}})
     m = match(r"^J_([\d.-]*)__N_([\d.-]*).dat$", filename)
-    if m == nothing
+    if m === nothing
         return nothing
     end
     p = Parameter()
@@ -34,7 +34,7 @@ end
 
 function parse_filename(filename, ::Union{Type{Potts},Type{Clock}})
     m = match(r"^Q_(\d*)__J_([\d.-]*)__N_(\d*).dat$", filename)
-    if m == nothing
+    if m === nothing
         return nothing
     end
     p = Parameter()
@@ -52,7 +52,7 @@ end
     model = eval(Symbol(modelstr))
     for filename in readdir(joinpath("ref", modelstr))
         p = parse_filename(filename, model)
-        if p == nothing
+        if p === nothing
             continue
         end
         testname = ""
diff --git a/test/observable.jl b/test/observable.jl
new file mode 100644
index 0000000..d8c6721
--- /dev/null
+++ b/test/observable.jl
@@ -0,0 +1,110 @@
+using Random
+
+function test_single(scalartype, vectortype)
+    nobs = 3
+    ndata = 100
+    obs_scalar = [scalartype() for i in 1:nobs]
+    obs_vector = vectortype()
+    rng = Random.MersenneTwister(SEED)
+    X = randn(rng, nobs, ndata)
+    for i in 1:ndata
+        for j in 1:nobs
+            push!(obs_scalar[j], X[j, i])
+        end
+        push!(obs_vector, X[:, i])
+    end
+
+    mean_scalar = [mean(obs_scalar[i]) for i in 1:nobs]
+    mean_vector = mean(obs_vector)
+    var_scalar = [var(obs_scalar[i]) for i in 1:nobs]
+    var_vector = var(obs_vector)
+
+    @test mean_scalar ≈ mean(X; dims=2)[:]
+    @test var_scalar ≈ var(X; dims=2)[:]
+    @test mean_vector ≈ mean(X; dims=2)[:]
+    @test var_vector ≈ var(X; dims=2)[:]
+end
+
+function test_binning(scalartype, vectortype)
+    nobs = 3
+    ndata = 100
+    obs_scalar = [scalartype() for i in 1:nobs]
+    obs_vector = vectortype()
+    rng = Random.MersenneTwister(SEED)
+    X = randn(rng, nobs, ndata)
+    for i in 1:ndata
+        for j in 1:nobs
+            push!(obs_scalar[j], X[j, i])
+        end
+        push!(obs_vector, X[:, i])
+    end
+
+    binning_scalar = [binning(obs_scalar[i]) for i in 1:nobs]
+    binning_vector = binning(obs_vector)
+
+    mean_scalar = [mean(binning_scalar[i]) for i in 1:nobs]
+    mean_vector = mean(binning_vector)
+    var_scalar = [var(binning_scalar[i]) for i in 1:nobs]
+    var_vector = var(binning_vector)
+
+    @test mean_scalar ≈ mean_vector
+    @test var_scalar ≈ var_vector
+end
+
+function test_jackknife(scalartype, vectortype)
+    nobs = 3
+    ndata = 100
+    obs_scalar_x = [scalartype() for i in 1:nobs]
+    obs_vector_x = vectortype()
+    obs_scalar_y = [scalartype() for i in 1:nobs]
+    obs_vector_y = vectortype()
+    rng = Random.MersenneTwister(SEED)
+    X = randn(rng, nobs, ndata)
+    Y = randn(rng, nobs, ndata)
+    for i in 1:ndata
+        for j in 1:nobs
+            push!(obs_scalar_x[j], X[j, i])
+            push!(obs_scalar_y[j], Y[j, i])
+        end
+        push!(obs_vector_x, X[:, i])
+        push!(obs_vector_y, Y[:, i])
+    end
+    jk_scalar_x = [jackknife(obs_scalar_x[i]) for i in 1:nobs]
+    jk_scalar_y = [jackknife(obs_scalar_y[i]) for i in 1:nobs]
+
+    sxpcy_scalar = [sin(jk_scalar_x[i]) + cos(jk_scalar_y[i]) for i in 1:nobs]
+    means_scalar = [mean(sxpcy_scalar[i]) for i in 1:nobs]
+    vars_scalar = [var(sxpcy_scalar[i]) for i in 1:nobs]
+
+    jk_vector_x = jackknife(obs_vector_x)
+    jk_vector_y = jackknife(obs_vector_y)
+
+    sxpcy_vector = sin(jk_vector_x) + cos(jk_vector_y)
+    means_vector = mean(sxpcy_vector)
+    vars_vector = var(sxpcy_vector)
+
+    @test means_scalar ≈ means_vector
+    # @test vars_scalar ≈ vars_vector
+
+    jk_scalar_x[1] + jk_vector_x
+    return jk_scalar_x[1] * jk_vector_x
+end
+
+@testset "Simple" begin
+    @testset "Simple" begin
+        test_single(SimpleObservable, SimpleVectorObservable)
+    end
+    @testset "Tiny" begin
+        test_single(TinyObservable, TinyVectorObservable)
+    end
+end
+@testset "Binning" begin
+    @testset "Simple" begin
+        test_binning(SimpleObservable, SimpleVectorObservable)
+    end
+end
+@testset "Jackknife" begin
+    @testset "Simple" begin
+        test_jackknife(SimpleObservable, SimpleVectorObservable)
+    end
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index ebb83f5..87bb875 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -9,7 +9,8 @@ const Therm = MCS
 const alpha = 0.001
 
 @testset begin
-    filenames = ["classical.jl",
+    filenames = ["observable.jl",
+                 "classical.jl",
                  "quantum.jl",
                  "checkpoint.jl"]
     for filename in filenames