データ処理の前処理って...

https://archive.ics.uci.edu/ml/datasets/water+treatment+plant
にある，water-treatment.data を読んで，データの様子を見る，というところまでやってみるというのを，
Julia だったらこうやるというところで...

using CSV, DataFrames

ヘッダーがないので，header=false を指定する。
欠損値は "?" ではいっているので，missingstring="?" とする。
日付データが dateformat="D-dd/mm/yy" の形で入っている（読み込んだ後は，年は 0090 とかになる。本当は 1990 年なんだろうけど），

df = CSV.read("water-treatment.data", DataFrame, header=false, missingstring="?", dateformat="D-dd/mm/yy");

先頭を見てみる。

first(df, 5)

5 rows × 39 columns (omitted printing of 30 columns)

Column1 Column2 Column3 Column4 Column5 Column6 Column7 Column8 Column9
Date… Int64? Float64? Float64 Int64? Int64? Int64? Float64? Float64?
1 0090-03-01 44101 1.5 7.8 missing 407 166 66.3 4.5
2 0090-03-02 39024 3.0 7.7 missing 443 214 69.2 6.5
3 0090-03-04 32229 5.0 7.6 missing 528 186 69.9 3.4
4 0090-03-05 35023 3.5 7.9 205 588 192 65.6 4.5
5 0090-03-06 36924 1.5 8.0 242 496 176 64.8 4.0

行数，列数を確認する。

size(df) # (527, 39)

第1列目（日付）でソートする（df は直接書き換えられるので，代入しなくてよい）。

sort!(df, :Column1)

527 rows × 39 columns (omitted printing of 30 columns)

先頭を見てみる。

first(df, 5)

Column1 Column2 Column3 Column4 Column5 Column6 Column7 Column8 Column9
Date… Int64? Float64? Float64 Int64? Int64? Int64? Float64? Float64?
1 0090-01-01 41230 0.35 7.6 120 344 136 54.4 4.5
2 0090-01-02 37386 1.4 7.9 165 470 170 76.5 4.0
3 0090-01-03 34535 1.0 7.8 232 518 220 65.5 5.5
4 0090-01-04 32527 3.0 7.8 187 460 180 67.8 5.2
5 0090-01-07 27760 1.2 7.6 199 466 186 74.2 4.5

using Plots

データが取られている日付の間隔をみるのだが，間隔は diff() だけで求まる。それをヒストグラムで表す。

histogram(diff(df[:, 1]))

縦軸を対数表示にしたいということなので，yaxis=:log を指定する。

histogram(diff(df[:, 1]), yaxis=:log)

各項目（各列）で欠損値がどれくらいあるかということで，まずは，項目名と欠損値率を表示する。

hcat(names(df), describe(df)[:, :nmissing] ./ size(df,1))

39×2 Array{Any,2}:
 "Column1"   0.0
 "Column2"   0.0341556
 "Column3"   0.0056926
 "Column4"   0.0
 "Column5"   0.0436433
 ⋮           
 "Column38"  0.0151803
 "Column39"  0.0588235

棒グラフで見てもよいが。

bar(describe(df)[:, :nmissing] ./ size(df,1))

Julia のグラフで添え字などの表示

旨くいかないので，「.texで保存し、LaTeXでコンパイルしてpdfファイルを得ようとしている」という記事を読んで，やってみた。

一応できているのかな？

違うところは，Julia のバージョン？ Plots のバージョンは書いてなかったのでわからない。

using Plots

f(x) = sin(x)/x

g(x) = atan(x)

function plot_test()
    x = collect(-4:0.01:4) .* pi
    y1 = f.(x)
    y2 = g.(x)

    plot(x,y1,
        xaxis=(raw"$k_x$"),
        xlim=(-4pi, 4pi),
        xticks=(collect(-4:4) .* pi, ["$i \\pi" for i =-4:4]),
        yaxis=raw"$E_k$",
        label=raw"$f_k$",
        legend=:topleft,
        line=:solid
    )
    plot!(x,y2,
        label=raw"$\sigma_z$",
        legend=:topleft,
        line=:dash
    )

    annotate!(6, -0.8, Plots.text(raw"$\Sigma_\alpha=0.1 \textrm{eV/\AA}$", 16, :blue, :center))
end
plot_test()

3 次元空間の単位ベクトルを表す乱数（その4）

2014年8月5日に「3 次元空間の単位ベクトルを表す乱数」，2014年8月6日に「3 次元空間の単位ベクトルを表す乱数（その2）」という記事を2本書いた。

そのときは，ノルムが1より大きいベクトルを捨てていなかった。

n=3 のとき，ノルムが1より大きいベクトルを捨てないと，

のように，横軸が 1 〜2の間の頻度が比較的近い（横軸1あたりでカクンと変化がある）。

ノルムが1より大きいベクトルを捨ててシミュレーションすると，

のように，頭が丸くなる。

シミュレーションプログラムは以下の通り。

using Plots
plotly()

function sim(n=50000)
    xyz = (rand(n, 3) .- 0.5) .* 2
    r = vec(sum(xyz.^2, dims=2))
    xyz = xyz[r .<= 1, :]
    r = r[r .<= 1]
    factor = sqrt.(r)
    xyz ./ factor
end

function PearsonProblem(n = 4; dR = 0.1, trial = 100000)
    xyz = sim(2n*trial) # 歩留まり考慮の水増し発注
    N = size(xyz, 1)
    println(N÷n)
    dst = []
    for i = 1:n:N÷n
        dist = sqrt(sum(sum(xyz[i:i+n-1, :], dims=1) .^ 2))
        append!(dst, dist)
    end
    return dst
end

dst = PearsonProblem(3, trial=1000000);
histogram(dst)

 

球面上で一様に分布する点を表す x, y, z 座標値を生成する

「Julia でビュフォンの針」を書いたのには，速度比較だけではなく，もう一つ別の目的があったからなのだ。

ビュフォンの針のシミュレーションで，原点からのベクトルを考えるときに，方向がランダムなベクトルを作るときに，「長さが1以上のベクトルは除く」という操作をした後のベクトルだけを使わないと方向がランダムなベクトルは作れないということだった。

3次元空間でも同じで，球面に一様分布する点を表す座標を求めるときにも同じなのだ。-1 <= x, y, z <= 1 の一様乱数 x, y, z を発生させ，原点を起点，x, y, z を終点とするベクトルの内，ノルムが 1 より大きいは捨てる。ノルムが 1 以下のベクトルはノルムが 1 になるように伸ばすそうすれば，ベクトルの終点は球面に一様分布する。

以下のプログラムでできる x, y, z を，3次元空間にプロットできる。

グラフをぐりぐり動かすと点が一様分布している（らしい）ことが確認できる。

using Plots
plotly()

n=10000
xyz = (rand(n, 3) .- 0.5) .* 2
r = vec(sum(xyz.^2, dims=2))
xyz = xyz[r .<= 1, :] #@@@
r = r[r .<= 1] #@@@
factor = sqrt.(r)
xyz = xyz ./ factor
p = scatter(xyz[:, 1], xyz[:, 2], xyz[:, 3], markersize=0.1)

#@@@ を付けた 2 行を省いて，わかりやすくするために点の数を 5 倍にして図を描いてみると以下のようになる。濃淡があるということは，球面上で点が一様分布していないということである。

Julia でビュフォンの針

以前に掲載したが，π を使わないシミュレーション

n = 100000000 で，

Python 35.806 秒
R      14.543 秒
Julia   5.836 秒

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

★★　Python

import numpy as np
def buffon(n = 100000000): # 取りあえずの試行回数（実際はこれより少なくなる）
    d = 10 # 平行線の間隔
    l = d / 2 # 針の長さ
    # 以下 5 行で，0 <= theta < π/2 の一様乱数をもとめ，その sine.theta を生成する
    xy = np.random.rand(2, n) # 一辺 1 の正方形内の一様乱数による座標点
    r = (xy**2).sum(axis=0) # 原点からの距離
    xy = xy[:, r <= 1] # 半径 1 の円の内部にある点に限る
    r = r[r <= 1] # 半径 1 の円の内部にある点に限る
    sine_theta = xy[1,] / np.sqrt(r) # 三角関数の定義により
    n = len(r) # sine.theta の個数
    y = np.random.rand(n)*l # 同じ個数の[0, l) の一様乱数（針の中心から，平行線までの距離）
    return 2 * l / d / np.mean(y <= l / 2 * sine_theta) # 条件を満たすものの逆数が π の推定値

from time import time
s = time(); print(buffon(100000000), time() - s) # 3.1409810174102484 35.805713176727295

★★　R

