text:
"""
    _wrap(::Type{T}, x...)

Wrap value(s) into a data object of the given type, for use with associative combinations.
"""
_wrap(::Type{T}, x::T) where {T} = x

"""
    _unfiltered(op) -> Bool

Returns true iff `_should_tick` will always return true.
"""
_unfiltered(::NodeOp) = false

"""
    _should_tick(op, data) -> Bool

This should be defined for any op that does not have `_unfiltered(op)` returning true.
The return value determines whether a knot should be emitted for this value.
"""
function _should_tick end

"""
    _combine(op, data_1, data_2) -> Data

This should be defined for all inception and windowed ops. Given two data objects, combine
them into a new data object.
"""
function _combine end

"""
    _update(op, data_1, x...) -> Data

Given a data object, and a new observation (of potentially multiple arguments), generate a
new data field.

If only creating an `InceptionOp`, it is sufficient to define this instead of `_combine`.
By default this will use `_wrap` and `_combine`.
"""
_update(op, data_1, x...) = _combine(op, data_1, _wrap(_data_type(op), x...))

"""
    _extract(op, data) -> value

This should be defined for all inception and windowed ops. Given some data object, it should
compute the appropriate output value for the node.
"""
function _extract end

"""
    _data_type(op)

Return the type of data used for the given op.
"""
function _data_type end

"""Unary operator accumulated from inception."""
abstract type UnaryInceptionOp{T,Data} <: UnaryNodeOp{T} end

"""Binary operator accumulated from inception."""
abstract type BinaryInceptionOp{T,Data,A} <: BinaryNodeOp{T,A} end

const InceptionOp{T,Data} = Union{UnaryInceptionOp{T,Data},BinaryInceptionOp{T,Data}}
_data_type(::InceptionOp{T,Data}) where {T,Data} = Data

always_ticks(op::InceptionOp) = _unfiltered(op)
time_agnostic(::InceptionOp) = true

mutable struct InceptionOpState{Data} <: NodeEvaluationState
    initialised::Bool
    data::Data
    # `data` will be uninitialised until the first call.
    InceptionOpState{Data}() where {Data} = new{Data}(false)
end

function create_operator_evaluation_state(
    ::Tuple{Node}, ::UnaryInceptionOp{T,Data}
) where {T,Data}
    return InceptionOpState{Data}()
end

function create_operator_evaluation_state(
    ::Tuple{Node,Node}, ::BinaryInceptionOp{T,Data}
) where {T,Data}
    return InceptionOpState{Data}()
end

function operator!(
    op::InceptionOp{T,Data}, state::InceptionOpState{Data}, x...
) where {T,Data}
    if !state.initialised
        state.data = _wrap(Data, x...)
        state.initialised = true
    else
        state.data = _update(op, state.data, x...)
    end
    return if always_ticks(op)
        # Deal with the case where we always emit.
        _extract(op, state.data)
    elseif _unfiltered(op) || _should_tick(op, state.data)
        Maybe(_extract(op, state.data))
    else
        Maybe{T}()
    end
end

"""
Windowed associative binary operator, potentially emitting early before the window is full.
"""
abstract type UnaryWindowOp{T,Data,EmitEarly} <: UnaryNodeOp{T} end

abstract type BinaryWindowOp{T,Data,EmitEarly,A} <: BinaryNodeOp{T,A} end

const WindowOp{T,Data,EmitEarly} = Union{
    UnaryWindowOp{T,Data,EmitEarly},BinaryWindowOp{T,Data,EmitEarly}
}
_data_type(::WindowOp{T,Data}) where {T,Data} = Data

"""
    _window(window_op) -> Int64

Return the window for the specified op.
The default implementation expects a field called `window` on the op structure.
"""
_window(op::WindowOp) = op.window

"""Whether or not this window op is set to emit with a non-full window."""
function _emit_early(::WindowOp{T,Data,true}) where {T,Data}
    return true
end
function _emit_early(::WindowOp{T,Data,false}) where {T,Data}
    return false
end

always_ticks(op::WindowOp) = _emit_early(op) && _unfiltered(op)
time_agnostic(::WindowOp) = true

mutable struct WindowOpState{Data} <: NodeEvaluationState
    window_state::FixedWindowAssociativeOp{Data}
end

function create_operator_evaluation_state(::Tuple{Node}, op::UnaryWindowOp)
    return create_operator_evaluation_state(op)
end

function create_operator_evaluation_state(::Tuple{Node,Node}, op::BinaryWindowOp)
    return create_operator_evaluation_state(op)
end

function create_operator_evaluation_state(op::WindowOp{T,Data}) where {T,Data}
    return WindowOpState{Data}(
        FixedWindowAssociativeOp{Data,(x, y) -> _combine(op, x, y)}(_window(op))
    )
end

function operator!(op::WindowOp{T,Data}, state::WindowOpState{Data}, x...) where {T,Data}
    update_state!(state.window_state, _wrap(Data, x...))
    if always_ticks(op)
        # Deal with the case where we always emit.
        return _extract(op, window_value(state.window_state))
    end

    ready = _emit_early(op) || window_full(state.window_state)
    if !ready
        return Maybe{T}()
    end

    data = window_value(state.window_state)
    return if _unfiltered(op) || _should_tick(op, data)
        Maybe(_extract(op, data))
    else
        Maybe{T}()
    end
end

# Sum, cumulative over time.
struct Sum{T} <: UnaryInceptionOp{T,T} end
_unfiltered(::Sum) = true
_combine(::Sum, x, y) = x + y
_extract(::Sum, data) = data
Base.show(io::IO, ::Sum{T}) where {T} = print(io, "Sum{$T}")
"""
    sum(x::Node) -> Node

Create a node which ticks when `x` ticks, with values of the cumulative sum of `x`.
"""
function Base.sum(x::Node)
    _is_constant(x) && return x
    return obtain_node((x,), Sum{value_type(x)}())
end

# Sum over fixed window.
struct WindowSum{T,EmitEarly} <: UnaryWindowOp{T,T,EmitEarly}
    window::Int64
end
_unfiltered(::WindowSum) = true
_combine(::WindowSum, x, y) = x + y
_extract(::WindowSum, data) = data
Base.show(io::IO, op::WindowSum{T}) where {T} = print(io, "WindowSum{$T}($(_window(op)))")

"""
    sum(x::Node, window::Int; emit_early::Bool=false) -> Node

Create a node of the rolling sum of `x` over the last `window` knots.

If `emit_early` is false, then the node returned will only start ticking once the window is
full. Otherwise, it will tick immediately with a partially-filled window.
"""
function Base.sum(x::Node, window::Int; emit_early::Bool=false)
    return obtain_node((x,), WindowSum{value_type(x),emit_early}(window))
end

# Product, cumulative over time.
struct Prod{T} <: UnaryInceptionOp{T,T} end
_unfiltered(::Prod) = true
_combine(::Prod, x, y) = x * y
_extract(::Prod, data) = data
Base.show(io::IO, ::Prod{T}) where {T} = print(io, "Prod{$T}")
"""
    prod(x::Node) -> Node

Create a node which ticks when `x` ticks, with values of the cumulative product of `x`.
"""
function Base.prod(x::Node)
    _is_constant(x) && return x
    return obtain_node((x,), Prod{value_type(x)}())
end

# Product over fixed window.
struct WindowProd{T,EmitEarly} <: UnaryWindowOp{T,T,EmitEarly}
    window::Int64
end
_unfiltered(::WindowProd) = true
_combine(::WindowProd, x, y) = x * y
_extract(::WindowProd, data) = data
Base.show(io::IO, op::WindowProd{T}) where {T} = print(io, "WindowProd{$T}($(_window(op)))")
"""
    prod(x::Node, window::Int; emit_early::Bool=false) -> Node

Create a node of the rolling product of `x` over the last `window` knots.

If `emit_early` is false, then the node returned will only start ticking once the window is
full. Otherwise, it will tick immediately with a partially-filled window.
"""
function Base.prod(x::Node, window::Int; emit_early::Bool=false)
    return obtain_node((x,), WindowProd{value_type(x),emit_early}(window))
end

# Mean, cumulative over time.
# In order to be numerically stable, use a generalisation of Welford's algorithm.
# Disable formatting: https://github.com/domluna/JuliaFormatter.jl/issues/480
#! format: off
const MeanData{T} = @NamedTuple{n::Int64, mean::T} where {T}
#! format: on
_wrap(::Type{MeanData{T}}, x) where {T} = MeanData{T}((1, x))
struct Mean{T} <: UnaryInceptionOp{T,MeanData{T}} end
_unfiltered(::Mean) = true
function _combine(state_a::MeanData{T}, state_b::MeanData{T}) where {T}
    na = state_a.n
    nb = state_b.n
    nc = na + nb
    return MeanData{T}((n=nc, mean=state_a.mean * (na / nc) + state_b.mean * (nb / nc)))
end
_combine(::Mean, x::MeanData, y::MeanData) = _combine(x, y)
_extract(::Mean, data::MeanData) = data.mean
Base.show(io::IO, ::Mean{T}) where {T} = print(io, "Mean{$T}")
"""
    mean(x::Node) -> Node

Create a node which ticks when `x` ticks, with values of the running mean of `x`.
"""
function Statistics.mean(x::Node)
    _is_constant(x) && return x
    T = output_type(/, value_type(x), Int)
    return obtain_node((x,), Mean{T}())
end

# Mean over fixed window.
struct WindowMean{T,EmitEarly} <: UnaryWindowOp{T,MeanData{T},EmitEarly}
    window::Int64
end
_unfiltered(::WindowMean) = true
_combine(::WindowMean, x::MeanData, y::MeanData) = _combine(x, y)
_extract(::WindowMean, data::MeanData) = data.mean
Base.show(io::IO, op::WindowMean{T}) where {T} = print(io, "WindowMean{$T}($(_window(op)))")
"""
    mean(x::Node, window::Int; emit_early::Bool=false) -> Node

Create a node of the rolling mean of `x` over the last `window` knots.

If `emit_early` is false, then the node returned will only start ticking once the window is
full. Otherwise, it will tick immediately with a partially-filled window.
"""
function Statistics.mean(x::Node, window::Int; emit_early::Bool=false)
    T = output_type(/, value_type(x), Int)
    return obtain_node((x,), WindowMean{T,emit_early}(window))
end

# Variance, cumulative over time.
# In order to be numerically stable, use a generalisation of Welford's algorithm.
# Disable formatting: https://github.com/domluna/JuliaFormatter.jl/issues/480
#! format: off
const VarData{T} = @NamedTuple{n::Int64, mean::T, s::T} where {T}
#! format: on
_wrap(::Type{VarData{T}}, x) where {T} = VarData{T}((1, x, 0))
struct Var{T,corrected} <: UnaryInceptionOp{T,VarData{T}} end
_should_tick(::Var, data::VarData) = data.n > 1
function _combine(state_a::VarData{T}, state_b::VarData{T}) where {T}
    na = state_a.n
    nb = state_b.n
    nc = na + nb

    μa = state_a.mean
    μb = state_b.mean
    μc = μa * (na / nc) + μb * (nb / nc)

    sa = state_a.s
    sb = state_b.s

    return VarData{T}((n=nc, mean=μc, s=(sa + sb) + nb * (μb - μa) * (μb - μc)))
end
_combine(::Var, x::VarData, y::VarData) = _combine(x, y)
_extract(::Var{T,true}, data::VarData) where {T} = data.s / (data.n - 1)
_extract(::Var{T,false}, data::VarData) where {T} = data.s / data.n
Base.show(io::IO, ::Var{T}) where {T} = print(io, "Var{$T}")
"""
    var(x::Node; corrected::Bool=true) -> Node

Create a node which ticks when `x` ticks, with values of the running variance of `x`.

This is equivalent to a sample variance over the `n` values of `x` observed at and before a
given time. If `corrected` is `true` (the default), we normalise by `n-1`, otherwise we
normalise by `n`.
"""
function Statistics.var(x::Node; corrected::Bool=true)
    _is_constant(x) && throw(ArgumentError("Cannot compute variance of constant $x"))
    T = output_type(/, value_type(x), Int)
    return obtain_node((x,), Var{T,corrected}())
end

# Variance over fixed window.
struct WindowVar{T,Corrected,EmitEarly} <: UnaryWindowOp{T,VarData{T},EmitEarly}
    window::Int64
end
_should_tick(::WindowVar, data::VarData) = data.n > 1
_combine(::WindowVar, x::VarData, y::VarData) = _combine(x, y)
_extract(::WindowVar{T,true}, data::VarData) where {T} = data.s / (data.n - 1)
_extract(::WindowVar{T,false}, data::VarData) where {T} = data.s / data.n
Base.show(io::IO, op::WindowVar{T}) where {T} = print(io, "WindowVar{$T}($(_window(op)))")
"""
    var(x::Node, window::Int; emit_early::Bool=false, corrected::Bool=true) -> Node

Create a node of the rolling variance of `x` over the last `window` knots.

If `emit_early` is false, then the node returned will only start ticking once the window is
full. Otherwise, it will tick immediately with a partially-filled window.

This is equivalent to a sample variance over the `n` values of `x` (capped at `window`),
observed at and before a given time. If `corrected` is `true` (the default), we normalise by
`n-1`, otherwise we normalise by `n`.
"""
function Statistics.var(x::Node, window::Int; emit_early::Bool=false, corrected::Bool=true)
    _is_constant(x) && throw(ArgumentError("Cannot compute variance of constant $x"))
    window >= 2 || throw(ArgumentError("Got window=$window, but should be at least 2"))
    T = output_type(/, value_type(x), Int)
    return obtain_node((x,), WindowVar{T,corrected,emit_early}(window))
end

# Standard deviation.
"""
    std(x::Node; corrected::Bool=true) -> Node

Create a node which ticks when `x` ticks, with values of the running standard deviation of
`x`.

This is equivalent to `sqrt(var(x; corrected))`.
"""
Statistics.std(x::Node; corrected::Bool=true) = sqrt(var(x; corrected))

"""
    std(x::Node, window::Int; emit_early::Bool=false, corrected::Bool=true) -> Node

Create a node of the rolling standard deviation of `x` over the last `window` knots.

If `emit_early` is false, then the node returned will only start ticking once the window is
full. Otherwise, it will tick immediately with a partially-filled window.

This is equivalent to `sqrt(var(x; emit_early, corrected))`.
"""
function Statistics.std(x::Node, window::Int; emit_early::Bool=false, corrected::Bool=true)
    return sqrt(var(x, window; emit_early, corrected))
end

# Covariance, cumulative over time.
# Disable formatting: https://github.com/domluna/JuliaFormatter.jl/issues/480
#! format: off
const CovData{T} = @NamedTuple{n::Int64, mean_x::T, mean_y::T, c::T} where {T}
#! format: on
_wrap(::Type{CovData{T}}, x, y) where {T} = CovData{T}((1, x, y, 0))
struct Cov{T,corrected,A} <: BinaryInceptionOp{T,CovData{T},A} end
_should_tick(::Cov, data::CovData) = data.n > 1
function _combine(state_a::CovData{T}, state_b::CovData{T}) where {T}
    na = state_a.n
    nb = state_b.n
    nc = na + nb

    μxa = state_a.mean_x
    μxb = state_b.mean_x
    μxc = μxa * (na / nc) + μxb * (nb / nc)

    μya = state_a.mean_y
    μyb = state_b.mean_y
    μyc = μya * (na / nc) + μyb * (nb / nc)

    ca = state_a.c
    cb = state_b.c
    # FIXME This is speculation - do the algebra to check this!
    cc = (ca + cb) + nb * (μxb - μxa) * (μyb - μyc)

    return CovData{T}((n=nc, mean_x=μxc, mean_y=μyc, c=cc))
