プログラミング書法 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 でグラフを描いたっていいと思う。