buffon = function(n = 1e+06) { # 取りあえずの試行回数（実際はこれより少なくなる）
 d = 10 # 平行線の間隔
 l = d/2 # 針の長さ
 # 以下 5 行で，0 <= theta < π/2 の一様乱数をもとめ，その sine.theta を生成する
 xy = matrix(runif(2 * n, 0, 1), 2) # 一辺 1 の正方形内の一様乱数による座標点
 r = colSums(xy^2) # 原点からの距離
 xy = xy[, r <= 1] # 半径 1 の円の内部にある点に限る
 r = r[r <= 1] # 半径 1 の円の内部にある点に限る
 sine.theta = xy[2, ]/sqrt(r) # 三角関数の定義により
 n = length(r) # sine.theta の個数
 y = runif(n, 0, l) # 同じ個数の[0, l) の一様乱数（針の中心から，平行線までの距離）
 2 * l/d/mean(y <= l/2 * sine.theta) # 条件を満たすものの逆数が π の推定値
}

system.time(print(buffon(100000000)))
# [1] 3.141178
#    ユーザ   システム       経過  
#     14.543      2.852     17.349 

★★　Julia

function buffon(n = 1000000) # 取りあえずの試行回数（実際はこれより少なくなる）
    # 以下 5 行で，0 <= theta < π/2 の一様乱数をもとめ，その sinθ を生成する
    xy = rand(2, n);  # 一辺 1 の正方形内の一様乱数による座標点
    r = vec(sum(xy .^ 2, dims=1)); # 原点からの距離
    xy = xy[:, r .<= 1]; # 半径 1 の円の内部にある点に限る
    r = r[r .<= 1]; # 半径 1 の円の内部にある点に限る
    sinθ = xy[2, :] ./ sqrt.(r); # 三角関数の定義により
    n = length(r) # sinθ の個数
    n / sum(rand(n) .<= (0.5 .* sinθ)) # 条件を満たすものの逆数が π の推定値
end

@time buffon(100000000) #   5.836088 seconds (30 allocations: 7.269 GiB, 6.80% gc time)
# 3.141511691253939

ピアソンの積率相関係数（計算，検定，信頼区間，パワーアナリシス）

using Statistics # cor()
using Distributions # t()
using Printf # @printf()

ピアソンの積率相関係数の計算，母相関係数=0 の検定，信頼区間

function cortest(x, y; conflevel=0.95)
    r = cor(x, y)
    n = length(x)
    t = abs(r) *sqrt(n - 2) / sqrt(1 - r^2)
    p = ccdf(TDist(n - 2), t) * 2
    ξr = atanh(r)
    Zα2 = cquantile(Normal(), (1 - conflevel) / 2)
    Z = ξr .+ [-1, 1] .* Zα2 ./ sqrt(n - 3)
    conf = tanh.(Z)
    @printf("n = %d,  r = %.3f\n", n, r)
    @printf("t = %.5g,  df = %d,  p value = %.5g\n", t, n - 2, p)
    @printf("%d%% CI = [%.3f, %.3f]\n", conflevel*100, conf[1], conf[2])
end

有意水準，検出力を指定して必要なサンプルサイズを求める（両側検定のとき）

function poweranalysis(r; alpah=0.05, power=0.8)
    Za = cquantile(Normal(), alpah/2)
    Zb = quantile(Normal(), power)
    z = atanh(r)
    n = ((Za + Zb) / z) ^ 2 + 3
    @printf("n = %d (%.3f)\n", ceil(Int, n), n)
end

chukan = [64,40,71,33,30,71,92,23,41,55,93,74];
kimatu = [55,52,76,24,48,87,100,30,35,67,86,81];

cortest(chukan, kimatu)

n = 12,  r = 0.919
t = 7.3927,  df = 10,  p value = 2.3348e-05
95% CI = [0.731, 0.978]

poweranalysis(0.919)

n = 7 (6.134)

プログラミング書法 R, Python, Julia （その 2）

問題: 文字列内の各アルファベットを10個後のものに変換せよ

小文字に限定されているのかな？

move_ten("abcdefghijklmnopqrstuvwxyz")
の解が
"klmnopqrstuvwxyzabcdefghij"

Python プログラム

解答例

def move_ten(st):
    result = ""
    for i in st:
        ordI = ord(i)
        if ordI >= 113:
            result += chr(ordI - 16)
        else:
            result += chr(ordI + 10)
    return result

def move_ten(st):
    return ''.join(chr(ord(i)+10) if i<'q' 
        else chr(ord(i)-16) for i in st)

R プログラム

move_ten=function(st)chartr("a-z","k-za-j",st)

Julia プログラム

function move_ten(st)
    join(i < 'q' ? Char(i+10) : Char(i-16) for i in st)
end

move_ten(st)=join(Char(i+(i<'q' ? 10 : -16)) for i in st)

move_ten(st)=join(Char(i+10-(i>'p')*26) for i in st)

R はそれ用の関数 chartr があるので，かなわないなあ。最短。

追記：

英大文字・小文字しか含まないなら，以下のようにすれば大文字にも対応

move_ten(st)=join(Char(i+10-((95&Int(i))>80)*26) for i in st)

move_ten("abcdefghijklmnopqrstuvwxyz") # "klmnopqrstuvwxyzabcdefghij"
move_ten("ABCDEFGHIJKLMNOPQRSTUVWXYZ") # "KLMNOPQRSTUVWXYZABCDEFGHIJ"

 

プログラミング書法 R, Python, Julia

ユークリッドの互除法

R の場合

gcd <- function(a, b) {
    if (!(a%%1 == 0 & b%%1 == 0 & a > 0 & b > 0)) { # && でないとよくない
        cat("入力が自然数じゃないのでやり直し")
    } else if (a < b) { # これを特別扱いする必要はない
        w <- a
        a <- b
        b <- w
        
        while (b != 0) { # 入れ替えたのだから重複して書く必要はない
            r <- a%%b
            a <- b
            b <- r
        }
        return(a)
    } else {
        while (b != 0) {
            r <- a%%b
            a <- b
            b <- r
        }
        return(a)
    }
}

gcd(2, 0) # 入力が自然数じゃないのでやり直し
gcd(0, 22) # 入力が自然数じゃないのでやり直し
gcd(2.3, 4) # 入力が自然数じゃないのでやり直し
gcd(2, 4.1) # 入力が自然数じゃないのでやり直し
gcd(50856, 96007) # 163
gcd(96007, 50856) # 163

gcd2 <- function(a, b) {
    if (a%%1 != 0 || b%%1 != 0 || a <= 0 || b <= 0) {
        return("入力が自然数じゃないのでやり直し")
    }
    while (b != 0) {
        r <- a%%b
        a <- b
        b <- r
    }
    return(a)
}

gcd2(2, 0) # 入力が自然数じゃないのでやり直し
gcd2(0, 22) # 入力が自然数じゃないのでやり直し
gcd2(2.3, 4) # 入力が自然数じゃないのでやり直し
gcd2(2, 4.1) # 入力が自然数じゃないのでやり直し
gcd2(50856, 96007) # 163
gcd2(96007, 50856) # 163

================================================================

Python の場合

def gcd(a, b):
  if not (a % 1 == 0 and b % 1 == 0 and a > 0 and b > 0): # not (x == y and ...) は書き換えできる
    print('入力が自然数じゃないのでやり直し')
  elif a < b: # 場合分けは不要
    w = a
    a = b
    b = w

    while not b == 0:
      r = a % b
      a = b
      b = r
    else:
      return(a) # Python の場合は return は関数ではないので (a) としなくてよい
  else:
    while not b == 0: # not b == 0 は b != 0 または単に b
      r = a % b
      a = b
      b = r
    else:
      return(a)

gcd(2, 0) # 入力が自然数じゃないのでやり直し
gcd(0, 22) # 入力が自然数じゃないのでやり直し
gcd(2.3, 4) # 入力が自然数じゃないのでやり直し
gcd(2, 4.1) # 入力が自然数じゃないのでやり直し
gcd(50856, 96007) # 163
gcd(96007, 50856) # 163

def gcd2(a, b):
  if a % 1 != 0 or b % 1 != 0 or a <= 0 or b <= 0:
    return '入力が自然数じゃないのでやり直し'
  while b:
    a, b = b, a % b # R にはない記述法
  return a

gcd2(2, 0) # 入力が自然数じゃないのでやり直し
gcd2(0, 22) # 入力が自然数じゃないのでやり直し
gcd2(2.3, 4) # 入力が自然数じゃないのでやり直し
gcd2(2, 4.1) # 入力が自然数じゃないのでやり直し
gcd2(50856, 96007) # 163
gcd2(96007, 50856) # 163

================================================================

Julia  の場合

function gcd2(a, b)
  (a%1 == 0 && b%1 == 0 && fa > 0 && b > 0) ||
    return "入力が自然数じゃないのでやり直し" # R, Python にはない記述法
  while b
    a, b = b, a % b # R にはない記述法
  end
  return a
end

gcd2(2, 0) # 入力が自然数じゃないのでやり直し
gcd2(0, 22) # 入力が自然数じゃないのでやり直し
gcd2(2.3, 4) # 入力が自然数じゃないのでやり直し
gcd2(2, 4.1) # 入力が自然数じゃないのでやり直し
gcd2(50856, 96007) # 163
gcd2(96007, 50856) # 163

 

機械学習にもの申す

機械学習における誤認識

> 数値の大きさとばらつきに大きな差があるため、同じ基準では評価できません。
> そのため、標準化を行い、スケールを揃える必要があります。