end
_combine(::Cov, x::CovData, y::CovData) = _combine(x, y)
_extract(::Cov{T,true}, data::CovData) where {T} = data.c / (data.n - 1)
_extract(::Cov{T,false}, data::CovData) where {T} = data.c / data.n
Base.show(io::IO, ::Cov{T}) where {T} = print(io, "Cov{$T}")
"""
    cov(x, y[, alignment]; corrected::Bool=true) -> Node

Create a node which ticks with values of the running covariance of `x` and `y`.

The specified `alignment` controls the behaviour when `x` and `y` tick at different times,
as per the documentation in [Alignment](@ref). When not specified, it defaults to
[`UNION`](@ref).

This is equivalent to a sample covariance over the `n` values of `(x, y)` pairs observed at
and before a given time. If `corrected` is `true` (the default), we normalise by `n-1`,
otherwise we normalise by `n`.
"""
function Statistics.cov(x, y, ::A; corrected::Bool=true) where {A<:Alignment}
    x = _ensure_node(x)
    y = _ensure_node(y)
    if _is_constant(x) && _is_constant(y)
        throw(ArgumentError("Cannot compute variance of constants $x and $y"))
    end
    T = output_type(/, output_type(*, value_type(x), value_type(y)), Int)
    return obtain_node((x, y), Cov{T,corrected,A}())
end
function Statistics.cov(x::Node, y::Node; corrected::Bool=true)
    return cov(x, y, DEFAULT_ALIGNMENT; corrected)
end
Statistics.cov(x::Node, y; corrected::Bool=true) = cov(x, y, DEFAULT_ALIGNMENT; corrected)
Statistics.cov(x, y::Node; corrected::Bool=true) = cov(x, y, DEFAULT_ALIGNMENT; corrected)

# Covariance over fixed window.
struct WindowCov{T,Corrected,EmitEarly,A} <: BinaryWindowOp{T,CovData{T},EmitEarly,A}
    window::Int64
end
_should_tick(::WindowCov, data::CovData) = data.n > 1
_combine(::WindowCov, x::CovData, y::CovData) = _combine(x, y)
_extract(::WindowCov{T,true}, data::CovData) where {T} = data.c / (data.n - 1)
_extract(::WindowCov{T,false}, data::CovData) where {T} = data.c / data.n
Base.show(io::IO, op::WindowCov{T}) where {T} = print(io, "WindowCov{$T}($(_window(op)))")

"""
    cov(x, y, window::Int[, alignment]; emit_early=false, corrected=true) -> Node

Create a node of the rolling covariance of `x` and `y` over the last `window` knots.

The specified `alignment` controls the behaviour when `x` and `y` tick at different times,
as per the documentation in [Alignment](@ref). When not specified, it defaults to
[`UNION`](@ref).

If `emit_early` is false, then the node returned will only start ticking once the window is
full. Otherwise, it will tick immediately with a partially-filled window.

This is equivalent to a sample covariance over the `n` values of `(x, y)` pairs observed at
and before a given time, with `n` capped at `window`. If `corrected` is `true` (the
default), we normalise by `n-1`, otherwise we normalise by `n`.
"""
function Statistics.cov(
    x, y, window::Int, ::A; emit_early::Bool=false, corrected::Bool=true
) where {A<:Alignment}
    window >= 2 || throw(ArgumentError("Got window=$window, but should be at least 2"))
    x = _ensure_node(x)
    y = _ensure_node(y)
    if _is_constant(x) && _is_constant(y)
        throw(ArgumentError("Cannot compute covariance of constants $x and $y"))
    end
    T = output_type(/, output_type(*, value_type(x), value_type(y)), Int)
    return obtain_node((x, y), WindowCov{T,corrected,emit_early,A}(window))
end
function Statistics.cov(
    x::Node, y::Node, window::Int; emit_early::Bool=false, corrected::Bool=true
)
    return cov(x, y, window, DEFAULT_ALIGNMENT; emit_early, corrected)
end
function Statistics.cov(
    x::Node, y, window::Int; emit_early::Bool=false, corrected::Bool=true
)
    return cov(x, y, window, DEFAULT_ALIGNMENT; emit_early, corrected)
end
function Statistics.cov(
    x, y::Node, window::Int; emit_early::Bool=false, corrected::Bool=true
)
    return cov(x, y, window, DEFAULT_ALIGNMENT; emit_early, corrected)
end

# Covariance matrix.
abstract type AbstractCovMatrixData{T} end

"""State for an n-dimensional covariance matrix."""
struct CovMatrixData{T} <: AbstractCovMatrixData{T}
    n::Int64
    μ::Vector{T}
    c::Matrix{T}
end

function _wrap(::Type{CovMatrixData{T}}, x::AbstractVector) where {T}
    return CovMatrixData{T}(1, Vector(x), zeros(eltype(x), size(x, 1), size(x, 1)))
end

"""State for an n-dimensional covariance matrix of statically known dimension."""
struct CovMatrixStaticData{T,N,L} <: AbstractCovMatrixData{T}
    n::Int64
    μ::SVector{N,T}
    c::SMatrix{N,N,T,L}
end

function _wrap(::Type{CovMatrixStaticData{T,N,L}}, x::SVector{N}) where {T,N,L}
    return CovMatrixStaticData{T,N,L}(1, x, zeros(SMatrix{N,N,T,L}))
end

function _combine(a::Data, b::Data) where {Data<:AbstractCovMatrixData}
    n = a.n + b.n
    μ = @. a.μ * (a.n / n) + b.μ * (b.n / n)
    # TODO wrap c in `Hermitian` / only update one triangle?
    c = (a.c .+ b.c) .+ b.n .* (b.μ .- a.μ) * (b.μ .- μ)'
    return Data(n, μ, c)
end

# Static covariance matrix node ops.
struct CovMatrixStatic{N,L,T,corrected} <:
       UnaryInceptionOp{SMatrix{N,N,T,L},CovMatrixStaticData{T,N,L}} end
_should_tick(::CovMatrixStatic, data::CovMatrixStaticData) = data.n > 1
_combine(::CovMatrixStatic, x::CovMatrixStaticData, y::CovMatrixStaticData) = _combine(x, y)
function _extract(::CovMatrixStatic{N,L,T,true}, data::CovMatrixStaticData) where {N,L,T}
    return data.c ./ (data.n - 1)
end
function _extract(::CovMatrixStatic{N,L,T,false}, data::CovMatrixStaticData) where {N,L,T}
    return data.c ./ data.n
end
function Base.show(io::IO, ::CovMatrixStatic{N,L,T}) where {N,L,T}
    return print(io, "CovMatrixStatic{$T,$(N)x$(N)}")
end

"""
    cov(x::Node{<:AbstractVector}; corrected::Bool=true) -> Node

Create a node which ticks with values of the running covariance of `x`.

!!! warning
    If the values of `x` change shape over time, this node will throw an exception during
    evaluation.
"""
function Statistics.cov(
    x::Node{<:StaticVector{N,T}}; corrected::Bool=true
) where {N,T<:Number}
    _is_constant(x) && throw(ArgumentError("Cannot compute covariance of constant $x"))
    Out = output_type(/, T, Int)
    return obtain_node((x,), CovMatrixStatic{N,N * N,Out,corrected}())
end

# An n-dimensional covariance matrix over a fixed window.
struct WindowCovMatrixStatic{N,L,T,Corrected,EmitEarly} <:
       UnaryWindowOp{SMatrix{N,N,T,L},CovMatrixStaticData{T,N,L},EmitEarly}
    window::Int64
end
_should_tick(::WindowCovMatrixStatic, data::CovMatrixStaticData) = data.n > 1
function _combine(::WindowCovMatrixStatic, x::CovMatrixStaticData, y::CovMatrixStaticData)
    return _combine(x, y)
end
function _extract(
    ::WindowCovMatrixStatic{N,L,T,true}, data::CovMatrixStaticData
) where {N,L,T}
    return data.c ./ (data.n - 1)
end
function _extract(
    ::WindowCovMatrixStatic{N,L,T,false}, data::CovMatrixStaticData
) where {N,L,T}
    return data.c ./ data.n
end
function Base.show(io::IO, op::WindowCovMatrixStatic{N,L,T}) where {N,L,T}
    return print(io, "WindowCovMatrixStatic{$T,$(N)x$(N)}($(_window(op)))")
end
"""
    cov(x::Node{<:AbstractVector}, window::Int; corrected::Bool=true) -> Node

Create a node of the rolling covariance of `x` over the last `window` knots.

If `emit_early` is false, then the node returned will only start ticking once the window is
full. Otherwise, it will tick immediately with a partially-filled window.

!!! warning
    If the values of `x` change shape over time, this node will throw an exception during
    evaluation.
"""
function Statistics.cov(
    x::Node{<:StaticVector{N,T}}, window::Int; emit_early::Bool=false, corrected::Bool=true
) where {N,T<:Number}
    _is_constant(x) && throw(ArgumentError("Cannot compute variance of constant $x"))
    window >= 2 || throw(ArgumentError("Got window=$window, but should be at least 2"))
    Out = output_type(/, T, Int)
    return obtain_node(
        (x,), WindowCovMatrixStatic{N,N * N,Out,corrected,emit_early}(window)
    )
end

# Non-static covariance matrix nodes.
struct CovMatrix{T,corrected} <: UnaryInceptionOp{Matrix{T},CovMatrixData{T}} end
_should_tick(::CovMatrix, data::CovMatrixData) = data.n > 1
_combine(::CovMatrix, x::CovMatrixData, y::CovMatrixData) = _combine(x, y)
_extract(::CovMatrix{T,true}, data::CovMatrixData) where {T} = data.c ./ (data.n - 1)
_extract(::CovMatrix{T,false}, data::CovMatrixData) where {T} = data.c ./ data.n
Base.show(io::IO, ::CovMatrix{T}) where {T} = print(io, "CovMatrix{$T,?x?}")

function Statistics.cov(
    x::Node{<:AbstractVector{T}}; corrected::Bool=true
) where {T<:Number}
    _is_constant(x) && throw(ArgumentError("Cannot compute covariance of constant $x"))
    Out = output_type(/, T, Int)
    return obtain_node((x,), CovMatrix{Out,corrected}())
end

struct WindowCovMatrix{T,Corrected,EmitEarly} <:
       UnaryWindowOp{Matrix{T},CovMatrixData{T},EmitEarly}
    window::Int64
end
_should_tick(::WindowCovMatrix, data::CovMatrixData) = data.n > 1
_combine(::WindowCovMatrix, x::CovMatrixData, y::CovMatrixData) = _combine(x, y)
_extract(::WindowCovMatrix{T,true}, data::CovMatrixData) where {T} = data.c ./ (data.n - 1)
_extract(::WindowCovMatrix{T,false}, data::CovMatrixData) where {T} = data.c ./ data.n
function Base.show(io::IO, op::WindowCovMatrix{T}) where {T}
    return print(io, "WindowCovMatrix{$T,?x?}($(_window(op)))")
end
function Statistics.cov(
    x::Node{<:AbstractVector{T}}, window::Int; emit_early::Bool=false, corrected::Bool=true
) where {T<:Number}
    _is_constant(x) && throw(ArgumentError("Cannot compute variance of constant $x"))
    window >= 2 || throw(ArgumentError("Got window=$window, but should be at least 2"))
    Out = output_type(/, T, Int)
    return obtain_node((x,), WindowCovMatrix{Out,corrected,emit_early}(window))
end

# Correlation, cumulative over time.

"""
    cor(x, y[, alignment]; corrected::Bool=true) -> Node

Create a node which ticks with values of the running correlation of `x` and `y`.

The specified `alignment` controls the behaviour when `x` and `y` tick at different times,
as per the documentation in [Alignment](@ref). When not specified, it defaults to
[`UNION`](@ref).

This is equivalent to a sample correlation over the `n` values of `(x, y)` pairs observed at
and before a given time.
"""
function Statistics.cor(x, y, alignment::Alignment)
    # Do alignment first so that calls to cov & std so that we only do aligning work once.
    # Calls to cov & std should then be (relatively) trivial.
    x, y = coalign(x, y, alignment)
    return cov(x, y) / (std(x) * std(y))
end
Statistics.cor(x::Node, y::Node) = cor(x, y, DEFAULT_ALIGNMENT)
Statistics.cor(x::Node, y) = cov(x, y, DEFAULT_ALIGNMENT)
Statistics.cor(x, y::Node) = cov(x, y, DEFAULT_ALIGNMENT)

# Correlation over fixed window.

"""
    cor(x, y, window::Int[, alignment]; emit_early=false) -> Node

Create a node of the rolling covariance of `x` and `y` over the last `window` knots.

The specified `alignment` controls the behaviour when `x` and `y` tick at different times,
as per the documentation in [Alignment](@ref). When not specified, it defaults to
[`UNION`](@ref).

If `emit_early` is false, then the node returned will only start ticking once the window is
full. Otherwise, it will tick immediately with a partially-filled window.

This is equivalent to a sample correlation over the `n` values of `(x, y)` pairs observed at
and before a given time, with `n` capped at `window`.
"""
function Statistics.cor(x, y, window::Int, alignment::Alignment; emit_early::Bool=false)
    x, y = coalign(x, y, alignment)
    return cov(x, y, window; emit_early) /
           (std(x, window; emit_early) * std(y, window; emit_early))
end
function Statistics.cor(x::Node, y::Node, window::Int; emit_early::Bool=false)
    return cor(x, y, window, DEFAULT_ALIGNMENT; emit_early)
end
function Statistics.cor(x::Node, y, window::Int; emit_early::Bool=false)
    return cor(x, y, window, DEFAULT_ALIGNMENT; emit_early)
end
function Statistics.cor(x, y::Node, window::Int; emit_early::Bool=false)
    return cor(x, y, window, DEFAULT_ALIGNMENT; emit_early)
end

# TODO for time windows, it is safe to dispatch on the same functions with a TimePeriod
#   In the future, if we want to support more general types of time, we could introduce a
#   TimeDelta wrapper, that could server to distinguish between an Int64 "time" and an Int64
#   "number of knots".

struct EmaData{T}
    weighted_sum::T
    weighted_count::Float64
end

_wrap(::Type{EmaData{T}}, x) where {T} = EmaData{T}(x, 1.0)
struct Ema{T} <: UnaryInceptionOp{T,EmaData{T}}
    α::Float64
    function Ema{T}(α::Float64) where {T}
        0 < α < 1 || throw(ArgumentError("α = $α is outside of the valid range."))
        return new{T}(α)
    end
end
_unfiltered(::Ema) = true
function _update(op::Ema{T}, a::EmaData{T}, x) where {T}
    return EmaData(x + (1 - op.α) * a.weighted_sum, 1 + (1 - op.α) * a.weighted_count)
end
_extract(::Ema, data::EmaData) = data.weighted_sum / data.weighted_count
Base.show(io::IO, ::Ema{T}) where {T} = print(io, "Ema{$T}")

