Listing

Stamm-Blatt Diagramm

# Copyright (C) 2017 Georg Hertkorn.
# Dieses Programm kann durch jedermann gemäß den Bestimmungen der
# Deutschen Freien Software Lizenz genutzt werden.
#
# Die Lizenz kann unter http://www.d-fsl.de abgerufen werden.
#
# This Program may be used by anyone in accordance with the terms of the
# German Free Software License
#
# The License may be obtained under http://www.d-fsl.org."


"""
    stemleaf(x)
    
    Generate a stem-leaf diagram of the data in `x`. 
"""
function stemleaf(x::Vector{T} where T<:Real)
      
    kmin = Int(div(minimum(x), 10))
    kmax = Int(div(maximum(x), 10))
    joffset = -kmin + 1

    res = [Array{Int,1}() for j in 1:kmax + joffset]
     
    for i in x
        j = Int(div(i, 10)) + joffset
        lsd = i % 10
        push!(res[j], lsd)
        #println("j = $j, lsd = $lsd")
    end

    j = kmin
    width = Int(floor(log(10, kmax))) + 1
    
    for r in res
        sj = lpad("$j", width, ' ')
        print(sj, " | ")
        # print(r)
        rord = sort(r)
        for i in rord
            print(i)
        end
        print('\n')
        j += 1
    end
end

Lineare Regression

# Copyright (C) 2017, 2018 Georg Hertkorn.
# Dieses Programm kann durch jedermann gemäß den Bestimmungen der
# Deutschen Freien Software Lizenz genutzt werden.
#
# Die Lizenz kann unter http://www.d-fsl.de abgerufen werden.
#
# This Program may be used by anyone in accordance with the terms of the
# German Free Software License
#
# The License may be obtained under http://www.d-fsl.org."

module DStLinReg

export residuals, vary, varab, varalpha, varbeta, cialpha, cibeta, nqvalues

import Distributions
Md = Distributions

using Statistics, LinearAlgebra


""" 
    residuals(alpha, beta, x, y)
    
Compute the residuals of `y` with respect to `alpha + beta * x`.

"""
function residuals(alpha::Real, beta::Real, x::T, y::T) where T<:AbstractArray
    return y .- (beta .* x .+ alpha)
end

""" 
    vary(alpha, beta, x, y)

Compute the variance of `y`.
"""
function vary(alpha::Real, beta::Real, x::T, y::T) where T<:AbstractArray
    n = length(x)
    res = residuals(alpha, beta, x, y)
    return dot(res, res) / (n - 2)
end


""" 
    varab(alpha, beta, x, y)

    Compute the variation of `alpha` and `beta`.
"""
function varab(alpha::Real, beta::Real, x::T, y::T) where T<:AbstractArray
    x2 = dot(x, x)
    mx = mean(x)
    n = length(x)
    varmod = vary(alpha, beta, x, y)
    vara = varmod * x2 / (n * (x2 - n * mx))
    varb = varmod / (x2 - n * mx)
    return (vara, varb)
end


2017-07-03
"""
    nqvalues(eps)
    
    Sort `eps` and return it together with quantiles for an NQ-plot.
"""
function nqvalues(eps::T) where T<:AbstractArray 
    epsord = sort(eps)
    dist = Md.Normal(0, 1)
    len = length(eps)
    q = [Md.quantile(dist, (i-0.5) / len) for i in 1:len]
    return (epsord, q)
end


"""
    varalpha(alpha, beta, x, y)

Compute the variance of `alpha`
"""
function varalpha(alpha::Real, beta::Real, x::T, y::T) where T<:AbstractArray
    n = length(x)
    x2 = dot(x, x)
    mx = mean(x)
    return vary(alpha, beta, x, y) * x2 / (n * (x2 - n * mx * mx))
end

"""
    varbeta(alpha, beta, x, y)

Compute the variance of `a`
"""
function varbeta(alpha::Real, beta::Real, x::T, y::T) where T<:AbstractArray
    n = length(x)
    x2 = dot(x, x)
    mx = mean(x)
    return vary(alpha, beta, x, y) / (x2 - n * mx * mx)
end

""" 
    cialpha(alpha, beta, x, y, perr)

Compute the confidence interval for `alpha`. 

`perr` is the probability of error.
"""
function cialpha(alpha::Real, beta::Real, x::T, y::T, perr::Real) where T<:AbstractArray
    n = length(x)
    dist = Md.TDist(n - 2)
    salpha = sqrt(varalpha(alpha, beta, x, y))
    hw = salpha * Md.quantile(dist, 1 - perr/2)              
    return (alpha - hw, alpha + hw)
end

""" 
    cibeta(beta, beta, x, y, perr)

Compute the confidence interval for `beta`. 

`perr` is the probability of error.
"""
function cibeta(alpha::Real, beta::Real, x::T, y::T, perr::Real) where T<:AbstractArray
    n = length(x)
    dist = Md.TDist(n - 2)
    sbeta = sqrt(varbeta(alpha, beta, x, y))
    hw = sbeta * Md.quantile(dist, 1 - perr/2)              
    return (beta - hw, beta + hw)
end


end