t分布線形回帰の実験と動画の作成¶

黒木玄

2017-10-06, 2018-03-25

t分布線形回帰の函数の定義¶

########### t-distribution linear regression library (Ver.0.1)

using PyPlot
using Distributions
using Optim

function hcatones(X)
    return hcat(ones(eltype(X),size(X,1)), X)
end

########### residual error

@inline function residualerror(a::AbstractArray{T,1}, X::AbstractArray{T,2}, Y::AbstractArray{T,1}) where T
    return Y .- X*a
end

@inline function residualerror(a, X, Y::AbstractArray{T,1}) where T
    residualerror(T.(a), reshape(T.(X),size(X,1),size(X,2)), Y)
end

function show_residualerror(a, X, Y)
    R = residualerror(a, X, Y)
    @printf("sum of squares of residual errors = %.2e", sum(R.^2))
    @printf(" -> sqrt = %.2e;", sqrt(sum(R.^2)))
    @printf(" mean = %.2e,", mean(R))
    @printf(" std = %.2e", std(R))
    println()
end

function sprint_residualerror(a)
    r = length(a)
    s = "y - (" * @sprintf("%.1e", a[1]) * " x[1]"
    for i in 2:r
        s *= " + " * @sprintf("%.1e", a[i]) * " x[$i]" 
    end
    return s * ")"
end

########## sd2hp

@inline function signeddistancetohyperplane(p::AbstractArray{T,1}, X::AbstractArray{T,2}, Y::AbstractArray{T,1}) where T
    return (X[:,1] + hcat(Y,X[:,2:end])*p)/norm(p)
end

@inline function signeddistancetohyperplane(p, X, Y::AbstractArray{T,1}) where T
    return signeddistancetohyperplane(T.(p), reshape(T.(X), size(X,1), size(X,2)), Y)
end

sd2hp(p,X,Y) = signeddistancetohyperplane(p,X,Y)

function p2a(p)
    a = similar(p)
    a[1] = -1.0/p[1]
    a[2:end] = -p[2:end]./p[1]
    return a
end

function a2p(a)
    p = similar(a)
    p[1] = -1.0/a[1]
    p[2:end] = a[2:end]./a[1]
    return p
end

function sprint_signeddistancetohyperplane(p)
    a = p2a(p)
    r = length(a)
    s = "[y - (" * @sprintf("%.1e", a[1]) * " x[1]"
    for i in 2:r
        s *= " + " * @sprintf("%.1e", a[i]) * " x[$i]" 
    end
    return s * @sprintf(")]/%.1e", norm(p)/p[1])
end

sprint_sd2hp(p) = sprint_signeddistancetohyperplane(p)

######### reg_normal

@inline rand_normal(sigma, n) = sigma * randn(n)
@inline pdf_normal(sigma, z; mu=0.0) = pdf(Normal(mu,sigma), z)
@inline logpdf_normal(sigma, z; mu=0.0) = logpdf(Normal(mu,sigma), z)

@inline function negloglik_normal(sigma, a, X, Y, errorfunc)
    if sigma > zero(sigma)
        return -sum(logpdf_normal.(sigma, errorfunc(a,X,Y)))
    else
        return Inf
    end
end

function solve_normal(X, Y;
        init = [1.0; zeros(size(X,2))],
        errorfunc = residualerror,
        solver = NelderMead(),
        options = Optim.Options())
    negloglik(x) = negloglik_normal(exp(x[1]), x[2:end], X, Y, errorfunc)
    init_tmp = init
    init_tmp[1] = log(init[1])
    return optimize(negloglik, init, solver, options)
end

function parse_normal(result)
    sigma = exp(result.minimizer[1])
    a = result.minimizer[2:end]
    aic = 2*result.minimum + 2*(1+length(a))
    return sigma, a, aic
end

reg_normal(X, Y; kwargs...) = parse_normal(solve_normal(X, Y; kwargs...))

