Julia 入門06 -- R と比較して（関数の作り方）

以下は，R 入門のために書かれた

https://ai-trend.jp/programming/r-beginner/r-beginner/

を，Julia ではこんな風になるよということで書いたものである。

R の入門というページは多く存在する。

R に入門しようとしている人 ... Julia にしませんか？

---


# 関数の作り方

Julia での関数の書き方は 2 通りある。

一つは，他の言語にも似ている方法である。

```Julia
function 関数名(引数)
    処理内容
end
```

もう一つは，数学などで関数を定義するのと同じような書式である。

```Julia
関数名(引数）= 戻り値の定義
```

これは，一行でなくても `begin - end` で囲めば何行でもよい。

```Julia
関数名(引数）= begin
        処理内容
    end
```

## `function - end` による関数定義


コラッツの予想とは，「ある数が偶数なら2で割り，奇数なら3倍して1を足す。できた数にこの操作を繰り返していくと必ずいつかは結果が1になる。」というもものである。これを行う関数を作ってみる。

関数を作る前に，まず，これを計算するプログラムを書いてみる。[^1]

[^1]: 当たり前だが，プログラムが書けないと関数も書けない。





```julia
n = 5
while n != 1
    if n % 2 == 0
        n = n ÷ 2
    else
        n = 3n + 1 # 3*n + 1 と同じ
    end
    println(n)
end
```

    16
    8
    4
    2
    1


確かに，最後に 1 になった。
別の値についても確かめるときには，1 行目の n = 5 を変えるだけでよい。
しかし，例えば REPL で実行しているときには `while - end` の部分は再入力する訳ではないが，操作が面倒だ。
この部分は 1 度定義すれば，しょっちゅう変える必要はないものである。そこで，この部分を関数 `function - end` で記述する。
同時に，毎回変わる部分（今回は変数 n の内容）は引数として書く。[^2]

[^2]: 引数はカンマで区切って複数個使える。


```julia
function collatz(n)
    while n != 1
        if n % 2 == 0
            n = n ÷ 2
        else
            n = 3n + 1 # 3*n + 1 と同じ
        end
        println(n)
    end
end
```




    collatz (generic function with 1 method)



関数定義が終わったら，さまざまな `n` について，関数を適用してみる。

`関数名(引数の実際の値)` のように書けばよい。[^2]

[^2]: `引数の実際の値` は数値でなくてもよい。式（数値も式であるが）を書くことができる。


```julia
collatz(3)
```

    10
    5
    16
    8
    4
    2
    1



```julia
collatz(11)
```

    34
    17
    52
    26
    13
    40
    20
    10
    5
    16
    8
    4
    2
    1


### `関数名(引数) begin - end` の形式で定義する

本来このような処理内容が複数行にわたる関数を書くための形式ではないが，3 番目の書き方を紹介するつなぎとして説明しておく。

以下のように書いて，同じように実行される。


```julia
collatz2(n) = begin
    while n != 1
        if n % 2 == 0
            n = n ÷ 2
        else
            n = 3n + 1 # 3*n + 1 と同じ
        end
        println(n)
    end
end
```




    collatz2 (generic function with 1 method)




```julia
collatz2(3)
```

    10
    5
    16
    8
    4
    2
    1


## `関数名(引数) = 処理内容` の形式で定義する

この形式の関数定義は，処理内容が 1 行（多くても数行）で書ける場合に適している。

数学で `f(x) = a * x^2 + b * x + c` のように書くのと同じである。


```julia
quadraticequation(a, b, c, x) = a * x^2 + b * x + c
println(quadraticequation(1, 2, 3, 0.5))
println(quadraticequation(4.2, 2.7, 3.3, 3.14))
```

    4.25
    53.188320000000004


なお，この関数定義を最初の形式で書いても，何の問題もなく，Julia ではプログラムはコンパイルされてから実行されるので，どちらで記述しても同じものになる。


```julia
function quadraticequation2(a, b, c, x)
    a * x^2 + b * x + c
end
println(quadraticequation(1, 2, 3, 0.5))
println(quadraticequation(4.2, 2.7, 3.3, 3.14))
```

    4.25
    53.188320000000004


## 関数にすべきもの

Julia ではどのようなプログラムでも，たとえ 1 度しか使われないものであっても，関数にすべきである。

Julia は，プログラムをコンパイルして高速に実行できる訳であるが，関数にしないとそのメリットを享受できない。

途中結果を出力しないで，n = 5000000000000000 についての実行時間の計測をしてみると 0.000087 seconds であった。


```julia
@time begin
    n = 5000000000000000
    while n != 1
        if n % 2 == 0
            n = n ÷ 2
        else
            n = 3n + 1 # 3*n + 1 と同じ
        end
    end
end
```

      0.000087 seconds (249 allocations: 3.891 KiB)


