ガウス過程と機械学習 第1章 lm.py version 01

20200504_1200.jl

# ==================================================================================================================================
# ISBN 978-4-06-152926-7
# ガウス過程と機械学習
# 
# [reference]
# http://chasen.org/~daiti-m/gpbook/python/lm.py
# http://chasen.org/~daiti-m/gpbook/data/nonlinear.dat
# https://nbviewer.jupyter.org/gist/genkuroki/a37894d5669ad13b4cd27da16096bfd2
# https://github.com/JuliaData/CSV.jl/issues/371                                <- `ignorerepeated = true`
# 
# [How to execute]
# Base.cd("?\\GitHub\\ISBN9784061529267\\julia\\chap01")
# Base.MainInclude.include("?\\GitHub\\ISBN9784061529267\\julia\\chap01\\lm_v03.jl")
# ==================================================================================================================================


module ISBN9784061529267

    import CSV    # v0.6.1
    import HTTP   # v0.8.14
    import Printf
    


    export main


    struct ObservedData
        size :: Core.Int64
        inpt :: Base.Vector{Core.Float64} # observed input
        otpt :: Base.Vector{Core.Float64} # observed output
    end


    struct RegressionDegree
        poly :: Core.Int64
        trig :: Core.Int64
    end


    struct RegressionWeight
        degree :: RegressionDegree
        value  :: Base.Vector{Core.Float64}
    end


    mutable struct EstimatedData
        size :: Core.Int64
        inpt :: Base.Vector{Core.Float64} # input to estimate
        otpt :: Base.Vector{Core.Float64} # estimated output
    end


    function read_data(file_path)

        # STEP.01
        # get the sample data from online
        # obj_CSV = CSV.read("../../downloaded/chap01/nonlinear.dat"; delim=' ', ignorerepeated = true, header=[:inpt, :otpt], footerskip=1) # access to local file
        obj_CSV = CSV.read(HTTP.get(file_path).body; delim=' ', ignorerepeated = true, header=[:inpt, :otpt], footerskip=1) # access to online file


        # STEP.02
        # convert the read data to the array on Julia
        data_size   = Base.size(obj_CSV.inpt, 1)
        buffer_inpt = Base.zeros(data_size)
        buffer_otpt = Base.zeros(data_size)

        for itr in 1:1:data_size
            buffer_inpt[itr] = obj_CSV.inpt[itr]
            buffer_otpt[itr] = obj_CSV.otpt[itr]
        end


        # STEP.END
        # return `ObservedData` instance which is generated from read data
        return ObservedData(data_size, buffer_inpt, buffer_otpt)

    end


    function compute_weight(degree::RegressionDegree, observed_data::ObservedData)

        # STEP.01
        # prepare the array to store the design matrix
        design_matrix = Base.ones(Core.Float64, observed_data.size, degree.poly + 2 * degree.trig + 1)


        # STEP.02
        # compute the polynominal part of the design matrix
        for itr_cl in 2:1:(degree.poly + 1)
            design_matrix[:, itr_cl] = design_matrix[:, itr_cl - 1] .* observed_data.inpt[:]
        end


        # STEP.03
        # compute the trigonometric function part of the design matrix
        for itr_cl in 1:1:degree.trig
            design_matrix[:, (degree.poly + 1) + (2 * itr_cl - 1)] = Base.cos.(itr_cl * observed_data.inpt)
            design_matrix[:, (degree.poly + 1) + (2 * itr_cl    )] = Base.sin.(itr_cl * observed_data.inpt)
        end



        # STEP.04
        # generate the transposed design matrix
        design_matrix_transposed = Base.transpose(design_matrix)


        # STEP.END
        # compute the weight of regression
        # return `RegressionWeight` instance
        return RegressionWeight( degree, Base.inv( design_matrix_transposed * design_matrix ) * design_matrix_transposed * observed_data.otpt )

    end


    function estimate_output(size::Core.Int64, observed_data::ObservedData, regression_weight::RegressionWeight)

        # STEP.01
        # compute the input
        buffer_extrema = Base.extrema(observed_data.inpt)
        buffer_input   = Base.range( buffer_extrema[1]; length = size, stop = buffer_extrema[2] )


        # STEP.02
        # compute the output (polynominal part)
        itr_degree    = regression_weight.degree.poly
        buffer_output = regression_weight.value[itr_degree] .+ regression_weight.value[itr_degree + 1] .* buffer_input

        while itr_degree > 1
            itr_degree    -= 1
            buffer_output  = buffer_input .* buffer_output .+ regression_weight.value[itr_degree]
        end


        # STEP.03
        # compute the output (trigonometric function part)
        for itr_degree in 1:1:regression_weight.degree.trig
            buffer_output += regression_weight.value[ (regression_weight.degree.poly + 1) + (2 * itr_degree - 1) ] * Base.cos.(itr_degree * buffer_input)
            buffer_output += regression_weight.value[ (regression_weight.degree.poly + 1) + (2 * itr_degree    ) ] * Base.sin.(itr_degree * buffer_input)
        end

        return EstimatedData(size, buffer_input, buffer_output)

    end


    function save_result(observed_data::ObservedData, regression_weight::RegressionWeight, estimated_data::EstimatedData)

        obj_txt_result = Base.open("result.csv", "w")


        Printf.@printf(obj_txt_result, "%s,%s\n", "degree", "weight")

        for itr in 1:1:(regression_weight.degree.poly + 2 * regression_weight.degree.trig + 1)
            Printf.@printf(obj_txt_result, "%d,%E\n", itr - 1, regression_weight.value[itr])
        end

        Printf.@printf(obj_txt_result, "\n")


        Printf.@printf(obj_txt_result, "%s,%s\n", "input", "output (observed)")

        for itr in 1:1:observed_data.size
            Printf.@printf(obj_txt_result, "%E,%E\n", observed_data.inpt[itr], observed_data.otpt[itr])
        end

        Printf.@printf(obj_txt_result, "\n")


        Printf.@printf(obj_txt_result, "%s,%s\n", "input", "output (estimated)")

        for itr in 1:1:estimated_data.size
            Printf.@printf(obj_txt_result, "%E,%E\n", estimated_data.inpt[itr], estimated_data.otpt[itr])
        end


        Base.close( obj_txt_result )

    end


    function main()

        # STEP.01
        # get the sample data from online
        obj_observed_data = read_data("http://chasen.org/~daiti-m/gpbook/data/nonlinear.dat")

        # STEP.02
        # compute the weight for regression
        obj_regression_weight = compute_weight( RegressionDegree(2, 1), obj_observed_data )

        # STEP.03
        # compute the estimated output using regression
        obj_estimated_data = estimate_output(101, obj_observed_data, obj_regression_weight)

        # STEP.04
        # save the computed result
        save_result(obj_observed_data, obj_regression_weight, obj_estimated_data)

    end

end


# ==================================================================================================================================


@time ISBN9784061529267.main()


# ==================================================================================================================================
# EOF
# ==================================================================================================================================

20200504_1200.plt

# ==================================================================================================================================
# ISBN 978-4-06-152926-7
# ガウス過程と機械学習
# 
# [reference]
# http://chasen.org/~daiti-m/gpbook/python/lm.py
# http://chasen.org/~daiti-m/gpbook/data/nonlinear.dat
# https://nbviewer.jupyter.org/gist/genkuroki/a37894d5669ad13b4cd27da16096bfd2
# https://github.com/JuliaData/CSV.jl/issues/371                                <- `ignorerepeated = true`
# ==================================================================================================================================

set datafile separator comma
set key on outside right center vertical Left reverse
set format x '%3.1f'
set format y '%3.1f'
set xlabel 'input'
set ylabel 'output'
set xzeroaxis linetype -1
set yzeroaxis linetype -1

plot 'result.csv' using 1:2 every ::1:1::1 with points title 'observed'   , \
     'result.csv' using 1:2 every ::1:2::2 with lines  title 'regression'

# ==================================================================================================================================
# EOF
# ==================================================================================================================================

result.csv