Julia に翻訳--184 冪曲線回帰 y = a * x^b

#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。

Python による統計処理　
http://aoki2.si.gunma-u.ac.jp/Python/PythonDir.tar.gz

冪曲線回帰 y = a * x**b
http://aoki2.si.gunma-u.ac.jp/Python/power.pdf

ファイル名: power.jl　　関数名: power

翻訳するときに書いたメモ

==========#

using Statistics, Printf, Plots

function power(x, y; printout=true, graph=true)
    f(x, a, b) = a * x ^ b
    model = "y = a * x^b"
    n = length(x)
    n >= 2 || error("sample size must be greater than 1")
    (all(x .> 0) && all(y .> 0)) || error("all x[i] && y[i] must be positive")
    a, b = sreg(log.(x), log.(y))
    a = exp(a)
    predicted = f.(x, a, b)
    resid = y .- predicted
    println(model)
    println("a = $a, b = $b")
    if printout
        @printf("%12s %12s %12s %12s\n", "x", "y", "pred.", "resid.")
        for i =1:n
            @printf("%12.6f %12.6f %12.6f %12.6f\n", x[i], y[i], predicted[i], resid[i])
        end
    end
    if graph
        plt = scatter(x, y, grid=false, tick_direction=:out, title="\$$model\$",
                      color=:blue, markerstrokecolor=:blue, label="")
        min, max = extrema(x)
        w = 0.05(max - min)
        x2 = range(min - w, max + w, length=1000)
        y2 = f.(x2, a, b)
        plot!(x2, y2, linewidth=0.5, color=:red, label="")
        display(plt)
    end
    Dict(:a => a, :b => b, :predicted => predicted, :resid => resid,
        :x => x, :y => y, :model => model)
end

function sreg(x, y)
    a = cov(x, y) / var(x)
    mean(y) - a * mean(x), a
end

x = collect(1:10);
y = [2, 16, 54, 128, 250, 432, 686, 1024, 1458, 2000];
power(x, y)
# y = a * x^b
# a = 1.9999999999999984, b = 2.9999999999999996
#            x            y        pred.       resid.
#     1.000000     2.000000     2.000000     0.000000
#     2.000000    16.000000    16.000000     0.000000
#     3.000000    54.000000    54.000000     0.000000
#     4.000000   128.000000   128.000000     0.000000
#     5.000000   250.000000   250.000000     0.000000
#     6.000000   432.000000   432.000000     0.000000
#     7.000000   686.000000   686.000000     0.000000
#     8.000000  1024.000000  1024.000000     0.000000
#     9.000000  1458.000000  1458.000000     0.000000
#    10.000000  2000.000000  2000.000000     0.000000