"""
    ema(x::Node, α::AbstractFloat) -> Node
    ema(x::Node, w_eff::Integer) -> Node

Create a node which computes the exponential moving average of `x`.

The decay is specified either by `α`, which should satisfy `0 < α < 1`, or by `w_eff`, which
should be an integer greater than 1. If the latter is specified, then we compute
`α = 2 / (w_eff + 1)`.

For internal state ``s_t``, with ``s_0 = 0``, and resulting EMA series ``m_t``, this has the
form:

```math
\\begin{aligned}
s_t &= s_{t-1} + (1 - \\alpha) x_t \\\\
m_t &= \\frac{\\alpha s_t}{1 - (1 - \\alpha)^t}.
\\end{aligned}
```

For further information, see the notational conventions and discussion on
[Wikipedia](https://en.wikipedia.org/wiki/Moving_average#Exponential_moving_average). Note
that this function implements the variant including the correction for the initial
convergence problem.
"""
function ema(x::Node, α::AbstractFloat)
    T = output_type(*, value_type(x), Float64)
    return obtain_node((x,), Ema{T}(α))
end
function ema(x::Node, w_eff::Integer)
    w_eff > 1 || throw(ArgumentError("w_eff = $w_eff is out of range"))
    return ema(x, 2 / (w_eff + 1))
end


text:
import NbodyGradient: InitialConditions

#========== Integrator ==========#
abstract type AbstractIntegrator end

"""
    Integrator{T<:AbstractFloat}

Integrator. Used as a functor to integrate a [`State`](@ref).

# Fields
- `scheme::Function` : The integration scheme to use.
- `h::T` : Step size.
- `t0::T` : Initial time.
- `tmax::T` : Duration of simulation.
"""
mutable struct Integrator{T<:AbstractFloat, schemeT} <: AbstractIntegrator
    scheme::schemeT
    h::T
    t0::T
    tmax::T
end

Integrator(h::T,t0::T) where T<:AbstractFloat = Integrator(ahl21!,h,t0,t0+h)
Integrator(scheme,h::Real,t0::Real,tmax::Real) = Integrator(scheme,promote(h,t0,tmax)...)

# Default to ahl21!
Integrator(h::T,t0::T,tmax::T) where T<:AbstractFloat = Integrator(ahl21!,h,t0,tmax)

#========== State ==========#
abstract type AbstractState end

"""
    State{T<:AbstractFloat} <: AbstractState

Current state of simulation.

# Fields (relevant to the user)
- `x::Matrix{T}` : Positions of each body [dimension, body].
- `v::Matrix{T}` : Velocities of each body [dimension, body].
- `t::Vector{T}` : Current time of simulation.
- `m::Vector{T}` : Masses of each body.
- `jac_step::Matrix{T}` : Current Jacobian.
- `dqdt::Vector{T}` : Derivative with respect to time.
"""
struct State{T<:AbstractFloat} <: AbstractState
    x::Matrix{T}
    v::Matrix{T}
    t::Vector{T}
    m::Vector{T}
    jac_step::Matrix{T}
    dqdt::Vector{T}
    jac_init::Matrix{T}
    xerror::Matrix{T}
    verror::Matrix{T}
    dqdt_error::Vector{T}
    jac_error::Matrix{T}
    n::Int64
    pair::Matrix{Bool}

    rij::Vector{T}
    a::Matrix{T}
    aij::Vector{T}
    x0::Vector{T}
    v0::Vector{T}
    input::Vector{T}
    delxv::Vector{T}
    rtmp::Vector{T}
end

"""
    State(ic)

Constructor for [`State`](@ref) type.

# Arguments
- `ic::InitialConditions{T}` : Initial conditions for the system.
"""
function State(ic::InitialConditions{T}) where T<:AbstractFloat
    x,v,jac_init = init_nbody(ic)
    n = ic.nbody
    xerror = zeros(T,size(x))
    verror = zeros(T,size(v))
    jac_step = Matrix{T}(I,7*n,7*n)
    dqdt = zeros(T,7*n)
    dqdt_error = zeros(T,size(dqdt))
    jac_error = zeros(T,size(jac_step))
    pair = zeros(Bool,n,n)

    rij = zeros(T,3)
    a = zeros(T,3,n)
    aij = zeros(T,3)
    x0 = zeros(T,3)
    v0 = zeros(T,3)
    input = zeros(T,8)
    delxv = zeros(T,6)
    rtmp = zeros(T,3)
    return State(x,v,[ic.t0],ic.m,jac_step,dqdt,jac_init,xerror,verror,dqdt_error,jac_error,ic.nbody,
    pair,rij,a,aij,x0,v0,input,delxv,rtmp)
end

function set_state!(s_old::State{T},s_new::State{T}) where T<:AbstractFloat
    s_old.t .= s_new.t
    s_old.x .= s_new.x
    s_old.v .= s_new.v
    s_old.jac_step .= s_new.jac_step
    s_old.xerror .= s_new.xerror
    s_old.verror .= s_new.verror
    s_old.jac_error .= s_new.jac_error
    s_old.dqdt .= s_new.dqdt
    s_old.dqdt_error .= s_new.dqdt_error
    return
end

"""Shows if the positions, velocities, and Jacobian are finite."""
Base.show(io::IO,::MIME"text/plain",s::State{T}) where {T} = begin
    println(io,"State{$T}:");
    println(io,"Positions  : ", all(isfinite.(s.x)) ? "finite" : "infinite!");
    println(io,"Velocities : ", all(isfinite.(s.v)) ? "finite" : "infinite!");
    println(io,"Jacobian   : ", all(isfinite.(s.jac_step)) ? "finite" : "infinite!");
    return
end

#========== Running Methods ==========#

"""
    (::Integrator)(s, time; grad=true)

Callable [`Integrator`](@ref) method. Integrate to specific time.

# Arguments
- `s::State{T}` : The current state of the simulation.
- `time::T` : Time to integrate to.

### Optional
- `grad::Bool` : Choose whether to calculate gradients. (Default = true)
"""
function (intr::Integrator)(s::State{T},time::T;grad::Bool=true) where T<:AbstractFloat
    t0 = s.t[1]

    # Calculate number of steps
    nsteps = abs(round(Int64,(time - t0)/intr.h))

    # Step either forward or backward
    h = intr.h * check_step(t0,time)

    # Calculate last step (if needed)
    #while t0 + (h * nsteps) <= time; nsteps += 1; end
    tmax = t0 + (h * nsteps)

    # Preallocate struct of arrays for derivatives
    if grad; d = Derivatives(T,s.n); end

    for i in 1:nsteps
        # Take integration step and advance time
        if grad
            intr.scheme(s,d,h)
        else
            intr.scheme(s,h)
        end
    end

    # Do last step (if needed)
    if nsteps == 0; hf = time; end
    if tmax != time
        hf = time - tmax
        if grad
            intr.scheme(s,d,hf)
        else
            intr.scheme(s,hf)
        end
    end

    s.t[1] = time
    return
end

"""
    (::Integrator)(s, N; grad=true)

Callable [`Integrator`](@ref) method. Integrate for N steps.

# Arguments
- `s::State{T}` : The current state of the simulation.
- `N::Int64` : Number of steps.

### Optional
- `grad::Bool` : Choose whether to calculate gradients. (Default = true)
"""
function (intr::Integrator)(s::State{T},N::Int64;grad::Bool=true) where T<:AbstractFloat
    s2 = zero(T) # For compensated summation

    # Preallocate struct of arrays for derivatives
    if grad; d = Derivatives(T,s.n); end

    # check to see if backward step
    if N < 0; intr.h *= -1; N *= -1; end
    h = intr.h

    for n in 1:N
        # Take integration step and advance time
        if grad
            intr.scheme(s,d,h)
        else
            intr.scheme(s,h)
        end
        s.t[1],s2 = comp_sum(s.t[1],s2,intr.h)
    end

    # Return time step to forward if needed
    if intr.h < 0; intr.h *= -1; end
    return
end

"""
    (::Integrator)(s; grad=true)

Callable [`Integrator`](@ref) method. Integrate from `s.t` to `tmax` -- specified in constructor.

# Arguments
- `s::State{T}` : The current state of the simulation.

### Optional
- `grad::Bool` : Choose whether to calculate gradients. (Default = true)
"""
(intr::Integrator)(s::State{T};grad::Bool=true) where T<:AbstractFloat = intr(s,s.t[1]+intr.tmax,grad=grad)

function check_step(t0::T,tmax::T) where T<:AbstractFloat
    if abs(tmax) > abs(t0)
        return sign(tmax)
    else
        if sign(tmax) != sign(t0)
            return sign(tmax)
        else
            return -1 * sign(tmax)
        end
    end
end

#========== Includes  ==========#
const ints = ["ahl21"]
for i in ints; include(joinpath(i,"$i.jl")); end

const ints_no_grad = ["ahl21","dh17"]
for i in ints_no_grad; include(joinpath(i,"$(i)_no_grad.jl")); end

ext: jl


text:
@testset "contexts.jl" begin
    @testset "PrefixContext" begin
        ctx = @inferred PrefixContext{:f}(
            PrefixContext{:e}(
                PrefixContext{:d}(
                    PrefixContext{:c}(
                        PrefixContext{:b}(PrefixContext{:a}(DefaultContext()))
                    ),
                ),
            ),
        )
        vn = VarName{:x}()
        vn_prefixed = @inferred DynamicPPL.prefix(ctx, vn)
        @test DynamicPPL.getsym(vn_prefixed) == Symbol("a.b.c.d.e.f.x")
        @test vn_prefixed.indexing === vn.indexing

        vn = VarName{:x}((1,))
        vn_prefixed = @inferred DynamicPPL.prefix(ctx, vn)
        @test DynamicPPL.getsym(vn_prefixed) == Symbol("a.b.c.d.e.f.x")
        @test vn_prefixed.indexing === vn.indexing
    end
end


text:
function Max()
    return Max(())
end

function value(obj::Max, arg0::jdouble, arg1::jdouble)
    return jcall(obj, "value", jdouble, (jdouble, jdouble), arg0, arg1)
end


text:
struct PSResult{T}
    value::T
end

struct PSError <: Exception
    message::String
end

# Abstract types used in other PS packages -
abstract type PSAbstractAsset end
abstract type PSAbstractAssetReturnModel end
abstract type PSAbstractAssetTreeNode end
abstract type PSAbstractAssetTree end
abstract type PSAbstractProblemModel end


text:
function time_stepping_works(arch, FT, Closure)
    model = Model(grid=RegularCartesianGrid(FT; size=(16, 16, 16), length=(1, 2, 3)),
                  architecture=arch, float_type=FT, closure=Closure(FT))
    time_step!(model, 1, 1)
    return true  # test that no errors/crashes happen when time stepping.
end

function run_first_AB2_time_step_tests(arch, FT)
    add_ones(args...) = 1.0
    model = Model(grid=RegularCartesianGrid(FT; size=(16, 16, 16), length=(1, 2, 3)),
                  architecture=arch, float_type=FT, forcing=ModelForcing(T=add_ones))
    time_step!(model, 1, 1)

    # Test that GT = 1 after first time step and that AB2 actually reduced to forward Euler.
    @test all(interior(model.timestepper.Gⁿ.u) .≈ 0)
    @test all(interior(model.timestepper.Gⁿ.v) .≈ 0)
    @test all(interior(model.timestepper.Gⁿ.w) .≈ 0)
    @test all(interior(model.timestepper.Gⁿ.T) .≈ 1.0)
    @test all(interior(model.timestepper.Gⁿ.S) .≈ 0)

    return nothing
end

"""
    This test ensures that when we compute w from the continuity equation that the full velocity field
    is divergence-free.
"""
function compute_w_from_continuity(arch, FT)
    Nx, Ny, Nz = 16, 16, 16
    Lx, Ly, Lz = 16, 16, 16

    grid = RegularCartesianGrid(FT; size=(Nx, Ny, Nz), length=(Lx, Ly, Lz))
    bcs = HorizontallyPeriodicSolutionBCs()

    u = FaceFieldX(FT, arch, grid)
    v = FaceFieldY(FT, arch, grid)
    w = FaceFieldZ(FT, arch, grid)
    div_u = CellField(FT, arch, grid)

    interior(u) .= rand(FT, Nx, Ny, Nz)
    interior(v) .= rand(FT, Nx, Ny, Nz)

    fill_halo_regions!(u.data, bcs.u, arch, grid)
    fill_halo_regions!(v.data, bcs.v, arch, grid)
    compute_w_from_continuity!((u=u.data, v=v.data, w=w.data), grid)

    fill_halo_regions!(w.data, bcs.w, arch, grid)
    velocity_div!(grid, u.data, v.data, w.data, div_u.data)

    # Set div_u to zero at the top because the initial velocity field is not divergence-free
    # so we end up some divergence at the top if we don't do this.
    interior(div_u)[:, :, Nz] .= zero(FT)

    min_div = minimum(interior(div_u))
    max_div = maximum(interior(div_u))
    sum_div = sum(interior(div_u))
    abs_sum_div = sum(abs.(interior(div_u)))
    @info "Velocity divergence after recomputing w ($arch, $FT): min=$min_div, max=$max_div, sum=$sum_div, abs_sum=$abs_sum_div"

    return all(isapprox.(interior(div_u), 0; atol=5*eps(FT)))
end

"""
    This tests to make sure that the velocity field remains incompressible (or divergence-free) as the model is time
    stepped. It just initializes a cube shaped hot bubble perturbation in the center of the 3D domain to induce a
    velocity field.
"""
function incompressible_in_time(arch, FT, Nt)
    Nx, Ny, Nz = 32, 32, 32
    Lx, Ly, Lz = 10, 10, 10

    model = Model(grid=RegularCartesianGrid(size=(Nx, Ny, Nz), length=(Lx, Ly, Lz)), architecture=arch, float_type=FT)

    grid = model.grid
    u, v, w = model.velocities.u, model.velocities.v, model.velocities.w

    div_u = CellField(FT, arch, model.grid)

    # Just add a temperature perturbation so we get some velocity field.
    @. model.tracers.T.data[8:24, 8:24, 8:24] += 0.01

    time_step!(model, Nt, 0.05)

    velocity_div!(grid, u, v, w, div_u)

    min_div = minimum(interior(div_u))
    max_div = minimum(interior(div_u))
    sum_div = sum(interior(div_u))
    abs_sum_div = sum(abs.(interior(div_u)))
    @info "Velocity divergence after $Nt time steps ($arch, $FT): min=$min_div, max=$max_div, sum=$sum_div, abs_sum=$abs_sum_div"

    # We are comparing with 0 so we use absolute tolerances. They are a bit larger than eps(Float64) and eps(Float32)
    # because we are summing over the absolute value of many machine epsilons. A better atol value may be
    # Nx*Ny*Nz*eps(FT) but it's much higher than the observed abs_sum_div.
    if FT == Float64
        return isapprox(abs_sum_div, 0; atol=5e-16)
    elseif FT == Float32
        return isapprox(abs_sum_div, 0; atol=1e-7)
    end
end