それを関数化して実行してみる。


```julia
function collatz4(n)
    while n != 1
        if n % 2 == 0
            n = n ÷ 2
        else
            n = 3n + 1 # 3*n + 1 と同じ
        end
    end
end
```




    collatz4 (generic function with 1 method)




```julia
実行時間は 0.000005 seconds である。

なお，Julia の JIT コンパイルはプログラムの初回実行時に行われるので，初回の実行時間は多めになるが，二回目以降は JIT コンパイルの効果が現れる。
```


```julia
n = 5000000000000000
@time collatz(n)
```

      0.000005 seconds


関数化するメリットとしては，一連の処理をまとめるという点があげられる。

例えば，前回取り上げた

`GLM による多変量解析`
https://qiita.com/StrawBerryMoon/items/d565f9ae892dcd940d3e

では，分析結果から必要なものを取り出すためにその都度いろいろな関数を使った。一々指定するのが面倒なので，以下のような関数を作っておけば，煩わしさが減る。

以下をファイル（たとえば `"linearregression.js"`）に保存しておき，使用したいときに `include("linearregression.jl")` で呼び出せばよい。ファイルの 1 行目に `using` を書いておけば，使用するときに一々指定する必要がない。


```julia
using GLM, Distributions, StatsBase
function linearregression(form, dataframe; digits=5)
    formatpval(p) = p < 0.00001 ? "< 0.00001" : " = $(string(round(p, digits=digits)))"
    results = lm(form, dataframe);
    dofe = Int(dof_residual(results))
    Se = sum(residuals(results) .^ 2)
    rse = round(sqrt(Se / dofe), digits=digits)
    dofres = Int(dof_residual(results))
    R2 = round(r2(results), digits=digits)
    R2adj = round(adjr2(results), digits=digits)
    loglik = round(loglikelihood(results), digits=digits)
    AIC = round(aic(results), digits=digits)
    y = results.model.rr.y;
    meany = sum(y) / nobs(results)
    St = sum((y .- meany).^2)
    Sr = St - Se
    doft = length(y) - 1
    dofr = doft - dofe
    MSr = Sr / dofr
    MSe = Se / dofe
    F = MSr / MSe
    p = ccdf(FDist(dofr, dofe), F)
    println(results)
    println("Residual standard error: $rse on $dofe degrees of freedom")
    println("Multiple R-squared: $R2,  Adjusted R-squared: $R2adj")
    println("log likelihood: $loglik,  AIC: $AIC")
    print("F-statistic: $(round(F, digits=digits)) ")
    println("on $(Int(dofr)) and $(Int(dofe)) DF,  p-value: $(formatpval(p))")
end
```




    linearregression (generic function with 1 method)



実際に使用するときには最初にファイルを，`include` する。


```julia
include("/Users/foo/Desktop/Julia/linearregression.jl")
```




    linearregression (generic function with 1 method)


まずは，使用するデータをデータフレームに読み込む。

```julia
using CSV, DataFrames, StringEncodings
df = CSV.File(open(read, "/Users/foo/Downloads/sample-data.csv", enc"cp932")) |> DataFrame
```

そして，実際に分析してみる。

書く内容は `lm()` を直接呼ぶときの `lm(@formula(肺活量 ~ 血圧 + 体重), df)` と全く同じである。

これだけで，自分が必要とする情報を全部コンパクトに表示できる。

```julia
linearregression(@formula(肺活量 ~ 血圧 + 体重), df)
```

    StatsModels.DataFrameRegressionModel{LinearModel{LmResp{Vector{Float64}}, DensePredChol{Float64, LinearAlgebra.Cholesky{Float64, Matrix{Float64}}}}, Matrix{Float64}}
    
    Formula: 肺活量 ~ 1 + 血圧 + 体重
    
    Coefficients:
    ─────────────────────────────────────────────────────
                  Estimate  Std.Error   t value  Pr(>|t|)
    ─────────────────────────────────────────────────────
    (Intercept)  3533.07    830.819     4.25251    0.0002
    血圧            -14.2164    5.57863  -2.54837    0.0168
    体重             30.7853    5.09687   6.04005    <1e-05
    ─────────────────────────────────────────────────────
    Residual standard error: 368.26786 on 27 degrees of freedom
    Multiple R-squared: 0.6748,  Adjusted R-squared: 0.65071
    log likelihood: -218.25207,  AIC: 444.50413
    F-statistic: 28.01303 on 2 and 27 DF,  p-value: < 0.00001

このようなものを一度作っておけば，その後の時間の節約になる。