function show_normal_residualerror(sigma, a, aic)
    println(@sprintf("AIC=%9.2e,  ", aic) * sprint_residualerror(a) * @sprintf(" ~ %.1e Normal(0,1)", sigma))
end

function show_normal_signeddistancetohyperplane(sigma, p, aic)
    println(@sprintf("AIC=%9.2e,  ", aic) * sprint_sd2hp(p) * @sprintf(" ~ %.1e Normal(0,1)", sigma))
end

show_normal_sd2hp(sigma, p, aic) = show_normal_signeddistancetohyperplane(sigma, p, aic)
show_normal(sigma, a, aic) = show_normal_residualerror(sigma, a, aic)

######### reg_t

@inline rand_t(rho, nu, n) = rho * rand(TDist(nu), n)
@inline pdf_t(rho, nu, z) = pdf(TDist(nu), z/rho)/rho
@inline logpdf_t(rho, nu, z) = logpdf(TDist(nu), z/rho) - log(rho)

@inline function negloglik_t(rho, nu, a, X, Y, errorfunc; rhomin = 0.0, numax = Inf)
    if rhomin < rho && zero(nu) < nu ≤ numax 
        return -sum(logpdf_t.(rho, nu, errorfunc(a,X,Y)))
    else
        return Inf
    end
end

### initial values are very important!
init_rho = 1.0
init_nu  = 1.0

function solve_t(X, Y;
        init = [init_rho; init_nu; zeros(size(X,2))],
        errorfunc = residualerror,
        rhomin = 1.0e-4,
        numax = 1.0e+4,
        solver = NelderMead(),
        options = Optim.Options())
    negloglik(x) = negloglik_t(exp(x[1]), exp(x[2]), x[3:end], X, Y, 
        errorfunc, rhomin=rhomin, numax=numax)
    init_tmp = init
    init_tmp[1] = log(init[1])
    init_tmp[2] = log(init[2])
    return optimize(negloglik, init, solver, options)
end

function parse_t(result)
    rho = exp(result.minimizer[1])
    nu = exp(result.minimizer[2])
    a = result.minimizer[3:end]
    aic = 2*result.minimum + 2*(2+length(a))
    return rho, nu, a, aic
end

reg_t(X, Y; kwargs...) = parse_t(solve_t(X, Y; kwargs...))

function show_t_residualerror(rho, nu, a, aic)
    println(@sprintf("AIC=%9.2e,  ", aic) * sprint_residualerror(a) * @sprintf(" ~ %.1e TDist(%.1e)", rho, nu))
end

function show_t_signeddistancetohyperplane(rho, nu, p, aic)
    println(@sprintf("AIC=%9.2e,  ", aic) * sprint_sd2hp(p) * @sprintf(" ~ %.1e TDist(%.1e)", rho, nu))
end

show_t_sd2hp(rho, nu, p, aic) = show_t_signeddistancetohyperplane(rho, nu, p, aic)
show_t(rho, nu, a, aic) = show_t_residualerror(rho, nu, a, aic)

show_t (generic function with 1 method)