"""
    tracer_conserved_in_channel(arch, FT, Nt)

Create a super-coarse eddying channel model with walls in the y and test that
temperature is conserved after `Nt` time steps.
"""
function tracer_conserved_in_channel(arch, FT, Nt)
    Nx, Ny, Nz = 16, 32, 16
    Lx, Ly, Lz = 160e3, 320e3, 1024

    α = (Lz/Nz)/(Lx/Nx) # Grid cell aspect ratio.
    νh, κh = 20.0, 20.0
    νv, κv = α*νh, α*κh

    model = ChannelModel(architecture = arch, float_type = FT,
                         grid = RegularCartesianGrid(size=(Nx, Ny, Nz), length=(Lx, Ly, Lz)),
                         closure = ConstantAnisotropicDiffusivity(νh=νh, νv=νv, κh=κh, κv=κv))

    Ty = 1e-4  # Meridional temperature gradient [K/m].
    Tz = 5e-3  # Vertical temperature gradient [K/m].

    # Initial temperature field [°C].
    T₀(x, y, z) = 10 + Ty*y + Tz*z + 0.0001*rand()
    set_ic!(model, T=T₀)

    Tavg0 = mean(interior(model.tracers.T))

    time_step!(model; Nt=Nt, Δt=10*60)

    Tavg = mean(interior(model.tracers.T))
    @info "Tracer conservation after $Nt time steps ($arch, $FT): ⟨T⟩-T₀=$(Tavg-Tavg0) °C"

    # Interestingly, it's very well conserved (almost to machine epsilon) for
    # Float64, but not as close for Float32... But it does seem constant in time
    # for Float32 so at least it is bounded.
    if FT == Float64
        return isapprox(Tavg, Tavg0; rtol=2e-14)
    elseif FT == Float32
        return isapprox(Tavg, Tavg0; rtol=2e-4)
    end
end

Closures = (ConstantIsotropicDiffusivity, ConstantAnisotropicDiffusivity,
            ConstantSmagorinsky, AnisotropicMinimumDissipation)

@testset "Time stepping" begin
    println("Testing time stepping...")

    for arch in archs, FT in float_types, Closure in Closures
        @test time_stepping_works(arch, FT, Closure)
    end

    @testset "2nd-order Adams-Bashforth" begin
        println("  Testing 2nd-order Adams-Bashforth...")
        for arch in archs, FT in float_types
            run_first_AB2_time_step_tests(arch, FT)
        end
    end

    @testset "Recomputing w from continuity" begin
        println("  Testing recomputing w from continuity...")
        for arch in archs, FT in float_types
            @test compute_w_from_continuity(arch, FT)
        end
    end

    @testset "Incompressibility" begin
        println("  Testing incompressibility...")
        for arch in archs, FT in float_types, Nt in [1, 10, 100]
            @test incompressible_in_time(arch, FT, Nt)
        end
    end

    @testset "Tracer conservation in channel" begin
        println("  Testing tracer conservation in channel...")
        for arch in archs, FT in float_types
            @test tracer_conserved_in_channel(arch, FT, 10)
        end
    end
end


text:
# Copyright (c) 2020 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

is_transition_matrix(A::AbstractMatrix) =
(size(A,1) == size(A,2)) && all(abs.(sum(A,dims=2) .- 1) .< 1e-6)

function t2tr2t2IJ(t2tr)
  Z = [[k .* ones(Int,length(t2tr[t][k])) =>sort(t2tr[t][k]) for k in sort([keys(t2tr[t])...])] for t in eachindex(t2tr)]
  t2IJ=[vcat(getindex.(Z[t],1)...) => vcat(getindex.(Z[t],2)...) for t in eachindex(Z)]
  return t2IJ
end

function logadd(x::T, y::T) where {T<:AbstractFloat}
  if isinf(x) return y end
  if isinf(y) return x end
  if x < y
    diff = x-y
    x = y
  else
    diff = y-x
  end
  return x + log1p(exp(diff))
end


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ def88812-6b99-11eb-07fb-3f961d7ff6ae
begin
	using Plots
	plotly()
end

# ╔═╡ c3fd3b6c-6c33-11eb-3efe-57b3d5456387
using PlutoUI

# ╔═╡ e67f8692-799c-11eb-31bb-09a7f37da454
using StatsPlots

# ╔═╡ 9642a69e-7988-11eb-069e-157feefb4eed
using RDatasets

# ╔═╡ 88d2e984-7984-11eb-3f4c-5115712ac1a2
using Distributions

# ╔═╡ 1457d8d0-7b4b-11eb-3b05-fdea9ad03b02
using StatsBase

# ╔═╡ f054b9aa-7e63-11eb-36dd-3bf1932030ac
using QuadGK

# ╔═╡ 22f18e82-6b7d-11eb-16eb-31758abdb1b0
md"# Notebook 5: _Measures of Spread_
"

# ╔═╡ 59607216-7b43-11eb-2cf7-3f74231a1839
md"
> 📕Reference material - Chapter 3: Using statistics to summarise data (Sections 3.5 and 3.6) from Sheldon Ross.
"

# ╔═╡ 29c96178-7981-11eb-332f-95d53eceb696
md"So far we have been looking at statistics which capture central tendency. Now we will move on to statistics that capture spread of data."

# ╔═╡ 052cbe72-7984-11eb-0b02-1dd45dc17b3b
md"👉Why are measures of spread important?"

# ╔═╡ baa51422-6b93-11eb-16b7-bf0872dc2bea
md"## Range"

# ╔═╡ 449528c8-7983-11eb-1148-4fdaa380dcfd
md"The first natural statistic for spread of data is the difference between two sets of statistics"

# ╔═╡ 4d443e46-7983-11eb-33e3-316d616ad44b
md"The simplest type of range is the difference between maximum and minimum values."

# ╔═╡ 9119e782-7984-11eb-25cb-471cefe4e7d5
range_(x) = maximum(x) - minimum(x)

# ╔═╡ 9a2fe03e-7988-11eb-09ea-b709a8cbd77c
iris = dataset("datasets", "iris")

# ╔═╡ d781854e-7988-11eb-385e-e3dedbaae06e
range_(iris.PetalLength)

# ╔═╡ a70ee822-7988-11eb-112e-5be62d6310c4
range_(iris.PetalWidth)

# ╔═╡ ba3e633c-7988-11eb-2550-4b645fb2439d
range_(iris.SepalWidth)

# ╔═╡ dde99fd6-7988-11eb-061a-15e9399613b6
md"Let us look at range calculations for random variates sampled from uniform distribution"

# ╔═╡ a7f67e52-7984-11eb-22b1-69c114b468fb
uniform_sample = rand(Uniform(0, 1), 10)

# ╔═╡ bd098168-7986-11eb-110c-c5d17c258584
md"
```math
x_i \sim U(0, 1),\ i = 1, 2, \ldots, n
```
"

# ╔═╡ cf8e1c18-7984-11eb-20f9-3ff941f8b457
range_(uniform_sample)

# ╔═╡ 04a88244-7985-11eb-328d-6f3878193236
md"Notice that in the case of a uniform distribution, we know the range of the _population_ which is the *upper bound* on the range of the _sample_. The two values would get closer as we increase the sample size."

# ╔═╡ aba0e4f6-7985-11eb-0713-172db02e7224
md"👉 Probability exercise: Given that we are drawing $n$ random numbers from the uniform random distribution [0, 1], what is the probability that the maximum value is more than some $v \in (0, 1)$."

# ╔═╡ d961ff42-7985-11eb-37d2-35af70b59310
md"
```math
\mathrm{Pr}\large(\max(x_1, x_2, \ldots, x_n) > v\large) =  1 - \prod_{i = 1}^n \mathrm{Pr}(x_i \leq v) = 1 - v^n
```
"

# ╔═╡ 6636f208-7987-11eb-16a3-17e36a925c97
md"Since $v < 1$, the above probability tends to 1 as we increase $n$, i.e., as we increase the sample size, the maximum value gets closer to 1."

# ╔═╡ 1b6955c0-7986-11eb-180f-6fabd0e6cf36
md"👉 Another probability exercise: Given that we are drawing $n$ random numbers from the uniform random distribution [0, 1], what is the probability of the minimum being lesser than $u \in (0, 1)$."

# ╔═╡ e7e6e6de-7987-11eb-2fc7-37fab2d3c05b
md"
```math
\mathrm{Pr}\large(\min(x_1, x_2, \ldots, x_n) < u\large) =  1 - \prod_{i = 1}^n \mathrm{Pr}(x_i \geq u) = 1 - (1 - u)^n
```
"

# ╔═╡ eba873f4-7988-11eb-0c30-372c31f4cad8
md"Heads up -> You will see in the assignment a question that asks you to compute the probability of range of a sample being less than $d$."

# ╔═╡ 0697dcea-7989-11eb-0ead-8d17b8f4b6cd
md"We can do the same for other distributions like the normal distribution."

# ╔═╡ 47ae6534-7984-11eb-3be1-3501475c06f9
normal_sample = rand(Normal(), 100)

# ╔═╡ 98f6a88c-7984-11eb-3882-d121f8488565
range_(normal_sample)

# ╔═╡ f6fa502e-7cce-11eb-028e-c790c42e9bc2
md"👉💡 Compute the probability of bounds on maximum and minimum values in a set of $n$ random variates from the standard normal distribution"

# ╔═╡ 554abfc4-7ccf-11eb-25fb-df246eee94a0
md"👉 The range in terms of extrema values is susceptible to outliers, how would you fix this?"

# ╔═╡ 4a80d04c-7989-11eb-109b-4789a2c1480a
md"## Interquartile Range"

# ╔═╡ 87d63e9a-798a-11eb-393f-e59d4eb48bea
md"As we have seen, outliers affect the quality of statistics and the range statistic which depends on maximum and minimum values is highly susceptible."

# ╔═╡ 98bc8d9a-798a-11eb-1fc3-0bd77bb3fe8c
md"One easy way of fixing this is to discard a part of the extreme values and then take an effective range in the remaining values. Or differently stated, we can find the difference between two percentiles."

# ╔═╡ bc770cc6-798a-11eb-0333-c746ba22dc7f
md"A common choice is the **interquartile range**, that computes the difference between the first quartile (25-percentile) and third quartile (75-percentile). This is usually denoted as $\mathrm{IQR} = Q_3 - Q_1$"

# ╔═╡ 08692f0e-798b-11eb-2b02-9548f1a59426
x = rand(1:10, 11)

# ╔═╡ b0ea2b14-798c-11eb-0a6c-6f06820a0fa1
sort(x)

# ╔═╡ 15a6fa1e-798d-11eb-187e-b93de64da00b
md"Median is at position 6 -> $(sort(x)[6])"

# ╔═╡ 1aa85896-798d-11eb-15ff-7f64835af8c8
md"Q1 is average of values as positions 3, 4 -> $((sort(x)[3] + sort(x)[4])/2)"

# ╔═╡ 24b05384-798d-11eb-346a-77856d43c699
md"Q3 is average of values as positions 8, 9 -> $((sort(x)[8] + sort(x)[9])/2)"

# ╔═╡ 70ef989a-798d-11eb-1899-f51410a810eb
md"IQR = Q3 - Q1 = $((sort(x)[8] + sort(x)[9])/2 - (sort(x)[3] + sort(x)[4])/2)"

# ╔═╡ 970787c2-798d-11eb-12c7-bfd48299c383
boxplot(x)

# ╔═╡ 4b5fb2f4-798b-11eb-1af8-33890aaf0ed4
md"The highlighted box captures the middle 50 percentile of data and its height gives the IQR"

# ╔═╡ 242401ba-798b-11eb-31ed-07d5767cf6f2
boxplot(iris.PetalWidth)

# ╔═╡ 8071095e-798b-11eb-3e86-0b137ae8360f
boxplot(iris.SepalWidth)

# ╔═╡ 8842e8aa-798b-11eb-1746-f36c13e7c8ba
md"Note how the range for sepal and petal widths were the same, but IQR for sepal width is lesser indicating a very different data distribution. Also notice how there are _outliers_ shows as individual points beyond the _whiskers_ of the box plot."

# ╔═╡ 5d79259e-798b-11eb-298a-b58a6288162e
md"👉 How would you decide on the classification of outliers?"

# ╔═╡ 2f1b4990-798e-11eb-31fb-dbc8934ccdfb
md"The usual choice for the lower whisker (also called a fence) is $\max(P_0, Q_1 - 1.5 \times \mathrm{IQR})$ where $P_0$ is the minimum value in the sample."

# ╔═╡ 42b201a8-799d-11eb-00ce-27483756dce2
md"Similarly, the usual choice for the upper whisker is $\min(P_{100}, Q_3 + 1.5 \times IQR)$, where $P_{100}$ is the maximum value in the sample."

# ╔═╡ 673294b4-799d-11eb-07ac-9dbd30c54541
md"This value 1.5 is usually a parameter that can be specified:"

# ╔═╡ 6fce44da-799d-11eb-1be4-8b0d2097ce9a
boxplot(iris.SepalWidth, range = 1) # default value is 1.5

# ╔═╡ 92d81aa8-799d-11eb-25e3-7b2d722ba0c2
md"So the boxplot shows us the IQR range as a highlighted box and the whiskers as a more elongated range that also depends on the IQR."

# ╔═╡ 1082bef2-7989-11eb-0137-f36a82244d4f
md"## Variance and Standard Deviation"

# ╔═╡ 81db5570-79b2-11eb-2dcb-5b908ad62b3f
md"
| Symbol      | Description |
| ----------- | ----------- |
| $\overline{x}$         | Sample mean (statistic)     |
| $\mu$         | Population mean (parameter)      |
| $s^2$      | Sample variance (statistic) |
| $\sigma^2$ | Population variance (parameter) |
| $s$ | Sample standard deviation (statistic) |
| $\sigma$ | Population standard deviation (parameter) |
|
"

# ╔═╡ 0c4ac2ae-79b3-11eb-09be-b385fd8dd389
md"Let us start with variance. We want to get some notion of spread of data."

# ╔═╡ 5ae7f7e6-79b4-11eb-1609-470cf5246eb7
arr = rand(7)

# ╔═╡ 95d0b9b0-79b4-11eb-2120-978ee85a75a3
begin
	scatter(arr, zeros(7))
	scatter!([mean(arr)], [0.1])
end

# ╔═╡ 40841800-79b7-11eb-2667-25b9a152246d
md"We can think of spread as how far are the points from the mean. We already know that $\sum_{i=1}^n (x_i - \bar{x}) = 0$."

# ╔═╡ 5aec1d14-79b7-11eb-0342-997f50dbe294
md"One way to avoid the cancelling is to take a square, i.e., $\sum_{i = 1}^n (x_i - \bar{x})^2$."

# ╔═╡ b420b066-79b7-11eb-16f6-9fa333021ea9
md"But since these are positive summands, this sum will depend on $n$. Instead we may want an average value that does not depend on $n$."

# ╔═╡ c8fa957e-79b7-11eb-355c-db21c14bccdc
md"Hence we have the following formula"

# ╔═╡ f84b2414-79b3-11eb-29d0-a3da4a9d41be
md"
```math
\sigma^2 = \dfrac{\sum_{i = 1}^N(X_i - \mu)^2}{N}
```
"

# ╔═╡ 5be6eb56-79b0-11eb-10a2-2bd5979f3a63
md"
```math
s^2 = \dfrac{\sum_{i = 1}^n(x_i - \overline{x})^2}{n - 1}
```
"

# ╔═╡ 1b53b1cc-79bb-11eb-16c2-a799ed6bafbd
md"The variance is an average of the square distance of elements from the mean."

# ╔═╡ 1056b4f2-79b3-11eb-09a6-9f7ffaf3b5e3
md"Notice the $n - 1$ in the formula for the sample variance. We will later see why it is $n - 1$ and not $n$. Quick trailer: Are the $n$ summands in the computation of $s^2$ independent?"