そのようなことは必要なし　後に示す
「特徴量」，「正解」という言葉を使うが，独立変数（説明変数），従属変数（目的変数）と呼ぼう。
少なくとも，「正解」というのは変だ。
重回帰分析をする際にさえもデータを標準化する（平均値=0，標準偏差=1）。
標準化データで「モデルを作成」し，「モデルを訓練」する。
言葉も変だが，重回帰分析プログラムはそれ自体が内部でデータを標準化して偏回帰係数と必要なら定数項を計算する。
前もって標準化したデータを使って重回帰すると標準化を 2 回することになる。
標準化は何回行ってもかまわないが，無駄であるだけでなく，そのあとのあらゆる計算処理を複雑にする。

Python で以下のように行われている。

#住宅価格以外のデータを特徴量Xに設定
X = df.iloc[:, 0:13].values
#住宅価格を正解yに設定
y = df['MEDV'].values
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.2, random_state=0)
#特徴量の標準化
sc = StandardScaler()
#訓練データで標準化モデルを作成し変換
X_train_std = sc.fit_transform(X_train)
#作成した標準化モデルでテストデータを変換
X_test_std = sc.transform(X_test)
#標準化された訓練データ
X_train_std[0]
#モデルの作成
model = LinearRegression()
#モデルの訓練
model.fit(X_train_std,y_train)
print('傾き： ' , model.coef_)
print('切片： ' ,  model.intercept_)
傾き：  [-0.97082019  1.05714873  0.03831099  0.59450642 -1.8551476   2.57321942
 -0.08761547 -2.88094259  2.11224542 -1.87533131 -2.29276735  0.71817947
 -3.59245482]
切片：  22.611881188118836

------------------------------------------------

統計学の観点からは，従属変数と独立変数を定めて，解を求める。それだけだ。

using DataFrames, CSV
df_train = CSV.read("train.csv", DataFrame);
df_test = CSV.read("test.csv", DataFrame);
names(df_train); # "CRIM", "ZN", "INDUS", "CHAS", "NOX", "RM", "AGE"
                 # "DIS", "RAD", "TAX", "PTRATIO", "B", "LSTAT", "MEDV"
using GLM
a = lm(@formula(MEDV ~ CRIM+ZN+INDUS+CHAS+NOX+RM+AGE+DIS+RAD+TAX+PTRATIO+B+LSTAT), df_train);
a

StatsModels.TableRegressionModel{LinearModel{GLM.LmResp{Array{Float64,1}},GLM.DensePredChol{Float64,LinearAlgebra.Cholesky{Float64,Array{Float64,2}}}},Array{Float64,2}}

MEDV ~ 1 + CRIM + ZN + INDUS + CHAS + NOX + RM + AGE + DIS + RAD + TAX + PTRATIO + B + LSTAT

Coefficients:
─────────────────────────────────────────────────────────────────────────────────
                    Coef.  Std. Error      t  Pr(>|t|)     Lower 95%    Upper 95%
─────────────────────────────────────────────────────────────────────────────────
(Intercept)   38.0917      5.52243      6.90    <1e-10   27.2342      48.9491
CRIM          -0.119443    0.0366674   -3.26    0.0012   -0.191534    -0.0473529
ZN             0.04478     0.0144348    3.10    0.0021    0.0164002    0.0731597
INDUS          0.00548526  0.0634058    0.09    0.9311   -0.119175     0.130145
CHAS           2.3408      0.902172     2.59    0.0098    0.567075     4.11453
NOX          -16.1236      4.21177     -3.83    0.0002  -24.4042      -7.84299
RM             3.70871     0.45754      8.11    <1e-14    2.80916      4.60826
AGE           -0.00312108  0.0143289   -0.22    0.8277   -0.0312926    0.0250504
DIS           -1.3864      0.213951    -6.48    <1e-9    -1.80704     -0.965756
RAD            0.244178    0.070156     3.48    0.0006    0.106247     0.38211
TAX           -0.0109896   0.00389829  -2.82    0.0051   -0.0186539   -0.00332534
PTRATIO       -1.04592     0.136975    -7.64    <1e-12   -1.31522     -0.77662
B              0.00811011  0.00295064   2.75    0.0063    0.00230896   0.0139113
LSTAT         -0.492793    0.0542393   -9.09    <1e-17   -0.599431    -0.386155
─────────────────────────────────────────────────────────────────────────────────

機械学習で得られた定数項と偏回帰係数が違うじゃないか。

傾き：  [-0.97082019  1.05714873  0.03831099  0.59450642 -1.8551476   2.57321942
 -0.08761547 -2.88094259  2.11224542 -1.87533131 -2.29276735  0.71817947
 -3.59245482]
切片：  22.611881188118836

まず，上に示した偏回帰係数 Coef. は，標準化しないデータ（元データ）に掛ける係数である。

予測値は

predict(a);

で得られる。その最初の方を見てみると

32.55692655238963
21.927094777364772
27.543825734416142
23.603188290149532
 6.5719096200611755
 :

である。

機械学習では，