negloglik(rho, nu) = -logpdf_t(rho, nu, 1.0)
figure()
rho = collect(logspace(-6,4,100))
nu  = collect(logspace(-4,6,100))
pcolormesh(rho, nu, log.(negloglik.(rho',nu)), cmap="rainbow")
colorbar()
xscale("log")
yscale("log")
xlabel("rho")
ylabel("nu")
grid("on", ls=":");

using PyCall
@pyimport matplotlib.animation as anim
function showgif(filename)
    open(filename) do f
        base64_video = base64encode(f)
        display("text/html", """<img src="data:image/gif;base64,$base64_video">""")
    end
end

showgif (generic function with 1 method)

######### trace_normal

function trace_normal(result, X, Y; paramkey = "centroid")
    trace  = result.trace
    niters = length(trace)

    value = Array{Float64}(niters+1)
    sigma = Array{Float64}(niters+1)
    a     = Array{Array{Float64}}(niters+1)

    sigma[1] = exp(result.initial_x[1])
    a[1]     = result.initial_x[2:end]
    value[1] = negloglik_normal(sigma[1], a[1], hcatones(X), Y, residualerror)
    
    sigma[2:end] = exp.([trace[i].metadata[paramkey][1] for i in 1:niters])
    a[2:end]     = [trace[i].metadata[paramkey][2:end] for i in 1:niters]
    value[2:end] = [trace[i].value for i in 1:niters]

    return niters, value, sigma, a
end
    
function plot1trace_normal(t, value, sigma, a, X, Y; 
        N = length(value), xmin=minimum(X), xmax=maximum(X), ymin=minimum(Y), ymax=maximum(Y))
    XX = hcatones(X)
    
    clf()
    
    ax1 = subplot2grid((11,21), (1,0), rowspan=10, colspan=10)
    scatter(X, Y, s=20, color="red", alpha=0.5)
    x = [xmin-1, xmax+1]
    plot(x, a[t+1][1] + a[t+1][2]*x)
    xlim(xmin, xmax)
    ylim(ymin, ymax)
    astr = @sprintf("%.2f", a[t+1][1])
    bstr = @sprintf("%.2f", a[t+1][2])
    title("($t)    y = $astr + $bstr x", fontsize=10)
    ax1[:grid]("on")

    ax2 = subplot2grid((11,21), (1,10), rowspan=10, colspan=10)
    plot(0:t, sigma[1:t+1])
    xlim(0,N)
    ylim(minimum(sigma[1:N]), maximum(sigma[1:N]))
    yscale("log")
    sigmastr = @sprintf("%.2e", sigma[t+1])
    title("sigma($t) = $sigmastr", fontsize=10)
    ax2[:grid]("on")
      
    #valuestr = @sprintf("%.2e", value[t+1])
    valuestr = "$(value[t+1])"
    suptitle("-LogLikelihood($t) = $valuestr", fontsize=11)
    tight_layout()
    
    plot()
end

######### trace_t
    
function trace_t(result, X, Y; paramkey = "centroid")
    trace  = result.trace
    niters = length(trace)

    value = Array{Float64}(niters+1)
    rho   = Array{Float64}(niters+1)
    nu    = Array{Float64}(niters+1)
    a     = Array{Array{Float64}}(niters+1)

    rho[1]   = exp(result.initial_x[1])
    nu[1]    = exp(result.initial_x[2])
    a[1]     = result.initial_x[3:end]
    value[1] = negloglik_t(rho[1], nu[1], a[1], hcatones(X), Y, residualerror)
    
    rho[2:end]   = exp.([trace[i].metadata[paramkey][1] for i in 1:niters])
    nu[2:end]    = exp.([trace[i].metadata[paramkey][2] for i in 1:niters])
    a[2:end]     = [trace[i].metadata[paramkey][3:end] for i in 1:niters]
    value[2:end] = [trace[i].value for i in 1:niters]

    return niters, value, rho, nu, a
end

function plot1trace_t(t, value, rho, nu, a, X, Y; 
        N = length(value), xmin=minimum(X), xmax=maximum(X), ymin=minimum(Y), ymax=maximum(Y))
    XX = hcatones(X)
    
    clf()
    
    ax1 = subplot2grid((11,21), (1,0), rowspan=10, colspan=10)
    scatter(X, Y, s=20, color="red", alpha=0.5)
    x = [xmin-1, xmax+1]
    plot(x, a[t+1][1] + a[t+1][2]*x)
    xlim(xmin, xmax)
    ylim(ymin, ymax)
    astr = @sprintf("%.2f", a[t+1][1])
    bstr = @sprintf("%.2f", a[t+1][2])
    title("($t)    y = $astr + $bstr x", fontsize=10)
    ax1[:grid]("on")

    ax2 = subplot2grid((11,21), (1,10), rowspan=5, colspan=10)
    plot(0:t, rho[1:t+1])
    xlim(0,N)
    ylim(minimum(rho[1:N]), maximum(rho[1:N]))
    yscale("log")
    rhostr = @sprintf("%.2e", rho[t+1])
    title("rho($t) = $rhostr", fontsize=10)
    ax2[:grid]("on")

    ax3 = subplot2grid((11,21), (6,10), rowspan=5, colspan=10)
    plot(0:t, nu[1:t+1])
    xlim(0,N)
    ylim(minimum(nu[1:N]), maximum(nu[1:N]))
    yscale("log")
    nustr = @sprintf("%.2e", nu[t+1])
    title("nu($t) = $nustr", fontsize=10)
    ax3[:grid]("on")
        
    #valuestr = @sprintf("%.2e", value[t+1])
    valuestr = "$(value[t+1])"
    suptitle("-LogLikelihood($t) = $valuestr", fontsize=11)
    tight_layout()
    
    plot()
end

plot1trace_t (generic function with 1 method)

テスト¶

X = collect(-5:5)
n = length(X)
XX = hcatones(X)

Y = 0.2*X
Y[1] -= -1.0
Y[n] -=  1.0

figure(figsize=(5,4))
grid("on", ls=":")
scatter(X, Y, s=15, color="red", alpha=0.5)

PyObject <matplotlib.collections.PathCollection object at 0x000000003423ED30>

normal regressions¶

sn¶

options = Optim.Options(store_trace=true, extended_trace=true)

Optim.Options{Float64,Void}(1.0e-32, 1.0e-32, 1.0e-8, 0, 0, 0, false, 0, 1000, true, false, true, 1, nothing, NaN)

@show sn = solve_normal(XX, Y, options=options)
@show sigma, b, aicn = parse_normal(sn)

sn = solve_normal(XX, Y, options=options) = Results of Optimization Algorithm
 * Algorithm: Nelder-Mead
 * Starting Point: [0.0,0.0,0.0]
 * Minimizer: [-1.1554232665897586,-1.2244629788883089e-5, ...]
 * Minimum: 2.898462e+00
 * Iterations: 95
 * Convergence: true
   *  √(Σ(yᵢ-ȳ)²)/n < 1.0e-08: true
   * Reached Maximum Number of Iterations: false
 * Objective Calls: 168
(sigma, b, aicn) = parse_normal(sn) = (0.3149242117922135, [-1.22446e-5, 0.109091], 11.796924901191222)

(0.3149242117922135, [-1.22446e-5, 0.109091], 11.796924901191222)

niters, value, sigma, b = trace_normal(sn, X, Y)
plot1frame(t) = plot1trace_normal(t, value, sigma, b, X, Y, xmin=-5.5, xmax=5.5, ymin=-1.0, ymax=1.0)
@show niters

N=niters
figure(figsize=(8,3.5))
plot1trace_normal(N, value, sigma, b, X, Y, N=N, xmin=-5.5, xmax=5.5, ymin=-1.0, ymax=1.0)
#plot1frame(N)

niters = 96

0-element Array{Any,1}

file = "sn_95.gif"
niters, value, sigma, b = trace_normal(sn, X, Y)

N = 95
frames = [zeros(Int64,10); 0:N; fill(N,10)]
interval = 100
plot1frame(t) = plot1trace_normal(t, value, sigma, b, X, Y, N=N, xmin=-5.5, xmax=5.5, ymin=-1.0, ymax=1.0)

fig = figure(figsize=(8,3.5))
@time myanim = anim.FuncAnimation(fig, plot1frame, frames=frames, interval=interval)
@time myanim[:save](file, writer="imagemagick")
sleep(0.2)
showgif(file)
clf()

  0.159650 seconds (53.02 k allocations: 2.795 MiB)
 43.671325 seconds (297.60 k allocations: 13.047 MiB, 0.02% gc time)

Table of Contents

t分布線形回帰の実験と動画の作成¶

t分布線形回帰の函数の定義¶

テスト¶

normal regressions¶

sn¶