# ╔═╡ b4645be0-7b49-11eb-1f50-9328aece496b
md"👉💡 What happens to the sample variance if all sample values are offset by a constant?"

# ╔═╡ 07703e92-7b4a-11eb-0aa6-7180b261728c
md"👉💡 What happens to the sample variance if all sample values are scaled by a constant?"

# ╔═╡ 7e72da22-7cd2-11eb-3de0-93e1435b1987
md"👉Reall that the function $\sum_i (x_i - u)^2$ is minimized at $u = \overline{x}$. Variance is the ratio of that minimum value to the size of the population."

# ╔═╡ a0b559f6-79bb-11eb-0bb5-eb2b8e341b3c
mean(arr)

# ╔═╡ 842b718a-79bb-11eb-15e3-fd52f3ed3be4
[arr (arr .- mean(arr)).^2]

# ╔═╡ 1d2029b0-7b44-11eb-2df3-cf8c03c6c009
variance_(arr) = sum((arr .- mean(arr)).^2)/(length(arr) - 1)

# ╔═╡ 369e4a72-7b44-11eb-2947-d198ea29598e
variance_(arr)

# ╔═╡ 549ce572-7b44-11eb-3ea9-3334c962ca90
mean_(arr) = sum(arr)/length(arr)

# ╔═╡ cc201ce0-7b44-11eb-37fc-358d02e2d8ed
md"The variance is a measure of the spread of the data while the mean is a central tendency. It is easy to construct two arrays with roughly the same mean but different variance."

# ╔═╡ 5a33dad6-7b44-11eb-36cb-034dc4851a01
arr_narrow = rand(-10:10, 10000)

# ╔═╡ 671c8c5a-7b44-11eb-0302-53101084c5a6
arr_wide = rand(-100:100, 10000)

# ╔═╡ 7195a696-7b44-11eb-1073-91f2060f0fd5
[mean_(arr_narrow) mean_(arr_wide)]

# ╔═╡ 2e18943a-7ff1-11eb-1d35-7f2019c6c06e
[variance_(arr_narrow) variance_(arr_wide)]

# ╔═╡ d4f04180-7b45-11eb-210a-75f898949c20
md"The following formula will come handy in many theorems later and is also useful if you are estimating variance by hand."

# ╔═╡ afde9040-801a-11eb-3907-4f9f38b66893
md"
```math
\sigma^2 = \dfrac{\sum_{i = 1}^N(X_i - \mu)^2}{N}
```
"

# ╔═╡ edf6aaae-7b45-11eb-3e5c-ff386f3ca855
md"
```math
\sigma^2 = \overline{X^2} - \mu^2
```
"

# ╔═╡ 9ebeb81e-7cd2-11eb-3842-67d666de0d56
md"👉 In class exercise. Try to work out the proof for this."

# ╔═╡ 176b9300-7b47-11eb-34f5-e9f401587fbb
md"Notice that these are population parameters and not sample statistics"

# ╔═╡ 10813aba-7b46-11eb-15ae-9d34277557aa
md"
```math
\sigma^2 = \dfrac{1}{N}\sum_{i = 1}^N (X_i - \mu)^2 = \dfrac{1}{N} \left(\sum_{i = 1}^N X_i^2+ N\mu^2 - 2\sum_{i = 1}^N X_i \mu\right)
```
```math
= \overline{X^2} + \mu^2 - \dfrac{2\mu}{N}\sum_{i = 1}^NX_i = \overline{X^2} - \mu^2
```
"

# ╔═╡ 522289b6-7b47-11eb-2e36-f1fcaad71e46
md"This is also written as $\sigma^2 = \mathbb{E}[X^2] - {\mathbb{E}[X]}^2$."

# ╔═╡ 09b45022-7b49-11eb-3a09-25e2c5f15d5a
md"Let us build some visual intuition for this formula"

# ╔═╡ dfaca6b8-7b47-11eb-3232-c395a968eab0
arr1 = [1, 3, 5]

# ╔═╡ ce181838-7b47-11eb-09fb-21ae2437e92f
begin
	scatter(arr1, zeros(length(arr1)), xlims = (0,27), label="Original elements")
	scatter!([mean(arr1)], [0.01], label="Mean")
	scatter!(arr1.^2, zeros(length(arr1)) .+ 0.1, label="Squares of elements")
	scatter!([mean(arr1)^2], [0.2], label="Square of mean")
	scatter!([mean(arr1.^2)], [0.2], label="Mean of squares")
end

# ╔═╡ c3f0eb9c-7e5c-11eb-11ab-b39896644d0e
md"Squaring followed by mean is larger than mean followed by squaring. The difference captures the variance in data."

# ╔═╡ 1cf14bf8-7b4a-11eb-057c-cd6c391812d7
md"The formula $\sigma^2 = \overline{X^2} - \mu^2$ can be written for the sample variance also as
```math
s^2 = \dfrac{n}{n - 1} \left(\overline{x^2} - \overline{x}^2\right)
```
"

# ╔═╡ 8fe26172-7b4a-11eb-3fec-25a023ad0c21
md"Notice that the difference between dividing by $n$ and $n- 1$ in the formula for the standard deviation is significant for small values of $n$. For large samples (say $n = 10,000$), the difference is negligible."

# ╔═╡ 413e485c-7b4b-11eb-212e-81acd087a606
md"The variance function is typically denoted as `var` and is available from the `StatsBase` package.)"

# ╔═╡ 17a48a4c-7b4b-11eb-0935-c78e3fbe8875
[var(arr) variance_(arr)]

# ╔═╡ 5f6eae0c-7b4b-11eb-1f34-891b85a352b3
md"Notice that by default, `var` is the sample variance (with the $n-1$ denominator). If you are computing population variance, then you can set the `corrected` argument to false."

# ╔═╡ 843168c4-7b4b-11eb-2d50-f50df5c611b5
var(arr, corrected=false)

# ╔═╡ 93ba3102-7b4b-11eb-1e40-4d22cef75d41
md"Most numerical packages have some way of specifying whether it is a sample statistic or population parameter. It is crucial that you see what the default setting is. This is a common gotcha."

# ╔═╡ 5508697c-7e5d-11eb-18c8-0f27ee5cb6be
md"Now a further complication arises when we compute variance of a sample given that we know the population mean."

# ╔═╡ 754bdb0e-7e5d-11eb-24ea-a15f33748d49
md"
If we did not know $\mu$:
```math
s^2 = \dfrac{\sum_{i = 1}^n(x_i - \overline{x})^2}{n - 1}
```
If we knew $\mu$:
```math
s^2 = \dfrac{\sum_{i = 1}^n(x_i - \mu)^2}{n}
```
"

# ╔═╡ a39d3d72-7e5d-11eb-3fa8-b57f67497390
md"So, if we are using the _noisy_ estimate of $\mu$ in the form of $\overline{x}$, then we have the _corrected_ denominator $n - 1$. But if we knew $\mu$, we use the standard denominator $n$. We will build more formalism to understand this later."

# ╔═╡ 8f75f92a-7ff1-11eb-07db-07a8b1ee8fd3
md"👉 Can you think of situations where we know $\mu$ and want to estimate $s^2$?"

# ╔═╡ c16cb486-7e5d-11eb-3931-616157e00d6d
md"👉 Can you devise a way to build computational intuition of the above formulae."

# ╔═╡ cb9c1654-7e5d-11eb-12a4-b1b814fb894a
arr_n = rand(Normal(), 10)

# ╔═╡ 222b0534-7e5e-11eb-07ee-599111a5a85e
mu = 0

# ╔═╡ 3c7d66e8-7e5e-11eb-156b-bf16ab97429b
σ2 = 1

# ╔═╡ 1ba1e9a6-7e5e-11eb-26c0-c3259ca71dfa
x_bar = sum(arr_n)/length(arr_n)

# ╔═╡ 5802e8a2-7e5e-11eb-0e69-2d5da350e203
s2_wo_σ2 = sum((arr_n .- x_bar).^2)/(length(arr_n) - 1)

# ╔═╡ 72d949fa-7e5e-11eb-3661-8f14b46a3336
s2_with_σ2 = sum((arr_n .- mu).^2)/(length(arr_n))

# ╔═╡ a1d2283a-7e5e-11eb-0c22-a91232575d2e
md"
**In summary**

For population parameter:
```math
\sigma^2 = \dfrac{\sum_{i = 1}^N(X_i - \mu)^2}{N}
```
For sample statistic if we did not know $\mu$:
```math
s^2 = \dfrac{\sum_{i = 1}^n(x_i - \overline{x})^2}{n - 1}
```
For sample statistic if we knew $\mu$:
```math
s^2 = \dfrac{\sum_{i = 1}^n(x_i - \mu)^2}{n}
```
"

# ╔═╡ 01c7f2d0-7e5d-11eb-0134-b527c31615cd
md"👉 How would you compute variance with frequency tables?"

# ╔═╡ eff8182c-7b4c-11eb-3608-077859fd9405
md"Instead of individual values, if we are given a frequency table where element $x_i$ occurs $f_i$ number of times, then the sample variance is given as
```math
s^2 = \dfrac{\sum_{i = 1}^n \left(f_i \times (x_i - \overline{x}^2)\right)}{\sum_{i = 1}^n f_i - 1}
```
"

# ╔═╡ 3795afc2-7e5f-11eb-04dd-650633462059
md"👉How would you compute variance from probability distribution functions"

# ╔═╡ 2f9ccf56-7b4d-11eb-0620-85cc73e63986
md"We can extend this to a probability distribution with a given pdf $f$. Notice that probability distribution characterises a population and not a sample"

# ╔═╡ 5723264a-7b4d-11eb-3d80-d728e4e125ea
md"
```math
\sigma^2 = \int_{-\infty}^{\infty} (x - \mu)^2 f(x) dx
```
where $\mu$ as we know is
```math
\mu = \int_{-\infty}^{\infty} xf(x) dx
```
"

# ╔═╡ 393eab66-7e60-11eb-138b-b7b5d7578b9e
md"The above formulation is related to the idea of moments of functions.

The $n$th moment of a function $f$ for a given constant $c$ is:
```math
M_n = \int_{-\infty}^{\infty} (x - c)^n f(x) dx
```
"

# ╔═╡ 50be8e64-7e60-11eb-0bde-770db17176c5
md"Mean is the first moment ($n = 1$) of the pdf for $c = 0$."

# ╔═╡ 8af58b82-7e60-11eb-345b-f5497a167d14
md"Variance is the second moment $(n = 2)$ of the pdf for $c = \mu$. This is usually called the corrected or central moment."

# ╔═╡ 5d793888-7b4d-11eb-1760-7d6e5a06d299
 md"👉 For the standard normal distribution, what is $\sigma^2$"

# ╔═╡ 45d7729e-7b4e-11eb-16b8-9fe470c4422a
md"
```math
\sigma^2 = \int_{-\infty}^\infty x^2 \dfrac{1}{\sqrt{2\pi}} e^{\frac{-x^2}{2}} dx
```
We will invoke _integration by parts_.
"

# ╔═╡ 8f23a268-7b52-11eb-14d3-a19f581d662e
md"
Let us set $u$ and $v$ as: $u = x, dv = x e^{\frac{-x^2}{2}} dx$.
We may set $v = - e^{\frac{-x^2}{2}}$.
```math
\dfrac{1}{\sqrt{2\pi}} \left(\int_{-\infty}^{\infty} x^2 e^{\frac{-x^2}{2}} dx\right) = \dfrac{1}{\sqrt{2\pi}}  \left(-xe^{\frac{-x^2}{2}} \Biggr|_{-\infty}^{\infty} + \int_{-\infty}^{\infty} - e^{\frac{-x^2}{2}} dx\right)
```
```math
= 0 + 1 = 1
```
For the first part, note that the $e^{-x^2/2}$ decays much faster than $x$ grows. So both limits are 0. For the second part note the famous integral $\int_{-\infty}^{\infty} e^{-x^2} dx = \sqrt{\pi}$.
"

# ╔═╡ e62703f0-7e60-11eb-061b-05e2ebbae775
md"👉Homework: Compute variance analytically for $\mathcal{N}(\mu, \sigma)$."

# ╔═╡ d06cb7e4-7b45-11eb-0ecb-3fb5bd004188
md"## Standard deviation"

# ╔═╡ e7ee031a-7b44-11eb-1109-6d05b4ddca2e
md"Notice the variance is the sum of square of deivations. It's unit is thus the square of the unit of the values in the sample. If we want a measure of spread as the same unit as the sample data, then we can use the standard deviation which is simply the square-root of the variance."

# ╔═╡ 2806eb6a-7b45-11eb-227a-93964c77e93b
std_(arr) = sqrt(variance_(arr))

# ╔═╡ 89e42a9c-7b44-11eb-315d-7974de3531d9
[std_(arr_narrow) std_(arr_wide)]

# ╔═╡ 4c62ba78-7e62-11eb-080e-9f49bb8a60fb
md"
For population parameter:
```math
\sigma = \sqrt{\dfrac{\sum_{i = 1}^n(X_i - \mu)^2}{N}}
```
For sample statistic if we did not know $\mu$:
```math
s = \sqrt{\dfrac{\sum_{i = 1}^n(x_i - \overline{x})^2}{n - 1}}
```
For sample statistic if we knew $\mu$:
```math
s = \sqrt{\dfrac{\sum_{i = 1}^n(x_i - \mu)^2}{n}}
```
"

# ╔═╡ 7953f36c-7e62-11eb-318c-c30dc4038629
md"Notice that it the standard deviation is **not** the average of the square-root of sum of squares. It is the square-root of the avearge of the sum of squares."

# ╔═╡ fc298e8e-7e62-11eb-1275-bfaaba9e6d13
md"Because the standard deviation is the same unit as the data, we use it for analysis more often than variance. We can compute it for the petal lengths from before."

# ╔═╡ 35da102c-7e61-11eb-3068-3f8961843a9f
var(iris.PetalLength)

# ╔═╡ 5d6f375c-7e61-11eb-28ab-4561d1cabdb9
mean(iris.SepalLength)

# ╔═╡ 501021f2-7e61-11eb-3205-03823f96e0ba
std(iris.PetalWidth)

# ╔═╡ 14bf236c-7e63-11eb-31c0-1b57f03ae4ff
md"Since the standard deviation is of the same unit as the values of the sample, we can use it to define intervals of values. For instance $[\mu - \sigma, \mu + \sigma]$ defines an interval of size $2\sigma$ centered at the mean."

# ╔═╡ 711c29be-7e65-11eb-0aad-8f82e60e2ece
md"👉 A good question to ask is how many points are within this interval. Notice that counting points is the forte of another statistic - inter-quartile range, but we would like to see if we can use mean and standard deviation for this."

# ╔═╡ 5ceecd68-7e63-11eb-20a2-cf9cdf7e1758
md"Let us consider the normal distribution and identify the percentage of values than lie in the range $[\mu - k\times \sigma, \mu + k\times \sigma]$ for different values of $k$"

# ╔═╡ fab8dd52-7e65-11eb-3456-8708ccd8cd4f
md"👉How would you do this computationally?"

# ╔═╡ 9c434dd6-7e63-11eb-3bba-736c1f794722
D_ = Beta(2, 5)

# ╔═╡ 7b34bb66-7e63-11eb-2e44-1b00e23874cb
xslider = @bind x_slider html"<input type=range min=1 max=5 step=1>"