#訓練データ、テストデータの住宅価格を予測
y_train_pred = model.predict(X_train_std)
y_train_pred
array([32.55692655, 21.92709478, 27.54382573, 23.60318829,  6.57190962,

ということであるから，予測値は同じである。

機械学習の結果示されたのは，統計学では「標準化偏回帰係数」と呼ばれるもので，
「独立変数を標準化して重回帰分析を行ったときの偏回帰係数」というものである。

上に示された偏回帰係数にそれぞれの変数の標準偏差を掛けたものに等しい。

coefficient = coef(a)

14-element Array{Float64,1}:
  38.09169492628502
  -0.11944344700245542
   0.044779951066505716
   0.005485261681822424
   2.3408036062417357
 -16.123604315423176
   3.708709012220785
  -0.0031210817807432118
  -1.38639737027842
   0.24417832698878614
  -0.010989636563082644
  -1.0459211887458149
   0.008110106932704319
  -0.4927927245046289

coefficient[1] は定数項である。
coefficient[2:14] が独立変数に対する偏回帰係数である。
df_train[:. 1:13]が独立変数データ行列である。
よって，

using Statistics
coefficient[2:14] .* std(Matrix(df_train[:, 1:13]), corrected=false, dims=1)'

13×1 Array{Float64,2}:
-0.970820190461867 1.0571487264361754 0.038310993617902085 0.5945064227496165 -1.8551475975879652 2.573219418816067 -0.08761546983261663 -2.880942593098789 2.1122454191874787 -1.8753313056828191 -2.2927673485303135 0.7181794734592682 -3.5924548228794757

となり，機械学習のときの結果と同じであることがわかる。

傾き：  [-0.97082019  1.05714873  0.03831099  0.59450642 -1.8551476   2.57321942
 -0.08761547 -2.88094259  2.11224542 -1.87533131 -2.29276735  0.71817947
 -3.59245482]

標準化偏回帰係数が必要なときには，このように追加の計算をするだけである。

次は，機械学習のときの
切片：  22.611881188118836
である。
これは簡単で，標準化しないままの従属変数の平均値である。

print(mean(df_train[:, 14]))
22.611881188118815

あれこれ書いてきたが，機械学習派は「予測値が同じならドッチでもよいじゃないか」というであろう。
しかし，重回帰分析の一つの目的である「予測」という観点からは，元のデータに掛ける偏回帰係数を示す方がよいだろう。

そうしないと，test データセットにおける予測値を計算するときに面倒な手続きが必要になる。
まず，test データセットの標準化をしなければならないが，その時に使う平均値と標準偏差は train データセットのものを使わないといけない。

#作成した標準化モデルでテストデータを変換
X_test_std = sc.transform(X_test)
y_test_pred = model.predict(X_test_std)
y_test_pred
array([24.88963777, 23.72141085, 29.36499868, 12.12238621, 21.44382254,
       19.2834443 , 20.49647539, 21.36099298, 18.8967118 , 19.9280658 ,
        5.12703513, 16.3867396 , 17.07776485,  5.59375659, 39.99636726,

しかし，下手の考えで標準化してからあれこれさらに余計なことなどしなくても，元の test データセットそのままで，以下のようにするだけで予測値が求まる。

predict(a, df_test)

24.889637772755968
23.72141085316374
29.364998677355622
12.122386208110687
21.44382254072625
19.28344430480102
20.496475391456304
21.360992984142296
 :

というように，以上のように，機械学習は無駄なことばかりやっている。
それ以上に問題なのは，「必要なことをやっていない」ということである。
それは，分析結果の吟味である。
機械学習は，予測値がどれくらい正確かダケを見ている。

> MSEとはモデルの評価指標の一つで、この値が低いほどモデルの性能が良いと言えます。
> 算出方法は、予測値と正解の平均二乗和誤差です。

#MSEの計算
MSE_train = np.mean((y_train_pred - y_train) ** 2)
MSE_test = np.mean((y_test_pred - y_test) ** 2)
print('MSE train: ', MSE_train)
print('MSE test:', MSE_test)

MSE train:  19.32647020358573 # mean(residuals(a).^2)
MSE test: 33.44897999767649 # mean((predict(a, df_test) .- df_test[:, 14]) .^ 2)

何と比べて MSE が小さいのか？
test と比べて小さい（逆に言えば test の MSE が大きい）
それはある意味仕方ない。サンプルサイズが train:test ≒ 4:1 なのだから。

他のモデルと比べて MSE が小さいといわなければあまり意味がないのではないか。

重回帰分析結果の評価でよく使われるのは，決定係数 R2 と自由度調整済み決定係数 adj. R2 である。
これは，従属変数の分散が，独立変数によってどの程度説明できるかを表すもので，0 〜 1 の値をとり，
1 に近いほど予測がうまくいっているということになる。

r2(a) # 0.7730135569264233
adjr2(a) # 0.765447342157304

であるから，まずますではあるが，あまりよくはない。

重回帰分析の観点からは，上の方に示した GLM での解析結果 Coefficients で，Pr(>|t|) の列は，
その「変数に対する偏回帰係数が0である」という帰無仮説に対する P 値である。この値がたとえば
0.05 より大きければ，「有意水準5%で帰無仮説は棄却される」。
逆に言えば，P > 0.05 なら，「その変数の偏回帰係数は 0 でないとはいえない」のだ。

INDUS では 0.9311, AGE では 0.8277 なので，この 2 変数は予測にはあまり役立っていないということである。

b = lm(@formula(MEDV ~ CRIM+ZN+CHAS+NOX+RM+DIS+RAD+TAX+PTRATIO+B+LSTAT), df_train);

r2(b) # 0.7729811407453248
adjr2(b) # 0.7666107135723619

効率的なデータクリーニング

取り込んだデータをクリーニングする方法を紹介
　ということで，延々と説明が始まる
データを読みこんだだけで、すぐに分析ができる訳ではない
　ごもっとも。ちゃんとしたデータなら読み込んですぐ分析に入れます。
実際に使えるデータに整備するためには、前準備として「データクリーニング」が必要
read_csv を使う
エラーが出た
Warning: Missing column names filled in: 'X3' [3]
Warning: 265 parsing failures.
データが読み取れているかどうか最初の方を表示させて確かめる
head を使う
1 行目 (Last Updated Date 2019-03-21…) は不要
LibreOffice（もしくは Excel）を使ってファイルの様子を確かめる
　いやいや，最初から LibreOffice/Excel で読んでみればよいと思う
　それより最初は，head  や cat でファイルの内容を表示させてみればよい
　LibreOffice/Excel で読むのに適しているかわかる
　場合によっては単なるエディタで読んだほうがよいファイルもある
1行目から4行目（黄色い部分）をスキップ skip して読み込む
read_csv の引数を追加
Warning: Missing column names filled in: 'X64' [64]
　LibreOffice/Excel で読んでいれば，不要な部分は目で見て簡単に削除できる
str()関数を使ってデータが読み取れたかどうか確かめる
必要な変数だけを選ぶ
内容が不明な変数（X64）の内容を確認する
変数名にスペースが入っている場合、スペースを取り去った変数名に変更する
変数名にスペースが入っていると、RStudio が認識しないため
　これも，LibreOffice/Excel で読んだ段階で不要な部分は目で見て修正・削除できる
paged_table() 関数を使って、読み込んで修正を加えたデータの様子を表示させる
　Web ページ上でも綺麗に表示できるよということなんだろうけど，LibreOffice/Excel の手軽さには負ける
ワイド形式なので、実際に分析する際には ロング形式に変換する必要がある
ロング形式に変換されたデータを使うと様々な分析が可能になる
　今回のような場合には不要。鶏を割くに牛刀をもってす...
　ロング形式だとかえって面倒な場合もある
yearの class を文字列データ (chracter) から数値データ (numeric) に変換する
例えば、日本と中国の○○の推移を可視化してみる
日本と中国のデータだけを抜き出しとデータフレーム名をつける
　この後 ggplot を使うためということだが，単純な作業を複雑に行うことになる

ああだこうだいうよりも，代替法を実際にやってみよう

head でファイルの先頭を表示させてみる。

先頭に余計なものが付いていること，各行の最後が "," なのでちょっと気を付けよう，欠損値がかなりあるなあということがわかる。

まあ，CSV としては比較的まともなので，LibreOffice/Excel で読んでいいだろう。

bash-3.2$ head wb_gdp_pc.csv 
"Data Source","World Development Indicators",

"Last Updated Date","2019-03-21",

"Country Name","Country Code","Indicator Name","Indicator Code","1960","1961","1962","1963","1964","1965","1966","1967","1968","1969","1970","1971","1972","1973","1974","1975","1976","1977","1978","1979","1980","1981","1982","1983","1984","1985","1986","1987","1988","1989","1990","1991","1992","1993","1994","1995","1996","1997","1998","1999","2000","2001","2002","2003","2004","2005","2006","2007","2008","2009","2010","2011","2012","2013","2014","2015","2016","2017","2018",
"Aruba","ABW","GDP per capita (current US$)","NY.GDP.PCAP.CD","","","","","","","","","","","","","","","","","","","","","","","","","","","6472.50202920407","7885.79654466735","9764.78997879329","11392.4558105764","12307.3117378314","13496.0033851411","14046.5039965619","14936.827039276","16241.0463247117","16439.3563609282","16586.0684357542","17927.7496352086","

LibreOffice で読んでみたら，こうだった。

要らないところ先頭4行，B，C，D，BK 列を削除し，A1 の変数名から空白を除去すると以下のようになる。

こぎれいなデータファイルになったので，neat.csv と名前を付けて保存しよう。

df = read.csv("neat.csv") で，何の問題なく読める。

今回のような，中国と日本の経時変化をグラフに描くということだけならば，LibreOffice で読んだ段階で，中国と日本の 2 行を選択して新たな csv ファイルに保存するということをすれば，あるいは，単に図を描くのなら，LibreOffice でグラフを描いたっていいと思う。

Julia でデータの可視化

R での丁寧な説明があるページを見ながら，Julia でやってみた。

元のページでは，dplyr を使ったり，ロングフォーマットに変換したりで，かえってわかりにくいように思えた。

なお，いつも使っているのは using Plots だが，日本語が化けるので，PyPlot を使うことにした。

それでも最初は化けたのだが，Python で matplotlib を使って日本語は使えていたので，一度だけ以下の 2 行で使用するフォントファミリーを指定してやることにより Julia でも問題なく日本語が使えた。

using PyCall
PyCall.PyDict(matplotlib["rcParams"])["font.family"] = "IPAMincho"

さて，以下がプログラム。ファイルを読んで，必要行を抜き出して，プロットする。実に簡明。

この通りをやってみようと思う人はいないと思うが，必要ならば wb_gdp_pc.csv の在処はググってみてください。

using CSV, DataFrames, PyPlot
wb_gdp = CSV.read("wb_gdp_pc.csv", DataFrame, header=5);
japan = Matrix(filter(1 => x -> x == "Japan", wb_gdp))[1, 5:end-2];
china = Matrix(filter(1 => x -> x == "China", wb_gdp))[1, 5:end-2];
year = collect(1960:2017);
plot(year, china, "r.-", year, japan, "c^--", markersize=2, linewidth=0.5)
title("日本と中国の 1 人当たりGDP の推移：1960−2017")
xlabel("Year")
ylabel("1 人あたり GDP(US\$)")
legend(["China", "Japan"])

こんな図が描けました。

追記：

これをやった後で，using Plots でも日本語が使えるようになっていることを確認した。

ただし，明朝体を指定しているのにゴシックしか出ない。なぜだ。

using Plots
plot(collect(1:10),fontfamily="IPAMincho",
    label="日本語",title="日本語書けた")

Julia で使える R のデータセット

以下のリストは，R のパッケージ（赤字）に含まれているデータセット（黒字）の一覧である。

使用法

dataset("パッケージ名", "データセット名")

使用例

using RDatasets, DataFrames
df1 = dataset("COUNT", "affairs")
df2 = dataset("datasets", "iris")

COUNT
 affairs, azdrg112, azpro, badhealth, fasttrakg, lbw, lbwgrp, loomis, mdvis,
 medpar, rwm, rwm5yr, ships, titanic, titanicgrp

Ecdat
 Accident, Airline, Airq, Benefits, Bids, BudgetFood, BudgetItaly, BudgetUK,
 Bwages, CPSch3, Capm, Car, Caschool, Catsup, Cigar, Cigarette, Clothing,
 Computers, Cracker, Crime, DM,  Diamond, Doctor, DoctorAUS, DoctorContacts,
 Earnings, Electricity, Fair, Fatality, Fishing, Forward, FriendFoe, Garch,
 Gasoline, Griliches, Grunfeld, HC, HI, Hdma, Heating,  Hedonic, Housing,
 Icecream, Journals, Kakadu, Ketchup, Klein, LaborSupply, Labour, MCAS,
 Males, Mathlevel, MedExp, Metal, Mode, ModeChoice, Mofa, Mroz, MunExp,
 NaturalPark, Nerlove, OFP, Oil, PSID, Participation, PatentsHGH, PatentsRD,
 Pound, Produc, RetSchool, SP500, Schooling, Somerville, Star, Strike,
 StrikeDur, StrikeNb, SumHes, Tobacco, Train, TranspEq, Treatment, Tuna,
 UnempDur, Unemployment, University, VietNamI, Wages, Wages1, Workinghours,
 Yen, Yogurt, incomeInequality

HSAUR
 BCG, BtheB, CYGOB1, Forbes2000, GHQ, Lanza, agefat, aspirin, birthdeathrates,
 bladdercancer, clouds, epilepsy, foster, heptathlon, mastectomy, meteo,
 orallesions, phosphate, pistonrings, planets, plasma, polyps, polyps3,
 pottery, rearrests, respiratory, roomwidth, schizophrenia, schizophrenia2,
 schooldays, skulls, smoking, students, suicides, toothpaste, voting, water,
 watervoles, waves, weightgain, womensrole

HistData
 Arbuthnot, Bowley, Cavendish, ChestSizes, CushnyPeebles, CushnyPeeblesN,
 Dactyl, DrinksWages, Fingerprints, Galton, GaltonFamilies, Guerry, Jevons,
 Langren.all, Langren1644, Macdonell, MacdonellDF, Michelson, MichelsonSets,
 Minard.cities, Minard.temp, Minard.troops, Nightingale, OldMaps, PearsonLee,
 PolioTrials, Prostitutes, Pyx, Quarrels, Snow.deaths, Snow.polygons,
 Snow.pumps, Snow.streets, Wheat.monarchs, Yeast, YeastD.mat, ZeaMays

ISLR
 Auto, Caravan, Carseats, College, Default, Hitters, OJ, Portfolio, Smarket,
 Wage, Weekly

KMsurv
 aids, alloauto, allograft, azt, baboon, bcdeter, bfeed, bmt, bnct, btrial,
 burn, channing, drug6mp, drughiv, hodg, kidney, kidrecurr, kidtran, larynx,
 lung, pneumon, psych, rats, std, stddiag, tongue, twins

MASS
 Aids2, Animals, Boston, Cars93, Cushings, DDT, GAGurine, Insurance, Melanoma,
 OME, Pima.te, Pima.tr, Pima.tr2, Rabbit, Rubber, SP500, Sitka, Sitka89, Skye,
 Traffic, UScereal, UScrime, VA, abbey, anorexia, bacteria, beav1, beav2, biopsy,
 birthwt, cabbages, caith, cats, cement, chem, coop, cpus, crabs, eagles, epil,
 farms, fgl, forbes, galaxies, gehan, genotype, geyser, gilgais, hills, housing,
 immer, leuk, mammals, mcycle, menarche, michelson, minn38, motors, muscle,
 newcomb, nlschools, npk, npr1, oats, painters, petrol, quine, road, rotifer,
 ships, shoes, shrimp, shuttle, snails, steam, stormer, survey, synth.te,
 synth.tr, topo, waders, whiteside, wtloss

SASmixed
 AvgDailyGain, Multilocation, SIMS

Zelig
 PErisk, SupremeCourt, Weimar, approval, bivariate, coalition, coalition2,
 eidat, free1, free2, grunfeld, hoff, homerun, immi1, immi2, immi3, immi4,
 immi5, immigration, klein, kmenta, macro, mexico, mid, newpainters, sanction,
 swiss, tobin, turnout, voteincome

adehabitatLT
 albatross, bear, buffalo, capreochiz, capreotf, hseal, ibex, ibexraw, mouflon,
 porpoise, puechcirc, rupicabau, teal, whale

boot
 acme, aids, aircondit, aircondit7, amis, aml, bigcity, brambles, breslow,
 calcium, cane, capability, catsM, cav, cd4, channing, city, claridge,
 cloth, co.transfer, coal, darwin, dogs, downs.bc, ducks, fir, frets, grav,
 gravity, hirose, islay, melanoma, motor, neuro, nitrofen, nodal, nuclear,
 paulsen, poisons, polar, remission, salinity, survival, tau, tuna, urine

car
 AMSsurvey, Adler, Angell, Anscombe, Baumann, Bfox, Blackmore, Burt, CanPop,
 Chile, Chirot, Cowles, Davis, DavisThin, Depredations, Duncan, Ericksen,
 Florida, Freedman, Friendly, Ginzberg, Greene, Guyer, Hartnagel, Highway1,
 Leinhardt, Mandel, Migration, Moore, Mroz, OBrienKaiser, Ornstein,
 Pottery, Prestige, Quartet, Robey, SLID, Sahlins, Salaries, Soils, States,
 Transact, UN, USPop, Vocab, WeightLoss, Womenlf, Wool

cluster
 agriculture, animals, chorSub, flower, plantTraits, pluton, ruspini,
 votes.repub, xclara

datasets
 BOD, CO2, Formaldehyde, HairEyeColor, InsectSprays, LifeCycleSavings,
 Loblolly, OrchardSprays, PlantGrowth, Puromycin, Theoph, Titanic,
 ToothGrowth, UCBAdmissions, USArrests, USJudgeRatings, USPersonalExpenditure,
 VADeaths, WorldPhones, airquality, anscombe, attenu, attitude, cars,
 chickwts, crimtab, esoph, euro, faithful, freeny, infert, iris, islands,
 longley, morley, mtcars, occupationalStatus, precip, pressure, quakes,
 randu, rivers, rock, sleep, stackloss, swiss, trees, volcano, warpbreaks,
 women 

gamair
 aral.bnd, aral, bird, blowfly, bone, brain, cairo, chicago, chl, co2s,
 coast, engine, gas, harrier, hubble, ipo, mack, mackp, med, meh, mpg,
 prostate, sitka, sole, sperm.comp1, sperm.comp2, stomata, swer, wesdr,
 wine

gap
 PD, aldh2, apoeapoc, cf, crohn, fa, fsnps, hla, lukas, mao, mhtdata, nep499

ggplot2
 diamonds, economics, midwest, movies, mpg, msleep, presidential, seals

lattice
 barley, environmental, ethanol, melanoma, singer

lme4
 Dyestuff, Dyestuff2, InstEval, Pastes, Penicillin, VerbAgg, cake, cbpp,
 grouseticks, sleepstudy

mgcv
 columb

mlmRev
 Chem97, Contraception, Early, Exam, Gcsemv, Hsb82, Mmmec, Oxboys, ScotsSec,
 bdf, egsingle, guImmun, guPrenat, star

nlreg
 C1, C2, C3, C4, M2, M4, chlorsulfuron, daphnia, helicopter, metsulfuron, ria

plm
 Cigar, Crime, EmplUK, Gasoline, Grunfeld, Hedonic, LaborSupply, Males,
 Produc, Snmesp, SumHes, Wages

plyr
 baseball

pscl
 AustralianElectionPolling, AustralianElections, EfronMorris, RockTheVote,
 UKHouseOfCommons, absentee, admit, bioChemists, ca2006, iraqVote,
 politicalInformation, presidentialElections, prussian, unionDensity,
 vote92

psych
 Bechtoldt.1, Bechtoldt.2, Bechtoldt, Dwyer, Gorsuch, Harman.5, Harman.8,
 Harman.political, Holzinger.9, Holzinger, Reise, Schmid, Thurstone.33,
 Thurstone, Tucker, affect, bfi, blot, burt, cities, cubits, cushny,
 epi.bfi, galton, heights, income, iqitems, msq, neo, peas, sat.act,
 withinBetween

quantreg
 Bosco, CobarOre, Mammals, barro, engel, uis

reshape2
 french_fries, smiths, tips

robustbase
 Animals2, CrohnD, NOxEmissions, SiegelsEx, aircraft, airmay, alcohol,
 ambientNOxCH, bushfire, carrots, cloud, coleman, condroz, cushny, delivery,
 education, epilepsy, exAM, hbk, heart, kootenay, lactic, milk, pension,
 phosphor, pilot, possumDiv, pulpfiber, radarImage, salinity, starsCYG,
 telef, toxicity, vaso, wagnerGrowth, wood

rpart
 car.test.frame, car90, cu.summary, kyphosis, solder, stagec

sandwich
 PublicSchools

sem
 Bollen, CNES, Klein, Kmenta, Tests

survival
 bladder, cancer, cgd, colon, heart, kidney, leukemia, logan, lung, mgus,
 nwtco, ovarian, pbc, rats, stanford2, tobin, veteran

vcd
 Arthritis, Baseball, BrokenMarriage, Bundesliga, Bundestag2005,
 Butterfly, CoalMiners, DanishWelfare, Employment, Federalist, Hitters,
 HorseKicks, Hospital, JobSatisfaction, JointSports, Lifeboats, NonResponse,
 OvaryCancer, PreSex, Punishment, RepVict, Saxony, SexualFun, SpaceShuttle,
 Suicide, Trucks, UKSoccer, VisualAcuity, VonBort, WeldonDice, WomenQueue

Julia のベクトル（その 2）

並べ替え

d = [3, 2, 5, 1, 4];
sort(d) # 1, 2, 3, 4, 5; 昇順ソート
sort(d, rev=true) # 5, 4, 3, 2, 1; 降順ソート

値のある場所のインデックス

一番大きい要素のインデックス，二番目に大きい要素のインデックス，...

sortperm(d) # 4, 2, 5, 3, 1; 昇順ソートするときに，順に取り出す要素の位置
d[sortperm(d)] # sortperm(d) を使ってソートするのはこういうこと
sortperm(d, rev=true) # 3, 5, 1, 2, 4; 降順ソートするときに，順に取り出す要素の位置
d[sortperm(d, rev=true)] # 5, 4, 3, 2, 1; sortperm(d) を使ってソートするのはこういうこと
sort!(d) # 1, 2, 3, 4, 5; 破壊的ソート。d がインプレースでソートされるので，代入しなくてよい。
d # 1, 2, 3, 4, 5; 内容が変わっている

統計関数

using Statistics
using Random;
Random.seed!(12345); # 乱数の種
x = randn(100);

最大値

maximum(x) # 2.1307390315181536; 最大値
argmax(x)  # 55; 最大値のあるインデックス
findmax(x) # (2.1307390315181536, 55)); 最大値とそのインデックス

最小値

minimum(x) # -2.3280152457586003; 最小値
argmin(x)  # 73; 最小値のあるインデックス
findmin(x) # (-2.3280152457586003, 73); 最小値とそのインデックス

最大値と最小値

maximu(), minimu() 両方やるよりも遅い。

extrema(x) # (-2.3280152457586003, 2.1307390315181536);

和

sum(x) # 0.9821131440360702

平均値

mean(x) # 0.009821131440360701

不偏分散，分散

var(x) # 0.9954011110241598
var(x, corrected=false) # 0.9854470999139182

不偏標準編纂，標準偏差

std(x) # 0.9976979056929807
std(x, corrected=false) # 0.9926968821921011

中央値

median(x) # 0.13578110308214625
median!(x) # 0.13578110308214625; 破壊的（インプレース）関数

要約統計量

describe(x)
# Summary Stats:
# Length:         100
# Missing Count:  0
# Mean:           0.009821
# Minimum:        -2.328015
# 1st Quartile:   -0.684464
# Median:         0.135781
# 3rd Quartile:   0.690881
# Maximum:        2.130739
# Type:           Float64

y = [2, 4, 1, 3, 5];

差

diff(y) # 2, -3, 2, 2

累和

cumsum(y) # 2, 6, 7, 10, 15
[sum(y[1:i]) for i = 1:length(y)] # 2, 6, 7, 10, 15
z = similar(y);
cumsum!(z, y); # 2, 6, 7, 10, 15; 破壊的（インプレース）関数
y; # 2, 4, 1, 3, 5
z; 2, 6, 7, 10, 15

積

y = [2, 4, 1, 3, 5];
prod(y) # 120

cumprod, cummin, cummax

[prod(y[1:i]) for i = 1:length(y)] # 2, 8, 8, 24, 120; cumprod() に相当
[minimum(y[1:i]) for i = 1:length(y)] # 2, 2, 1, 1, 1; cummin() に相当
[maximum(y[1:i]) for i = 1:length(y)] # 2, 4, 4, 4, 5; cummax() に相当

不偏共分散，共分散

a = [1, 3, 4, 2, 6];
b = [3, 2, 5, 4, 3];
cov(a, b) # 0.15000000000000002; 不偏共分散
cov(a, b, corrected=false) # 0.12000000000000002; 共分散

ピアソンの積率相関係数

cor(a, b) # 0.06839411288813299; ピアソンの積率相関係数
A = randn(100, 3);
cor(A) # ピアソンの積率相関係数行列
# 3×3 Array{Float64,2}:
#  1.0         0.0227498   0.18433
#  0.0227498   1.0        -0.0643971
#  0.18433    -0.0643971   1.0

ノルム

using LinearAlgebra
norm(a) # 8.12403840463596; ノルム

内積

dot(a, b) # 55; 内積
sum(a .* b) # 55; 内積

外積

a = [1, 2, 3];
b = [10, 20, 30, 40];
a * b' # 外積
# 3×4 Array{Int64,2}:
#  10  20  30   40
#  20  40  60   80
#  30  60  90  120

using StatsBase  # geomean, harmmean, quantile, corspearman, corkendall

分位数（パーセンタイル値）

quantile(x)
  # -2.3280152457586003; 最小値
  #  -0.6844644232186423; 25%タイル値
  #   0.13578110308214625; 中央値（50%タイル値）
  #   0.690881410484019; 75%タイル値
  #   2.1307390315181536; 最大値
quantile(x, 0.00) # -2.3280152457586003; 最小値
quantile(x, 0.25) # -0.6844644232186423; 25%タイル値
quantile(x, 0.50) # 0.13578110308214625; 中央値（50%タイル値）
quantile(x, 0.75) # 0.690881410484019; 75%タイル値
quantile(x, 1.00) # 2.1307390315181536; 最大値
quantile!(x, 0.5) # 0.13578110308214625; 破壊的（インプレース）関数

幾何平均

geomean([2, 1, 3, 2, 1]) # 1.6437518295172258; 幾何平均

調和平均

harmmean([2, 1, 3, 2, 1]) # 1.5; 調和平均

スピアマンの順位相関係数

a = [1, 3, 4, 2, 6];
b = [3, 2, 5, 4, 3];
corspearman(a, b) # 0.10259783520851543; スピアマンの順位相関係数
A = randn(10, 3);
corspearman(A); # スピアマンの順位相関係数行列
# 3×3 Array{Float64,2}:
#   1.0         -0.430303   0.00606061
#  -0.430303     1.0       -0.139394
#   0.00606061  -0.139394   1.0

ケンドールの順位相関係数

corkendall(a, b) # 0.10540925533894598; ケンドールの順位相関係数
corkendall(A); # ケンドールの順位相関係数行列
# 3×3 Array{Float64,2}:
#   1.0        -0.288889    0.0222222
#  -0.288889    1.0        -0.0222222
#   0.0222222  -0.0222222   1.0

集合

v = [2, 2, 4, 4, 4, 4, 1, 3, 3, 3, 5, 5, 5, 5, 5];
A = Set(v); # 4, 2, 3, 5, 1; 集合（要素の順序はない。重複もない。）
B = Set([1, 3, 5, 7, 9]); 7, 9, 3, 5, 1
union(A, B) # 7, 4, 9, 2, 3, 5, 1; 和集合
intersect(A, B) # 3, 5, 1; 積集合
setdiff(A, B) # 4, 2; 差集合
symdiff(A, B) # 7, 4, 9, 2; 対称差集合
A == B # false; 等しい
1 in A # true; 集合に含まれる
0 in A # false
Set() # (); 空集合
empty(A) # (); 空集合にする

要素の置換

v = [2, 4, 1, 3, 5]; # 2, 4, 1, 3, 5
v[[1, 3, 5]] = [10, 30, 50];
v; # 10, 4, 30, 3, 50

ベクトルの結合

u = [1, 2, 3];
v = [4, 5];
w = [6, 7, 8, 9];
vcat(u, v, w); # 1, 2, 3, 4, 5, 6, 7, 8, 9
append!(u, v); # インプレースで結合
append!(u, w);
u; # 1, 2, 3, 4, 5, 6, 7, 8, 9

要素・ベクトルの挿入

a = [1, 2, 3];
pushfirst!(a, 0) # 0, 1, 2, 3; インプレースで先頭に挿入

a = [1, 2, 3];
push!(a, 99) # 1, 2, 3, 99; インプレースで末尾に挿入（追加）

a = [1, 2, 3];
insert!(a, 2, 99) # 1, 99, 2, 3; インプレースでインデックス位置に挿入

u = [1, 2, 3];
v = [4, 5];
[insert!(u, 3, v2) for v2 in reverse(v)]
u; # 1, 2, 4, 5, 3

ベクトルの要素の削除

v = [1, 2, 3, 4, 5];
popfirst!(v) # 1; インプレースで先頭の要素を削除して，返す
v; # 2, 3, 4, 5

v = [1, 2, 3, 4, 5];
pop!(v) # 5; インプレースで最後の要素を削除して，返す
v; # 1, 2, 3, 4

v = [1, 2, 3, 4, 5];
deleteat!(v, 3); # インプレースでインデックス位置の要素を削除(delete at)
v; # 1, 2, 4, 5

empty!(v); # 空ベクトルにする
v; # Int64[]

要素ごとの演算

v = [2, 4, 1, 3, 5];
w = [1, 3, 5, 7, 9];
v + w;  # 3, 7, 6, 10, 14
v .+ w; # 3, 7, 6, 10, 14
v - w;  # 1, 1, -4, -4, -4
v .- w; # 1, 1, -4, -4, -4
v .* w; # 2, 12, 5, 21, 45
v * w; # エラー　MethodError: no method matching *(::Array{Int64,1}, ::Array{Int64,1})
v ./ w; # 2.0, 1.3333333333333333, 0.2, 0.42857142857142855, 0.5555555555555556
v / w # 注！ 5 x 5 配列
# 5×5 Array{Float64,2}:
#  0.0121212   0.0363636  0.0606061  0.0848485  0.109091
#  0.0242424   0.0727273  0.121212   0.169697   0.218182
#  0.00606061  0.0181818  0.030303   0.0424242  0.0545455
#  0.0181818   0.0545455  0.0909091  0.127273   0.163636
#  0.030303    0.0909091  0.151515   0.212121   0.272727
v .+ 1; # 3, 5, 2, 4, 6
v + 1; # エラー MethodError: no method matching +(::Array{Int64,1}, ::Int64)
v * 2; # 4, 8, 2, 6, 10
v .* 2; # 4, 8, 2, 6, 10
v .^ 2; # 4, 16, 1, 9, 25
v ^ 2; # エラー MethodError: no method matching ^(::Array{Int64,1}, ::Int64)
1 ./ v # 0.5, 0.25, 1.0, 0.333, 0.2
1 / v # 0.0363636  0.0727273  0.0181818  0.0545455  0.0909091 ?? どういう計算？
v .< 3 # 1, 0, 1, 0, 0
mod.(v, 2) # 0, 0, 1, 1, 1
mod.(v, 2) .== 0; # 1, 1, 0, 0, 0
sqrt.(v); # 1.41421, 2.0, 1.0, 1.73205, 2.23607
exp.(v); # 7.38906, 54.5982, 2.71828, 20.0855, 148.413

定数

円周率

pi # π = 3.1415926535897... 円周率
π # π = 3.1415926535897...; # \pi<tab> で入力
typeof(π) # Irrational{:π}
1pi # 3.141592653589793
typeof(1pi) # Float64
BigFloat(pi) # 3.141592653589793238462643383279502884197169399375105820974944592307816406286198
BigFloat(1) + pi # 4.141592653589793238462643383279502884197169399375105820974944592307816406286198
BigFloat(1) + Float64(pi) # 4.141592653589793115997963468544185161590576171875
1 + pi # 4.141592653589793
Float64(pi) # 3.141592653589793

1024 ビットにする

setprecision(BigFloat, 1024) do
    BigFloat(pi)
end
# 3.14159265358979323846264338...24586997

ネイピア数

ℯ # ℯ = 2.7182818284590...; \euler<tab> で入力; ネイピア数
typeof(ℯ) # Irrational{:ℯ}
BigFloat(ℯ) # 2.718281828459045235360287471352662497757247093699959574966967627724076630353555
1ℯ # 2.718281828459045

 

Julia のベクトル（その 1）

ベクトル

Julia では，配列（特に1次元のときはベクトル，2次元のときは行列とも呼ぶ）は [ ] で囲んで使用する。

v = [1, 2, 3, 4, 5];
v # 1, 2, 3, 4, 5

対話モードで実行する場合（RERL; Run Evaluate Print Loop と呼ばれる），行の最後に ';' が付いているとその行の評価結果を表示しない。C 言語などのように，行末に必ず必要な ';' ではない。

typeof(v) # Array{Int64,1}

typeof 関数は，オブジェクトの型を表示する。上のように小数点なしの数を ',' で区切って定義すると整数ベクトルになる。Array は配列，Int64 は要素が整数，1 は次元数が 1 すなわちベクトルであることを意味する。

ちなみに，配列の次元サイズは size() で求めることができる。

size(v) # (5,)

次元が 1 のときは，上の例のように示される。
最後の ',' は (5) だと数値の 5 がカッコでくくられたもの（数式）と区別がつかないのでこのように示される。
つまり，この場合の '( )' は数式の中での '( )' ではなく，これ全体がタプル（複数の要素のまとまり）という型であることを意味している。

要素に一つでも小数点がついていると実数ベクトルになる。

v2 = [1.0, 2, 3, 4, 5]; # 1.0, 2.0, 3.0, 4.0, 5.0

整数ベクトルを実数ベクトルにするときは，以下のように型を指定して変換することができる。

v3 = Array{Float64,1}(v); # 1.0, 2.0, 3.0, 4.0, 5.0

Float64 で，実数型であることを示している。

typeof(v3) # Array{Float64,1}

または，

v4 = convert.(Float64, v); # 1.0, 2.0, 3.0, 4.0, 5.0

convert は最初の引数で示した型に変換する関数。'(' の前の '.' は，第 2 引数が複数の要素を持つ場合，それぞれの要素に作用することを指示する。もし，この '.' がないとエラーになる。

空白で区切ると横ベクトルになる。

v5 = [10 20 30 40 50];
# 1×5 Array{Int64,2}:
#  10  20  30  40  50

しかしこれは 1 行 n 列の配列（行列）である。typeof() で型を見ると次元数が 2 になっていることで確認できる。

typeof(v5) # Array{Int64,2}
size(v5)   # (1, 5)

v の定義のときのように ',' で区切って定義した縦ベクトルを n 行 1 列の配列（行列）にするには，配列の形を変える関数 reshape を使って次のようにする。

reshape(v, 5, 1) # 1, 2, 3, 4, 5

第 2，第 3 引数の 5, 1 は 結果として返す配列の第 1，第 2 次元の大きさを示すので，5 行 1 列の配列になる。
一般的には以下のようにするとよい。length(v) は v の長さ，':' は 2 次元配列の場合は，1 次元目の大きさが決まれば自動的に決まるので，特に示す必要がないときに指定するとよい。

reshape(v, length(v), :)

逆に ' ' で区切った 1 行 n 列の配列をベクトルにしたいときには vec 関数を使う。

vec([10 20 30 40 50]) # 10, 20, 30, 40, 50

要素が文字列である配列（ベクトル，行列）も同じように定義できる。typeof 関数で示される String は要素が文字型であることを示している。

v8 = ["A", "B", "C", "D", "E"]; # "A", "B", "C", "D", "E"
typeof(v8) # Array{String,1}

Julia では，文字列は長さが 1 であっても，二重引用符で囲まなければならない。
以下のように，1 文字を一重引用符で囲むと文字列型ではなく文字型 Char になる。逆に言えば，長さが2以上の文字列を一重引用符で囲むとエラーになる。

v7 = ['1', 'A', 'b'];
# 3-element Array{Char,1}:
#  '1': ASCII/Unicode U+0031 (category Nd: Number, decimal digit)
#  'A': ASCII/Unicode U+0041 (category Lu: Letter, uppercase)
#  'b': ASCII/Unicode U+0062 (category Ll: Letter, lowercase)
typeof(v7) # Array{Char,1}

規則性のあるベクトル

等差数列

collect(1:1:10) # 1 〜 10，初項 1，公差 1; 第2引数が公差
collect(1:10) # 1 〜 10，初項 1，公差 1; 第2引数が1の場合は省略可
collect(10:-1:1) # 10 〜 1，初項 10，公差 -1; 降順の場合も公差は明示する
collect(1:5:20) # 1, 6, 11, 16，初項 1，公差 5; 第2引数が公差
collect(0:0.2:3) # 初項 0, 公差 0.2の数列; 最終項は 3

collect は要素を実際にメモリー上に確保する。
for i = 1:10 のような使い方（イテレータ）のときには，メモリー上に確保しない。

range(1, 10, step=2) # 1:2:9; イテレータを作成する
range(1, 10, length=20) # 1.0:0.47368421052631576:10.0
collect(range(1, 10, length=3)); # 1.0, 5.5, 10.0

繰り返し

repeat([0], 5) # 0 が 5個; repeat(0, 5) はエラーになる
repeat([0, 1, 2], 5) # 0, 1, 2 が 5 回繰り返されるベクトル
repeat(["foo", "bar", "baz"], 5) # "foo", "bar", "baz" が 5 回繰り返されるベクトル

A = repeat([10 20 30], inner=(1,2,1)) # 10  10  20  20  30  30
B = repeat([10 20 30], inner=(1,4,1)) # 10  10  10  10  20  20  20  20  30  30  30  30

inner を指定する repeat 関数，repeat(v, inner=(1, m, 1)) は 1 × (length(v)×m) × 1 の3次元配列である。

size(A) # (1, 6, 1)
size(B) # (1, 12, 1)

vec 関数で，ベクトルに変換できる。

vec(A) # 10  10  20  20  30  30
vec(B) # 10  10  10  10  20  20  20  20  30  30  30  30

ゼロベクトル

zeros(5) # 実数 0.0 が 5 個のベクトル
zeros(Int64, 5) # 整数 0 が 5 個; 64 ビット環境では zeros(Int, 5) と同じ

1 ベクトル

ones(5) # 実数 1.0 が 5 個のベクトル; fill(1.0, 5) と同じ
ones(Int64, 5) # 整数 1 が 5 個; 64 ビット環境では zeros(Int, 5) と同じ

任意の要素の繰り返しを持つベクトル

fill(0, 5) # zeros(Int64, 5) と同じ
fill(1, 5) # ones(Int64, 5) と同じ
fill(0.0, 5) # zeros(5) と同じ
fill(1.0, 5) # ones(1) と同じ
fill("a", 5) # "a" が 5 個のベクトル

乱数ベクトル

using Random;
Random.seed!(12345); # 乱数の種の設定。本当に乱数が必要なときは使わない。

以下では，途中で乱数を使わないときには，示されたのと同じ結果になる。

無作為抽出

ベクトルから要素をランダムに取り出す。

rand(1:5, 10) # 3, 4, 1, 2, 2, 5, 5, 4, 5, 2; 要素 1〜5 から 10 個を無作為抽出したベクトル
rand(["foo", "bar", "baz"], 5) # "baz", "bar", "bar", "foo", "bar"

乱数

一様乱数を要素とする，長さ 3 のベクトル

rand(3) # 0.36007770681119133, 0.41511087872028374, 0.8325976767993883

正規乱数を要素とする，長さ  3 のベクトル

randn(3) # -0.55971756622102, -1.4990404969589635, 1.6739451586233716

要素の取得

Julia では，R と同様先頭要素のインデックスは 1 である。

v11 = [10, 20, 30, 40, 100]; # 10, 20, 30, 40, 100
v11[1] # 10
v11[begin] # 10; 'begin' はインデックス 1 を表す
v11[end] # 100; 'end' は最後のインデックスを表す
v11[end-2] # 30; 'begin', 'end' を算術演算してインデックスにできる
v11[[2,3,end]] # 20, 30, 100 # ベクトルで指定して要素を取り出すことができる
v11[[1,2,1,3,5,5]] # 10, 20, 10, 30, 100, 100; 何回でもどこからでも取り出せる
v11[1:2:end] # 10, 30, 100; イテレータで要素を取り出すことができる
v11[v11 .> 15] # 20, 30, 40, 100; 条件式（比較演算）を満たす要素を取り出せる

'.>' の '.' は複数の要素を持つオブジェクトの各要素に適用するために必要である。
なければエラーになる。

BitArray 型のベクトル

v11 .> 15 # BitArray 型のベクトル [0, 1, 1, 1, 1]

配列のインデックスに指定すると，1 に対応するオブジェクトの要素が取り出される。

v11[15 .< v11 .< 35] # 20, 30; 条件式（論理演算））を満たす要素を取り出せる

'15 .< v11 .< 35' は '15 .< v11 .& v11 .< 35' と同じではない。

以下の 2 つは 同じことを表していそうに思うが，結果が違う。

v11[15 .< v11 .& v11 .< 35] # 20, 30; 両方の条件を満たす要素を取り出す。 & は and
v11[v11 .> 15 .& v11 .< 35] # 20, 30, 40, 100; 両方の条件を満たす要素を取り出す。 & は and

'15 .< v11 .< 35' は '15 .< v11 .& v11 .< 35' と同じになる。'.&' とすることに注意。
'15 .< v11' と 'v11 .> 15' は同じ意味なので，
'v11[15 .< v11 .& v11 .< 35]' と 'v11[v11 .> 15 .& v11 .< 35]' の結果が違うのは一見納得ができない。

実は，ここでの '&' はビット演算子であり，比較演算子 '.<'，'.>' より優先順位が高い。
'v11 .& v11' は '[10, 20 ,30, 40, 100]'，で
'v11[15 .< [10, 20 ,30, 40, 100] .< 35]' が '[20, 30]' になる。

一方，'15 .& v11' は '[10, 4, 14, 8, 4]' で，
'v11[v11 .> [10, 4, 14, 8, 4] .< 35]' が '20, 30, 40, 100' になる。当然，予期せぬ間違った答えである。

以下のように，比較演算を先に行うようにカッコを補えば，正しい答えになる。

v11[(v11 .> 15) .& (v11 .< 35)] # 20, 30

また，'v11[15 .< v11 .& v11 .< 35]' もプログラムを書いた意図からは 'v11[(15 .< v11) .& (v11 .< 35)]' とすべきである。

以下のような単純な論理式の例だとわかりやすいかもしれない。

1 < 3 & 4 < 5 # false
(1 < 3) & (4 < 5) # true

1 > 3 | 5 < 6 # false
(1 > 3) | (5 < 6) # true

Python の場合も R の場合も，論理演算子のが出現する可能性のあるところにビット演算子も出現する可能性はないので，Julia のときのような誤解は生じないのである。

R にはビット演算はないので，'&' は論理演算子として解釈されるので，比較演算をカッコでくくらなくても正しい答えになる。

using RCall
R"""
v11 = c(10, 20, 30, 40, 100)
v11[v11 >= 15 & v11 <= 35] # 20, 30
"""

要素への代入

インデックスを指定してスカラー，ベクトルを代入することができる。

v = collect(1.:10);
v[1] = 10;
v[2:5] = [20, 30, 40 ,50];
v[9] = NaN;
v # 10.0  20.0  30.0  40.0  50.0  6.0  7.0  8.0  NaN  10.0

脇道

欠損値は 'missing' で表すが，上のようなベクトル v の要素として代入はできない。

'v[8] = missing' はエラーになる。

'missing' も許す型（たとえば 'Union{Missing, Int64}'）のオブジェクトならば，可能である。

x = Array{Union{Missing, Int64},1}([3,2,1,4,3,5,missing,6]);
x[4] = missing;
x #  3, 2, 1, missing, 3, 5, missing, 6

たとえば，以下のようにすると 'missing' を初期値として持つ長さ 10 のベクトルを定義できる。

y = Array{Union{Missing, Float64},1}(undef, 10)

1, 4, 7, 10 番目の要素を二乗する。

for i =1:3:length(y)
  y[i] = i^2
end
y; # 1, 2, 1, 16, 5

'missing' を取り除いたベクトルを求める。

z = collect(skipmissing(y)); # 1, 4, 1, 16, 5
using Statistics # 単に mean を使うだけなのにこれが必要
mean(z) # 41.5; missing のないベクトルの平均値
mean(y) # missing
mean(skipmissing(y)) # 41.5; missing を除いたベクトルの平均値

Julia に NA はない。

NA; # UndefVarError: NA not defined

Julia の NaN は？

NaN # NaN
NaN == missing # missing
NaN * 10.3 # NaN
NaN * Inf # NaN
NaN & missing # MethodError: no method matching &(::Float64, ::Missing)

ベクトルの属性

a = [1, 2, 3, 4, 5];

長さ

length(a) # 5

検索

要素はあるか？

2 in a # true
9 in a # false

何個あるか？

b = [0, 1, 2, 1, 2, 3, 4, 3, 5, 2, 2, 1, 1, 5, 0];

count の第一引数は関数である。
比較演算式，論理演算式も関数なので，count の第 1 引数になれる。

5 < 10 # true
<(5, 10) # true

count(isequal(1), b) # 4; 1 と等しいか？
count(b .== 1) # 4; b は 1 と等しいか？
count(isequal.(b, 1)) # 4; b と 1 が等しいか？ 等しければ 1 等しくなければ 0
count(iseven, b) # 7; 偶数か？
count(mod.(b, 2) .== 0) # 7; 偶数か？
count(isless.(b, 2)) # 6; b は 2 より小さいか？　isgreater() はない
count(b .< 2) # 6; b は 2 より小さいか？

最大値，最小値は何か？どこにあるか？

argmax(b) # 9; 最大値が最初に現れるインデックス
b[argmax(b)] # 5; 最大値
maximum(b) # 5; 最大値
argmin(b) # 1; 最小値が最初に現れるインデックス
b[argmin(b)] # 0; 最小値
minimum(b) # 0; 最小値

ある値がどこにあるか？

findall 関数の第 1 引数は，関数である。論理演算子も関数である。λ関数（無名関数）でもよい。

findall(==(1), b); # 2, 4, 12, 13; 1 と等しい要素のインデックス
findall(==(b[argmax(b)]), b) # 9, 14; 最大値が現れる全てのインデックス
findall(==(minimum(b)), b) # 9, 14; 最大値が現れる全てのインデックス
findall(==(b[argmin(b)]), b) # 1, 15; 最小値が現れる全てのインデックス
findall(==(minimum(b)), b) # 1, 15; 最小値が現れる全てのインデックス

findall(isodd, b) # 2, 4, 6, 8, 9, 12, 13, 14; 要素が奇数のインデックス
b[findall(isodd, b)] # 1, 1, 3, 3, 5, 1, 1, 5; その要素の列挙
findall(x -> 2 < x < 4, b) # 6, 8; 要素を x として 2 < x < 4 を満たす要素のインデックス（λ関数）
b[findall(x -> 2 < x < 4, b)] # 3, 3
c = findall(x -> sqrt(x) == floor(Int, sqrt(x)), b) # 1, 2, 4, 7, 12, 13, 15; 要素が平方数であるインデックス
b[c] # 0, 1, 1, 4, 1, 1, 0; その要素

Julia で注意しておかないといけないこと（その 1）

以下の 2 つの論理式は true にはならない

1 < 3 & 4 < 5 # false
1 > 3 | 5 < 6 # false

理由を知りたい人は，ずっと下にスクロールする

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

以下のように書かなければならない

(1 < 3) & (4 < 5) # true
(1 > 3) | (5 < 6) # true