# ╔═╡ 660f6682-7e63-11eb-396d-77ef6e6f1de2
begin
	plot(x->x, x->pdf(D_, x), mean(D_) - x_slider*std(D_), mean(D_) + x_slider*std(D_), fill=(0, :orange),fillalpha=0.5, label="", title="Within [μ - " * string(x_slider) * "σ, μ +" * string(x_slider) * "σ]: " * string(floor(quadgk(x -> pdf(D_, x), (mean(D_) - x_slider*std(D_), mean(D_) + x_slider*std(D_))...)[1] * 10000)/100) * "%")
	plot!(x->x, x->pdf(D_, x), mean(D_) - 5 * std(D_), mean(D_) + 5 * std(D_), label="", line=4)
end

# ╔═╡ 624943a8-7e66-11eb-12b0-6f9a43405d16
md"There are two major implications of this:

1. If we have a window of size $2\sigma$ around $\mu$, more than 95% of the samples would lie within this.
2. If we have a sample value that is greater than $\mu + 2\sigma$, we can conclude that it has a low probability of being sampled.
"

# ╔═╡ 8848bee0-7e6a-11eb-0cb3-63001bba774d
md"👉 Examples?"

# ╔═╡ 91787884-7e6a-11eb-3821-a11457bb4bfa
md"
1. If I am a shoe retailer, then I would stock shoe sizes in the range $[\mu - 2\sigma, \mu + 2\sigma]$.
2. If I see a student scoring significantly less than $\mu - 2\sigma$, I could infer that there was some specific issue which may rquire me to support the student with."

# ╔═╡ db73f682-7e6a-11eb-2890-adea119688ab
md"👉How do we compare the standard deviations of different variables?"

# ╔═╡ d4666004-7e6a-11eb-1502-41d2bdf3df30
std(iris.PetalLength)

# ╔═╡ fec9f55c-7e6a-11eb-337c-13273fff4450
std(iris.PetalWidth)

# ╔═╡ 01e8107c-7e6b-11eb-1e0b-fd862028439a
md"👉Can we conclude that petal length has a higher *spread* than petal width?"

# ╔═╡ 258692fe-7e6b-11eb-3401-e9e8d96483fe
md"No, the mean values themselves may be different."

# ╔═╡ 2b9af64e-7e6b-11eb-11be-714f382922dc
md"A common approach is to compute the ratio of standard deviation to mean. This is called the **coefficient of variation**."

# ╔═╡ 0fc0b998-7e6b-11eb-34be-8be4f0e4fe08
std(iris.PetalLength)/mean(iris.PetalLength)

# ╔═╡ 190b6e28-7e6b-11eb-069c-dbaa847caa19
std(iris.PetalWidth)/mean(iris.PetalWidth)

# ╔═╡ 3e31d85c-7e6b-11eb-0cb8-6d171248d19d
md"So it turns out that coefficient of variation is higher for petal width, contrary to our earlier guess."

# ╔═╡ 2cc216d4-7e6c-11eb-09c8-ddaa8b9b05b5
md"👉What happens to the standard deviation if we have a constant shift of all values and a constant scaling of all values?"

# ╔═╡ 99172a5c-7e6c-11eb-39db-b76d5eb6403f
md"👉Do you see a relation between the median/IQR and mean/standard-deviation?"

# ╔═╡ e700935c-7e6c-11eb-3c5f-79ee52818f27
md"👉 When will you use which?"

# ╔═╡ ec198862-7e6c-11eb-2368-cb1d288a537a
md"If the data is skewed, better to use median and IQR. This is usually the case with economic data like house prices and income. For more symmetric data, mean and standard deviation are used. This is usually the case with experimental data."

# ╔═╡ 0121fce4-7e6d-11eb-00a3-77cef143d4e0
md"Notice: Measure of central tendency required computation of spread to check if it was accurate enough. But now, we seem to also depend on **skew**. How is that measured? Is it a measure of spread?"

# ╔═╡ 272137ca-7e6d-11eb-1ccb-491acdf93e80
md"## Skew"

# ╔═╡ 5d67d74a-6b93-11eb-0df5-01899af4b4ab
md"Let us revisit our earlier plots on how the mean, median, and mode are ordered."

# ╔═╡ a53d2cd8-7e6e-11eb-1265-c3df598243b6
function dist_plot(density, xlim)
	plot(xlim[1]:(xlim[2]-xlim[1])/100:xlim[2], [pdf(density, x) for x in xlim[1]:(xlim[2]-xlim[1])/100:xlim[2]], label=false, line=3)
	d_mean = mean(density)
	d_mode = mode(density)
	d_median = median(density)
	plot!([d_mean, d_mean], [0, pdf(density, d_mean)], label="Mean", line=(4, :dash, :green))
	plot!([d_mode, d_mode], [0, pdf(density, d_mode)], label="Mode", line=(1, :red))
	plot!([d_median, d_median], [0, pdf(density, d_median)], label="Median", line=(3, :dot, :orange))
end

# ╔═╡ add3881e-7e74-11eb-24d0-c3909978a450
md"👉Think of two distributions with the same mean, but different standard deviation"

# ╔═╡ bbf9ee7a-7e6e-11eb-06d9-21df405b04a2
dist_plot(Normal(0, 1), [-5, 5])

# ╔═╡ ba25c50c-7e74-11eb-248f-95e03691af90
dist_plot(Normal(0, 2), [-5, 5])

# ╔═╡ c2e4c756-7e74-11eb-330c-db0b6ea51886
md"👉Think of two distributions with the same mean and standard deviation, but clearly different characteristics"

# ╔═╡ c8c0d22e-7e6e-11eb-1e6d-37acc66bdadb
dist_plot(Beta(2, 5), [0, 1])

# ╔═╡ ca2cce42-7e6e-11eb-1b9a-8958575b107e
dist_plot(Beta(5, 2), [0, 1])

# ╔═╡ 1024a3ec-7e70-11eb-3b96-391f330cfb32
[std(Beta(2, 5)) std(Beta(5,2))]

# ╔═╡ 0e46c762-7e70-11eb-2650-135cff7436b8
md"Notice that the standard deviation for the two distributions is the same. The mean can be made the same by simply adding an offset to all values. Then, how are the two distributions to be distinguished?"

# ╔═╡ 30f7ecf0-7e70-11eb-13d5-75ace1ff42cd
md"We need a metric for skew"

# ╔═╡ 371558a2-7e70-11eb-1661-430f14e76302
md"Pearson's coefficient of skew: (recall that it was Pearson who proposed the definition of standard deviation too)
```math
\mathrm{Skew} = \dfrac{\mathrm{Mean} - \mathrm{Mode}}{\mathrm{Standard\  deviation}}
```
He also proposed
```math
\mathrm{Skew} = \dfrac{3(\mathrm{Mean} - \mathrm{Median})}{\mathrm{Standard\  deviation}}
```
The constant $3$ occurs because of sufficiently skewed data, mode $\approx$ $3 \times $ Median - $2 \times$ Mean
"

# ╔═╡ 907116e8-7e70-11eb-333f-290f0a7d80e8
skew1(d_) = (mean(d_) - mode(d_))/std(d_)

# ╔═╡ 4766721c-7e71-11eb-0dd7-69fdac58cc30
skew1.([Beta(2, 5) Beta(5, 2)])

# ╔═╡ 661feb52-7e71-11eb-1bf0-5b4a4f4e1c02
skew2(d_) = 3*(mean(d_) - median(d_))/std(d_)

# ╔═╡ 719b971a-7e71-11eb-1a2d-e1b67112c405
skew2.([Beta(2, 5) Beta(5, 2)])

# ╔═╡ 56e6b8f2-7e71-11eb-1203-75faf28f79e0
md"Now we have a measure to distinguish between these two distributions"

# ╔═╡ b3419d86-7e71-11eb-1a63-87efbffc5d90
md"If we change the two parameters of the Beta distribution, we can get different skews. Having greater difference between the parameters leads to higher absolute values of skewness."

# ╔═╡ c7d99398-7e71-11eb-2a80-e1ee835a675d
skew1.([Beta(1.1, 5) Beta(5, 1.1)])

# ╔═╡ cc00e586-7e71-11eb-05b9-c92bdd03dc26
dist_plot(Beta(1.1, 2), [0, 1])

# ╔═╡ e3b93278-7e71-11eb-01fb-c1fb19c1af97
md"Symmetric distributions have 0 skew"

# ╔═╡ 8c02ca6a-7e71-11eb-0533-f56d312294ac
skew1.([Normal() Uniform() SymTriangularDist()])

# ╔═╡ d9d0162a-7e72-11eb-143c-a7a0fdaa02e4
md"We can also calculate the skew using a moment-based approach."

# ╔═╡ 629dae04-7e73-11eb-1075-69cffc768b82
md"
For continuous pdf $f$:
```math
M_3 = \dfrac{1}{\sigma^3}\int_{-\infty}^{\infty} (x - \mu)^3 f(x) dx
```
And for finite sums:
```math
M_3 = \dfrac{\sum_{i = 1}^N (X_i - \mu)^3}{N\sigma^3}
```
"

# ╔═╡ 6e49f3b0-7ff3-11eb-00b5-2157914eb9ed
md"Notice that we are scaling the moment by $\sigma^3$."

# ╔═╡ da31badc-7e73-11eb-34a0-75d05bc41291
md"For samples when we do not know the population mean and standard deviation, we again need a *correction*:
```math
\mathrm{Skew} = \dfrac{n}{(n - 1)(n - 2)}\dfrac{\sum_{i = 1}^n (x_i - \overline{x})^3}{s^3}
```
"

# ╔═╡ 0d7fbfbe-7e91-11eb-14eb-ed9b68f7b02e
md"👉 What is the intuition of sign of skew?"

# ╔═╡ 6be14370-7e70-11eb-36c1-97d6b5ab2e2a
md"![](https://www.oreilly.com/library/view/clojure-for-data/9781784397180/graphics/7180OS_01_180.jpg)"

# ╔═╡ 377f3020-7e74-11eb-22d1-3bf79ca13e7a
skew_moment(x) = length(x)/((length(x) - 1) * (length(x) - 2)) * sum((x .- mean(x)).^3) / std(x)^3

# ╔═╡ 649cd9e0-7e74-11eb-0997-4f7390f13936
skew_moment(iris.PetalWidth)

# ╔═╡ 89350110-7e74-11eb-0cb1-6bb00c61e12b
histogram(iris.PetalWidth, bins=20)

# ╔═╡ b39bb5ac-7e81-11eb-3da9-fb2f6ece1f8a
diamonds = dataset("ggplot2", "diamonds");

# ╔═╡ b5d346f4-7e81-11eb-15a4-657b672f23a5
skew_moment(diamonds.Price)

# ╔═╡ bf5aaed6-7e81-11eb-3036-cbe315f2a59a
density(diamonds.Price, line=4)

# ╔═╡ fa35967e-7e81-11eb-05bc-ab321c3635c0
md"Building an intuition for skew - What values are moderately and what are highly skewed?"

# ╔═╡ 2c5cdce8-7e82-11eb-3f31-eb1c12783dad
betaslider = @bind beta_slider html"<input type=range min=1 max=10 step=0.1>"

# ╔═╡ 7f9a4eb2-7e83-11eb-2667-7d3c451bfc09
D_skew = Weibull(beta_slider, 1)

# ╔═╡ 70a7f836-7e82-11eb-367d-4125cd359a72
skewness(D_skew)

# ╔═╡ 40d3b36a-7e82-11eb-140c-9757638394e7
dist_plot(D_skew, [0, 3])

# ╔═╡ c080c712-7e83-11eb-1143-158e83f386e2
md"Rule of thumb:
$0 \leq |\mathrm{skewness}| < 0.5$ is approximately symmetric.
$0.5 \leq |\mathrm{skewness}| < 1$ is moderatedly skewed.
$|\mathrm{skewness}| > 1$ is highly skewed.
"

# ╔═╡ 2b426de4-7e84-11eb-277f-41bf0b1ea003
md"#### Kurtosis"

# ╔═╡ 330784ce-7e84-11eb-3513-374b23f1c836
md"👉Can you think of two distributions with the same mean, standard deviation, and skew, but still being different?"

# ╔═╡ 42f59704-7e84-11eb-1410-e3edb4a255ad
dist_plot(Normal(0, sqrt(5/3)), [-6, 6])

# ╔═╡ 22147248-7e85-11eb-1f56-17469bc96208
dist_plot(TDist(5), [-6, 6])

# ╔═╡ 4b9a1c3a-7e85-11eb-36b8-9b92628a80eb
mean.([Normal(0, sqrt(5/3)) TDist(5)])

# ╔═╡ 5bb290de-7e85-11eb-16a9-45daa0302e09
std.([Normal(0, sqrt(5/3)) TDist(5)])

# ╔═╡ 6023d310-7e85-11eb-0a26-9d85c0398e18
skewness.([Normal(0, sqrt(5/3)) TDist(5)])

# ╔═╡ 97d9cc9c-7e85-11eb-1e65-6d5d3fb245a7
md"👉What is different about the two distributions?"

# ╔═╡ 9f460e70-7e85-11eb-2c83-776eb9a8a401
begin
	plot(-6:0.1:6, [pdf(Normal(0, sqrt(5/3)), x) for x in -6:0.1:6], label="Normal", line=3)
	plot!(-6:0.1:6, [pdf(TDist(5), x) for x in -6:0.1:6], label="T Distribution", line=3)
end

# ╔═╡ edb98cd0-7e85-11eb-0c7a-57d0fef4c7b3
md"We have taken some area from the _shoulders_ of the normal distribution and added it to the center. We have made the distribution more **peaky** from the normal distribution to the student distribution."

# ╔═╡ dee27332-7e8c-11eb-3896-2566e284cbcd
md"Though the peaky middle is visually arresting, note that the tails of the student distribution are heavier than the tails of the normal distribution."

# ╔═╡ fafcb580-7e93-11eb-139c-5574a72195f1
begin
	plot(5:0.1:10, [pdf(Normal(0, sqrt(5/3)), x) for x in 5:0.1:10], label="Normal", line=3)
	plot!(5:0.1:10, [pdf(TDist(5), x) for x in 5:0.1:10], label="T Distribution", line=3)
end

# ╔═╡ 1328a088-7e94-11eb-1893-8307ec0a0354
md"It is this **heavy tail** behavior that we are trying to capture with kurtosis. More on this shortly."

# ╔═╡ 272b2e92-7e86-11eb-1905-b3c59f94c4b8
md"This heavy tail behavior can be quantified with the fourth centralised and scaled moment."

# ╔═╡ 67d8500a-7e86-11eb-0cfb-916b7df799eb
md"
For continuous pdf $f$:
```math
M_4 = \dfrac{1}{\sigma^4}\int_{-\infty}^{\infty} (x - \mu)^4 f(x) dx
```
And for finite sums:
```math
M_4 = \dfrac{\sum_{i = 1}^N (X_i - \mu)^4}{N\sigma^4}
```
"

# ╔═╡ 06c81c4a-7e87-11eb-2b8b-49c3e2185e5d
md"This measure is called **Kurtosis** which comes for the Greek root for arched or bulging."

# ╔═╡ 699f856c-7e94-11eb-017f-b57a7bddc047
md"👉 Why does the kurtosis capture tails better than the values close to mean?"

# ╔═╡ 44f4d7ac-7e95-11eb-0bd1-d3eb98605f56
md"This has been a confusion in general and is well documented in [this](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4321753/#:~:text=The%20incorrect%20notion%20that%20kurtosis,rest%20once%20and%20for%20all.) paper helpfully titled _Kurtosis as Peakedness, 1905 – 2014. R.I.P._"

# ╔═╡ f08b0aae-7e87-11eb-2044-f59130ffc70a
kurtosis_(d) = round(quadgk(x -> (x - mean(d))^4 * pdf(d, x), -100, 100)[1] / var(d)^2 * 1000)/1000

# ╔═╡ 7aa15fca-7e8c-11eb-3fd3-93d2e6eef233
kurtosis_.([Normal(0, sqrt(5/3)) TDist(5)])

# ╔═╡ d98262b2-7e8c-11eb-3bf3-5fbf2fa3aaae
md"The higher value of kurtosis, the greater the concentration in the middle and at the tails."

# ╔═╡ 2385d402-7e8d-11eb-07bf-17752bfcc9b6
md"👉 For the skew we had 0 denoting symmetric distribution. How would you define a simliar benchmark for kurtosis?"

# ╔═╡ 5b8f4b08-7e8d-11eb-0039-f5a7fa194ba1
md"Use the kurtosis of the standard normal distribution as a baseline."

# ╔═╡ a615ea5c-7e8c-11eb-2346-a9a9b35a584a
kurtosis_(Normal())

# ╔═╡ 86d649e2-7e8d-11eb-3b53-47bb7f2c6a57
md"Since 3 is an artibrary number, it is a convention to subtract it from the computation of kurtosis so that the **excess kurtosis** for the standard normal distribution is 0. Positive values of excess kurtosis denote heavier tails while negative values denote lighter tails."

# ╔═╡ ce4c4f58-7e8d-11eb-3332-ab9b8d8148b9
md"The `kurtosis` function from `StatsBase.jl` computes the excess kurtosis by default"

# ╔═╡ c9b9ebe2-7e8d-11eb-0890-b3e23ceb3a83
kurtosis(Normal(0, sqrt(5/3)))

# ╔═╡ adc658c8-7e90-11eb-32d9-eb4c98405392
md"For samples when we do not know the population mean and standard deviation, we again need a *correction*:
```math
\mathrm{Kurtosis} = \dfrac{n(n+1)}{(n - 1)(n - 2)(n-3)}\dfrac{\sum_{i = 1}^n (x_i - \overline{x})^4}{s^4} - \dfrac{3(n - 1)^2}{(n - 2)(n - 3)}
```
"

# ╔═╡ 469a8e50-7e93-11eb-21b0-d1df5dfb4b0e
md"The second term is to ensure that we are computing the excess kurtosis"

# ╔═╡ 68e0762c-6c33-11eb-3d59-37d618e5c795
hint(text, title) = Markdown.MD(Markdown.Admonition("hint", title, [text]))

# ╔═╡ 1167a94a-7984-11eb-3af5-c1290bd2781c
hint(md"Measures of spread describe data. They also tell us how useful measures of central tendency are. For instance if a distribution of a data is very wide, then the mean may not be a useful representation of all values in the sample.", "Measures of spread in relation to measures of central tendency")

# ╔═╡ c2f6ff4e-7fef-11eb-0d7e-bd34ab2132af
hint(md"The probability that the maximum of $n$ variates is more than $v$ is given by $(1 - (\mathrm{cdf}(v))^n)$", "Bound on maxima")

# ╔═╡ 502ef336-7b4e-11eb-1984-4fb3fc2696a9
hint(md"$\int u dv = uv - \int v du$", "What is integration by parts")

# ╔═╡ 3b9ea1c0-7e6c-11eb-1cbf-dfe3d0fbc9c1
hint(md"
* Constant offset -> No change in standard deviation.
* Constant scaling by $c$ -> standard deviation scales by $c$.
", "Effect of transformations on standard deviation")

# ╔═╡ a641d362-7e6c-11eb-3a34-e727f9afbaa7
hint(md"The IQR is the measure of spread when the measure of central tendency is median. The IQR gives the interval around the median which contains exactly 50% of the values.", "Median and IQR")

# ╔═╡ bb96eaa6-7e6c-11eb-3251-7dfe84995533
hint(md"The standard deviation is the measure of spread when the measure of central tendency is the mean. For the normal distribution, an interval of length twice the standard deviation around the mean contains about 68% of the values.", "Mean and Standard deviation")

# ╔═╡ 0970b240-7e73-11eb-2ee9-f391e64d584b
hint(md"So far: $M_0$ is 1. $M_1$ is $\mu$. $M_2$ about $\mu$ is $\sigma^2$.", "The story of moments so far")

# ╔═╡ 3fe26ca6-7e91-11eb-1382-dd95809c4187
hint(md"Break up $\sum (X_i - \mu)^3$ into two parts: values $X_i$ that are less than $\mu$ and contribute a negative term and the values $X_i$ that are larger than $\mu$ and contribute a positive term. Negative skew implies that the former outweigh the latter.", "Sign of skew")

# ╔═╡ 7c0b35ac-7e94-11eb-19b8-3fa36803e662
hint(md"
```math
M_4 = \dfrac{\sum_{i = 1}^N (X_i - \mu)^4}{N\sigma^4} = \dfrac{N\sum_{i = 1}^N (X_i - \mu)^4}{\left(\sum_{i = 1}^N (X_i - \mu)^2\right)^2}
```
The outliers have an outsided role in the numerator and outweight the denominator (recall our difference between $\mathbb{E}[X^2]$ and $\mathbb{E}[X]^2$). Hence, the Kurtosis is decided by heavy tails and not by peaky centers.", "Kurtosis and peaks")

# ╔═╡ Cell order:
# ╟─22f18e82-6b7d-11eb-16eb-31758abdb1b0
# ╟─59607216-7b43-11eb-2cf7-3f74231a1839
# ╠═def88812-6b99-11eb-07fb-3f961d7ff6ae
# ╠═c3fd3b6c-6c33-11eb-3efe-57b3d5456387
# ╠═e67f8692-799c-11eb-31bb-09a7f37da454
# ╟─29c96178-7981-11eb-332f-95d53eceb696
# ╟─052cbe72-7984-11eb-0b02-1dd45dc17b3b
# ╟─1167a94a-7984-11eb-3af5-c1290bd2781c
# ╟─baa51422-6b93-11eb-16b7-bf0872dc2bea
# ╟─449528c8-7983-11eb-1148-4fdaa380dcfd
# ╟─4d443e46-7983-11eb-33e3-316d616ad44b
# ╠═9119e782-7984-11eb-25cb-471cefe4e7d5
# ╠═9642a69e-7988-11eb-069e-157feefb4eed
# ╠═9a2fe03e-7988-11eb-09ea-b709a8cbd77c
# ╠═d781854e-7988-11eb-385e-e3dedbaae06e
# ╠═a70ee822-7988-11eb-112e-5be62d6310c4
# ╠═ba3e633c-7988-11eb-2550-4b645fb2439d
# ╟─dde99fd6-7988-11eb-061a-15e9399613b6
# ╠═88d2e984-7984-11eb-3f4c-5115712ac1a2
# ╠═a7f67e52-7984-11eb-22b1-69c114b468fb
# ╟─bd098168-7986-11eb-110c-c5d17c258584
# ╠═cf8e1c18-7984-11eb-20f9-3ff941f8b457
# ╟─04a88244-7985-11eb-328d-6f3878193236
# ╟─aba0e4f6-7985-11eb-0713-172db02e7224
# ╟─d961ff42-7985-11eb-37d2-35af70b59310
# ╟─6636f208-7987-11eb-16a3-17e36a925c97
# ╟─1b6955c0-7986-11eb-180f-6fabd0e6cf36
# ╟─e7e6e6de-7987-11eb-2fc7-37fab2d3c05b
# ╟─eba873f4-7988-11eb-0c30-372c31f4cad8
# ╟─0697dcea-7989-11eb-0ead-8d17b8f4b6cd
# ╠═47ae6534-7984-11eb-3be1-3501475c06f9
# ╠═98f6a88c-7984-11eb-3882-d121f8488565
# ╟─f6fa502e-7cce-11eb-028e-c790c42e9bc2
# ╟─c2f6ff4e-7fef-11eb-0d7e-bd34ab2132af
# ╟─554abfc4-7ccf-11eb-25fb-df246eee94a0
# ╟─4a80d04c-7989-11eb-109b-4789a2c1480a
# ╟─87d63e9a-798a-11eb-393f-e59d4eb48bea
# ╟─98bc8d9a-798a-11eb-1fc3-0bd77bb3fe8c
# ╟─bc770cc6-798a-11eb-0333-c746ba22dc7f
# ╠═08692f0e-798b-11eb-2b02-9548f1a59426
# ╠═b0ea2b14-798c-11eb-0a6c-6f06820a0fa1
# ╟─15a6fa1e-798d-11eb-187e-b93de64da00b
# ╟─1aa85896-798d-11eb-15ff-7f64835af8c8
# ╟─24b05384-798d-11eb-346a-77856d43c699
# ╟─70ef989a-798d-11eb-1899-f51410a810eb
# ╠═970787c2-798d-11eb-12c7-bfd48299c383
# ╟─4b5fb2f4-798b-11eb-1af8-33890aaf0ed4
# ╠═242401ba-798b-11eb-31ed-07d5767cf6f2
# ╠═8071095e-798b-11eb-3e86-0b137ae8360f
# ╟─8842e8aa-798b-11eb-1746-f36c13e7c8ba
# ╟─5d79259e-798b-11eb-298a-b58a6288162e
# ╟─2f1b4990-798e-11eb-31fb-dbc8934ccdfb
# ╟─42b201a8-799d-11eb-00ce-27483756dce2
# ╟─673294b4-799d-11eb-07ac-9dbd30c54541
# ╠═6fce44da-799d-11eb-1be4-8b0d2097ce9a
# ╟─92d81aa8-799d-11eb-25e3-7b2d722ba0c2
# ╟─1082bef2-7989-11eb-0137-f36a82244d4f
# ╟─81db5570-79b2-11eb-2dcb-5b908ad62b3f
# ╟─0c4ac2ae-79b3-11eb-09be-b385fd8dd389
# ╠═5ae7f7e6-79b4-11eb-1609-470cf5246eb7
# ╠═95d0b9b0-79b4-11eb-2120-978ee85a75a3
# ╟─40841800-79b7-11eb-2667-25b9a152246d
# ╟─5aec1d14-79b7-11eb-0342-997f50dbe294
# ╟─b420b066-79b7-11eb-16f6-9fa333021ea9
# ╟─c8fa957e-79b7-11eb-355c-db21c14bccdc
# ╟─f84b2414-79b3-11eb-29d0-a3da4a9d41be
# ╟─5be6eb56-79b0-11eb-10a2-2bd5979f3a63
# ╟─1b53b1cc-79bb-11eb-16c2-a799ed6bafbd
# ╟─1056b4f2-79b3-11eb-09a6-9f7ffaf3b5e3
# ╟─b4645be0-7b49-11eb-1f50-9328aece496b
# ╟─07703e92-7b4a-11eb-0aa6-7180b261728c
# ╟─7e72da22-7cd2-11eb-3de0-93e1435b1987
# ╠═a0b559f6-79bb-11eb-0bb5-eb2b8e341b3c
# ╠═842b718a-79bb-11eb-15e3-fd52f3ed3be4
# ╠═1d2029b0-7b44-11eb-2df3-cf8c03c6c009
# ╠═369e4a72-7b44-11eb-2947-d198ea29598e
# ╠═549ce572-7b44-11eb-3ea9-3334c962ca90
# ╟─cc201ce0-7b44-11eb-37fc-358d02e2d8ed
# ╠═5a33dad6-7b44-11eb-36cb-034dc4851a01
# ╠═671c8c5a-7b44-11eb-0302-53101084c5a6
# ╠═7195a696-7b44-11eb-1073-91f2060f0fd5
# ╠═2e18943a-7ff1-11eb-1d35-7f2019c6c06e
# ╟─d4f04180-7b45-11eb-210a-75f898949c20
# ╟─afde9040-801a-11eb-3907-4f9f38b66893
# ╟─edf6aaae-7b45-11eb-3e5c-ff386f3ca855
# ╟─9ebeb81e-7cd2-11eb-3842-67d666de0d56
# ╟─176b9300-7b47-11eb-34f5-e9f401587fbb
# ╟─10813aba-7b46-11eb-15ae-9d34277557aa
# ╟─522289b6-7b47-11eb-2e36-f1fcaad71e46
# ╟─09b45022-7b49-11eb-3a09-25e2c5f15d5a
# ╠═dfaca6b8-7b47-11eb-3232-c395a968eab0
# ╠═ce181838-7b47-11eb-09fb-21ae2437e92f
# ╟─c3f0eb9c-7e5c-11eb-11ab-b39896644d0e
# ╟─1cf14bf8-7b4a-11eb-057c-cd6c391812d7
# ╟─8fe26172-7b4a-11eb-3fec-25a023ad0c21
# ╟─413e485c-7b4b-11eb-212e-81acd087a606
# ╠═1457d8d0-7b4b-11eb-3b05-fdea9ad03b02
# ╠═17a48a4c-7b4b-11eb-0935-c78e3fbe8875
# ╟─5f6eae0c-7b4b-11eb-1f34-891b85a352b3
# ╠═843168c4-7b4b-11eb-2d50-f50df5c611b5
# ╟─93ba3102-7b4b-11eb-1e40-4d22cef75d41
# ╟─5508697c-7e5d-11eb-18c8-0f27ee5cb6be
# ╟─754bdb0e-7e5d-11eb-24ea-a15f33748d49
# ╟─a39d3d72-7e5d-11eb-3fa8-b57f67497390
# ╟─8f75f92a-7ff1-11eb-07db-07a8b1ee8fd3
# ╟─c16cb486-7e5d-11eb-3931-616157e00d6d
# ╠═cb9c1654-7e5d-11eb-12a4-b1b814fb894a
# ╠═222b0534-7e5e-11eb-07ee-599111a5a85e
# ╠═3c7d66e8-7e5e-11eb-156b-bf16ab97429b
# ╠═1ba1e9a6-7e5e-11eb-26c0-c3259ca71dfa
# ╠═5802e8a2-7e5e-11eb-0e69-2d5da350e203
# ╠═72d949fa-7e5e-11eb-3661-8f14b46a3336
# ╟─a1d2283a-7e5e-11eb-0c22-a91232575d2e
# ╟─01c7f2d0-7e5d-11eb-0134-b527c31615cd
# ╟─eff8182c-7b4c-11eb-3608-077859fd9405
# ╟─3795afc2-7e5f-11eb-04dd-650633462059
# ╟─2f9ccf56-7b4d-11eb-0620-85cc73e63986
# ╟─5723264a-7b4d-11eb-3d80-d728e4e125ea
# ╟─393eab66-7e60-11eb-138b-b7b5d7578b9e
# ╟─50be8e64-7e60-11eb-0bde-770db17176c5
# ╟─8af58b82-7e60-11eb-345b-f5497a167d14
# ╟─5d793888-7b4d-11eb-1760-7d6e5a06d299
# ╟─45d7729e-7b4e-11eb-16b8-9fe470c4422a
# ╟─502ef336-7b4e-11eb-1984-4fb3fc2696a9
# ╟─8f23a268-7b52-11eb-14d3-a19f581d662e
# ╟─e62703f0-7e60-11eb-061b-05e2ebbae775
# ╟─d06cb7e4-7b45-11eb-0ecb-3fb5bd004188
# ╟─e7ee031a-7b44-11eb-1109-6d05b4ddca2e
# ╠═2806eb6a-7b45-11eb-227a-93964c77e93b
# ╠═89e42a9c-7b44-11eb-315d-7974de3531d9
# ╟─4c62ba78-7e62-11eb-080e-9f49bb8a60fb
# ╟─7953f36c-7e62-11eb-318c-c30dc4038629
# ╟─fc298e8e-7e62-11eb-1275-bfaaba9e6d13
# ╠═35da102c-7e61-11eb-3068-3f8961843a9f
# ╠═5d6f375c-7e61-11eb-28ab-4561d1cabdb9
# ╠═501021f2-7e61-11eb-3205-03823f96e0ba
# ╟─14bf236c-7e63-11eb-31c0-1b57f03ae4ff
# ╟─711c29be-7e65-11eb-0aad-8f82e60e2ece
# ╟─5ceecd68-7e63-11eb-20a2-cf9cdf7e1758
# ╟─fab8dd52-7e65-11eb-3456-8708ccd8cd4f
# ╠═f054b9aa-7e63-11eb-36dd-3bf1932030ac
# ╠═9c434dd6-7e63-11eb-3bba-736c1f794722
# ╠═7b34bb66-7e63-11eb-2e44-1b00e23874cb
# ╠═660f6682-7e63-11eb-396d-77ef6e6f1de2
# ╟─624943a8-7e66-11eb-12b0-6f9a43405d16
# ╟─8848bee0-7e6a-11eb-0cb3-63001bba774d
# ╟─91787884-7e6a-11eb-3821-a11457bb4bfa
# ╟─db73f682-7e6a-11eb-2890-adea119688ab
# ╠═d4666004-7e6a-11eb-1502-41d2bdf3df30
# ╠═fec9f55c-7e6a-11eb-337c-13273fff4450
# ╟─01e8107c-7e6b-11eb-1e0b-fd862028439a
# ╟─258692fe-7e6b-11eb-3401-e9e8d96483fe
# ╟─2b9af64e-7e6b-11eb-11be-714f382922dc
# ╠═0fc0b998-7e6b-11eb-34be-8be4f0e4fe08
# ╠═190b6e28-7e6b-11eb-069c-dbaa847caa19
# ╟─3e31d85c-7e6b-11eb-0cb8-6d171248d19d
# ╟─2cc216d4-7e6c-11eb-09c8-ddaa8b9b05b5
# ╟─3b9ea1c0-7e6c-11eb-1cbf-dfe3d0fbc9c1
# ╟─99172a5c-7e6c-11eb-39db-b76d5eb6403f
# ╟─a641d362-7e6c-11eb-3a34-e727f9afbaa7
# ╟─bb96eaa6-7e6c-11eb-3251-7dfe84995533
# ╟─e700935c-7e6c-11eb-3c5f-79ee52818f27
# ╟─ec198862-7e6c-11eb-2368-cb1d288a537a
# ╟─0121fce4-7e6d-11eb-00a3-77cef143d4e0
# ╟─272137ca-7e6d-11eb-1ccb-491acdf93e80
# ╟─5d67d74a-6b93-11eb-0df5-01899af4b4ab
# ╠═a53d2cd8-7e6e-11eb-1265-c3df598243b6
# ╟─add3881e-7e74-11eb-24d0-c3909978a450
# ╠═bbf9ee7a-7e6e-11eb-06d9-21df405b04a2
# ╠═ba25c50c-7e74-11eb-248f-95e03691af90
# ╟─c2e4c756-7e74-11eb-330c-db0b6ea51886
# ╠═c8c0d22e-7e6e-11eb-1e6d-37acc66bdadb
# ╠═ca2cce42-7e6e-11eb-1b9a-8958575b107e
# ╠═1024a3ec-7e70-11eb-3b96-391f330cfb32
# ╟─0e46c762-7e70-11eb-2650-135cff7436b8
# ╟─30f7ecf0-7e70-11eb-13d5-75ace1ff42cd
# ╟─371558a2-7e70-11eb-1661-430f14e76302
# ╠═907116e8-7e70-11eb-333f-290f0a7d80e8
# ╠═4766721c-7e71-11eb-0dd7-69fdac58cc30
# ╠═661feb52-7e71-11eb-1bf0-5b4a4f4e1c02
# ╠═719b971a-7e71-11eb-1a2d-e1b67112c405
# ╟─56e6b8f2-7e71-11eb-1203-75faf28f79e0
# ╟─b3419d86-7e71-11eb-1a63-87efbffc5d90
# ╠═c7d99398-7e71-11eb-2a80-e1ee835a675d
# ╠═cc00e586-7e71-11eb-05b9-c92bdd03dc26
# ╟─e3b93278-7e71-11eb-01fb-c1fb19c1af97
# ╠═8c02ca6a-7e71-11eb-0533-f56d312294ac
# ╟─d9d0162a-7e72-11eb-143c-a7a0fdaa02e4
# ╟─0970b240-7e73-11eb-2ee9-f391e64d584b
# ╟─629dae04-7e73-11eb-1075-69cffc768b82
# ╟─6e49f3b0-7ff3-11eb-00b5-2157914eb9ed
# ╟─da31badc-7e73-11eb-34a0-75d05bc41291
# ╟─0d7fbfbe-7e91-11eb-14eb-ed9b68f7b02e
# ╟─3fe26ca6-7e91-11eb-1382-dd95809c4187
# ╟─6be14370-7e70-11eb-36c1-97d6b5ab2e2a
# ╠═377f3020-7e74-11eb-22d1-3bf79ca13e7a
# ╠═649cd9e0-7e74-11eb-0997-4f7390f13936
# ╠═89350110-7e74-11eb-0cb1-6bb00c61e12b
# ╠═b39bb5ac-7e81-11eb-3da9-fb2f6ece1f8a
# ╠═b5d346f4-7e81-11eb-15a4-657b672f23a5
# ╠═bf5aaed6-7e81-11eb-3036-cbe315f2a59a
# ╟─fa35967e-7e81-11eb-05bc-ab321c3635c0
# ╠═2c5cdce8-7e82-11eb-3f31-eb1c12783dad
# ╠═7f9a4eb2-7e83-11eb-2667-7d3c451bfc09
# ╠═70a7f836-7e82-11eb-367d-4125cd359a72
# ╠═40d3b36a-7e82-11eb-140c-9757638394e7
# ╟─c080c712-7e83-11eb-1143-158e83f386e2
# ╟─2b426de4-7e84-11eb-277f-41bf0b1ea003
# ╟─330784ce-7e84-11eb-3513-374b23f1c836
# ╠═42f59704-7e84-11eb-1410-e3edb4a255ad
# ╠═22147248-7e85-11eb-1f56-17469bc96208
# ╠═4b9a1c3a-7e85-11eb-36b8-9b92628a80eb
# ╠═5bb290de-7e85-11eb-16a9-45daa0302e09
# ╠═6023d310-7e85-11eb-0a26-9d85c0398e18
# ╟─97d9cc9c-7e85-11eb-1e65-6d5d3fb245a7
# ╠═9f460e70-7e85-11eb-2c83-776eb9a8a401
# ╟─edb98cd0-7e85-11eb-0c7a-57d0fef4c7b3
# ╟─dee27332-7e8c-11eb-3896-2566e284cbcd
# ╠═fafcb580-7e93-11eb-139c-5574a72195f1
# ╟─1328a088-7e94-11eb-1893-8307ec0a0354
# ╟─272b2e92-7e86-11eb-1905-b3c59f94c4b8
# ╟─67d8500a-7e86-11eb-0cfb-916b7df799eb
# ╟─06c81c4a-7e87-11eb-2b8b-49c3e2185e5d
# ╟─699f856c-7e94-11eb-017f-b57a7bddc047
# ╟─7c0b35ac-7e94-11eb-19b8-3fa36803e662
# ╟─44f4d7ac-7e95-11eb-0bd1-d3eb98605f56
# ╠═f08b0aae-7e87-11eb-2044-f59130ffc70a
# ╠═7aa15fca-7e8c-11eb-3fd3-93d2e6eef233
# ╟─d98262b2-7e8c-11eb-3bf3-5fbf2fa3aaae
# ╟─2385d402-7e8d-11eb-07bf-17752bfcc9b6
# ╟─5b8f4b08-7e8d-11eb-0039-f5a7fa194ba1
# ╠═a615ea5c-7e8c-11eb-2346-a9a9b35a584a
# ╟─86d649e2-7e8d-11eb-3b53-47bb7f2c6a57
# ╟─ce4c4f58-7e8d-11eb-3332-ab9b8d8148b9
# ╠═c9b9ebe2-7e8d-11eb-0890-b3e23ceb3a83
# ╟─adc658c8-7e90-11eb-32d9-eb4c98405392
# ╟─469a8e50-7e93-11eb-21b0-d1df5dfb4b0e
# ╠═68e0762c-6c33-11eb-3d59-37d618e5c795


text:
module REoptPlots

using PlotlyJS

export
    plot_violin_compare_soc_efc,
    plot_electric_dispatch

include("dispatch.jl")
include("battery_degradation.jl")

end # module


text:
using Base.Iterators
using ArgParse

function readInput(filename::String)
    open(filename) do file
        data = read(file, String)

        rows = split(data, "\n")
        matrix = zeros(Int, length(rows), length(rows[1]))
        for (i, row) in enumerate(rows)
            row = split(row, "")
            row_int = zeros(Int, size(row))
            row_int[row .== "."] .= 0
            row_int[row .== ">"] .= 1
            row_int[row .== "v"] .= 2
            matrix[i, :] = row_int
        end
        return matrix
    end
end

function one_step(seafloor)
    seafloor_next = copy(seafloor)
    x, y = size(seafloor_next)
    moved_already = falses(size(seafloor_next))
    for i in 1:x
        for j in 1:y
            if seafloor_next[i, j] == 1 && !moved_already[i, j]
                j_next = j + 1
                if j_next > y
                    j_next = 1
                end
                if seafloor_next[i, j_next] == 0 && !moved_already[i, j_next]
                    seafloor_next[i, j_next] = 1
                    seafloor_next[i, j] = 0

                    moved_already[i, j_next] = true
                    moved_already[i, j] = true
                end
            end
        end
    end

    moved_already .= false
    for i in 1:x
        for j in 1:y
            if seafloor_next[i, j] == 2 && !moved_already[i, j]
                i_next = i + 1
                if i_next > x
                    i_next = 1
                end

                if seafloor_next[i_next, j] == 0 && !moved_already[i_next, j]
                    seafloor_next[i_next, j] = 2
                    seafloor_next[i, j] = 0

                    moved_already[i_next, j] = true
                    moved_already[i, j] = true
                end
            end
        end
    end
    return seafloor_next
end

function part1(seafloor)
    for i in countfrom(1)
        next = one_step(seafloor)
        all(seafloor .== next) && return i
        seafloor = next
    end
end

function parse_cli()
    s = ArgParseSettings()
    @add_arg_table s begin
        "--example"
            help = "Calculate example"
            action = :store_true
    end
    return parse_args(s)
end

args = parse_cli();
if args["example"] || true
    seafloor = readInput("example.txt");

    n_steps = part1(seafloor);
    @assert n_steps == 58;
end

seafloor = readInput("input.txt");

n_steps = part1(seafloor);
println(n_steps);


text:
function Project(Fun,CG,Param)
OrdPoly=CG.OrdPoly;
nz=Param.Grid.nz;
p=zeros(CG.NumG,nz,1);
fLoc=zeros(OrdPoly+1,OrdPoly+1);
X = Param.cache.X
JC = Param.cache.JC
for iz=1:nz
  for iF=1:Param.Grid.NumFaces
    for j=1:OrdPoly+1
      for i=1:OrdPoly+1
        x=0.5*(X[i,j,1,:,iF,iz]+X[i,j,2,:,iF,iz]);
        det=JC[i,j,iF,iz];
        fLoc[i,j]=Fun(x,Param)*det;
      end
    end
    p[CG.Faces[iF].Glob,iz]=p[CG.Faces[iF].Glob,iz]+reshape(fLoc,(OrdPoly+1)*(OrdPoly+1),1);
  end
end
p=p./CG.M;
return p
end


text:
using Distributions, Random
using Test, LinearAlgebra, PDMats

p  = 3
n1 = p + abs(10randn())
n2 = p + 1 + abs(10randn())
B  = rand( InverseWishart(p + 2, Matrix(1.0I, p, p)) )
b  = Matrix((n2 / n1) * I, 1, 1)

G = MatrixFDist(n1, n2, B)
F = MatrixFDist(n1, n2, b)
f = FDist(n1, n2)

@testset "MatrixFDist promotion during construction" begin
    R = MatrixFDist(Float32(n1), BigFloat(n2), Matrix{Float64}(B))
    @test partype(R) == BigFloat
end

@testset "MatrixFDist construction errors" begin
    @test_throws ArgumentError MatrixFDist(1, n2, B)
    @test_throws ArgumentError MatrixFDist(n1, 1, B)
    @test_throws ArgumentError mean(MatrixFDist(n1, p, B))
end

@testset "MatrixFDist params" begin
    df1, df2, PDB = params(G)
    @test df1 == n1
    @test df2 == n2
    @test PDB.mat == B
end

@testset "MatrixFDist dim" begin
    @test dim(G) == p
    @test dim(F) == 1
end

@testset "MatrixFDist size" begin
    @test size(G) == (p, p)
    @test size(F) == (1, 1)
end

@testset "MatrixFDist rank" begin
    @test rank(G) == p
    @test rank(F) == 1
    @test rank(G) == rank(rand(G))
    @test rank(F) == rank(rand(F))
end

@testset "MatrixFDist insupport" begin
    @test insupport(G, rand(G))
    @test insupport(F, rand(F))

    @test !insupport(G, rand(G) + rand(G) * im)
    @test !insupport(G, randn(p, p + 1))
    @test !insupport(G, randn(p, p))
end

@testset "MatrixFDist logpdf" begin
    Z = rand(F)
    z = Z[1, 1]
    @test logpdf(F, Z) ≈ logpdf(f, z)
end

@testset "MatrixFDist cov and var" begin
    @test vec(var(G)) ≈ diag(cov(G))
    @test reshape(cov(G), p, p, p, p) ≈ cov(G, Val(false))
    @test cov(F)[1] ≈ var(f)
    @test_throws ArgumentError cov(MatrixFDist(n1, p + 1, B))
end

@testset "MatrixFDist sample moments" begin
    @test isapprox(mean(rand(G, 100000)), mean(G) , atol = 0.1)
    @test isapprox(cov(hcat(vec.(rand(G, 100000))...)'), cov(G) , atol = 0.1)
end

@testset "MatrixFDist conversion" for elty in (Float32, Float64, BigFloat)
    Gel1 = convert(MatrixFDist{elty}, G)
    Gel2 = convert(MatrixFDist{elty}, G.W, n2, G.logc0)

    @test Gel1 isa MatrixFDist{elty, Wishart{elty, PDMat{elty,Array{elty,2}}}}
    @test Gel2 isa MatrixFDist{elty, Wishart{elty, PDMat{elty,Array{elty,2}}}}
    @test partype(Gel1) == elty
    @test partype(Gel2) == elty
end