R言語

ウロウロしてたらpythonの歴史に触れている方がいた。 なぜ機械学習にPythonが選ばれるのか

Rに追いつけ、追い越せが年表に良く表れているな。思う所pythonは力任せにぶん回す機械学習に進出したんだな。対してRは学術的な匂いがする。まあ、棲み分け、使い分けをすれば良い事なんで、これ以上は黙っていよう。

前々回にやったpythonで回帰をウブでも出来るようにする。ウブはほとんでパッケージが空なので、まずはpip3から入れる

curl https://bootstrap.pypa.io/get-pip.py -o get-pip.py
sudo python3 get-pip.py
pip3 install -U pip          ## for update

後はpip3を使って適当にインストールしてく。

sakae@ub:/tmp$ sudo pip3 install scipy
Collecting scipy
  Downloading scipy-1.5.2-cp38-cp38-manylinux1_x86_64.whl (25.7 MB)
     |████████████████████████████████| 25.7 MB 302 kB/s
Collecting numpy>=1.14.5
  Downloading numpy-1.19.1-cp38-cp38-manylinux2010_x86_64.whl (14.5 MB)
     |████████████████████████████████| 14.5 MB 12.0 MB/s
Installing collected packages: numpy, scipy
Successfully installed numpy-1.19.1 scipy-1.5.2

最後にipython3を趣味でいれた。いざ実行

TypeError: Couldn't find foreign struct converter for 'cairo.Context'

描画用の窓は出て来るものの、こんなエラーも出て来た。

sakae@ub:/tmp$ sudo pip3 install pyqt5

ぐぐったら、これも入れておけとの事。これ、真面目にDesktop環境を作らなかった弊害で、GKT3用の環境が出来ていなかったものと思わる。

/usr/local/lib/python3.8/dist-packages/matplotlib/mpl-data/matplotlibrc

## ***************************************************************************
## * BACKENDS                                                                *
## ***************************************************************************
## The default backend.  If you omit this parameter, the first working
## backend from the following list is used:
##     MacOSX Qt5Agg Gtk3Agg TkAgg WxAgg Agg

ここまではいいんだけど、ipython3を起動して、ひょいと図を書こうとすると、いちいち　importの呪文が必要。面倒嫌い。起動時に必要な物は自動読み込みして欲しい。pytthonで、初期化ファイルってどうやるの？

export PYTHONSTARTUP="/home/sakae/.pythonrc"

環境変数を経由して初期化ファイルを指定するようだ。自由度が高いな、pythonらしくないぞ、と、嫌味を言ってみる。

sakae@ub:~$ cat $PYTHONSTARTUP
import matplotlib.pyplot as plt
import numpy as np

後は、こんな具合に、よく利用されてるのを登録しておけば良い。 注意として、余り登録しすぎると、起動が遅くなるってのがあるな。そんな場合は、起動時に、-E オプションを付加し、全てのPYTHON* 環境変数を無視すれば良い。 See Also コマンドラインと環境

(setq python-shell-interpreter "ipython"
       python-shell-interpreter-args "--simple-prompt -i")

• 折角ipythonも入れたし、統計実習環境も完成してるんで、emacsから操れるように設定。コードを書いた後、C-c C-p すればipythonが起動してくる。後は、C-c C-c でbufferを評価するも佳し、C-c C-l しても良い。シンボルにカーソルを置きC-c C-fで調べるなんてのも有り。

ipython側に居るならば、 run hoge.py とかして、実行も出来る。

前回図書館から借りてきたEXCELで統計ってやつは、土台のEXCELがオイラーの感性にマッチせず、ろくに読まずに返却しちゃった。代わりに、日本統計学会公式認定　統計検定2級対応　統計学基礎 なる本が有ったので借りた。年寄りでも地方税をしっかりと徴収されるんで、元を取らねばと、さもしい根性です。おかげで、現役時代に多かった書籍代がめっきり少なくなりましたよ。

でも、あの手この手で、金を出させようと言う試みが多いですなあ。資格試験なんてのもその一つ。リナの経験十分有りますって、面接時に力説しても、試験官は困ってしまう。そんな時、検定１級ですとか言えば、まあ信用されるかな。

ボラクルの検定とかM$の検定とかその最たるもの。時代を反映してパイソン検定とかもあるな。 そして、統計検定 なるもののある。

受験料が微妙な値段。それより資料に使う書籍代とかで儲ける積りなのかしらん。 申し訳程度に、過去問が掲載されてた。コンピュータを使って、技能検定ってしないの？

まあ、いいや。

上記の本の付録が公開されてた。本文中に乗っている図とかはRを使っているので、それのコードやらが載ってるんだ。

その前に、Rの使い方が軽く説明されてた。

> bld <- read.csv(file.choose(), header=F)
Enter file name: hoge.csv
Error in file(file, "rt") : cannot open the connection
In addition: Warning message:
In file(file, "rt") :
  cannot open file 'hoge.csv': No such file or directory
> bld <- read.csv(file.choose(), header=F)
Enter file name: current.csv
> head(bld)
        V1  V2 V3 V4
1 12010104 118 73 57
2 12010121 109 70 71
:

read.csvで読み込む時、ファイル名を直接指定せず、ユーザーにその場で問う事が出来るぞう。 知らなかったぞ。

例の血圧データ。測定年月日時は、面倒なので、ひとまとめにしてる。けど、色々分析する場合、年・月・日・時に分解しておいた方が便利だろう。Rの中でやる？　そんな面倒をかけるのは、Windowsを使っている人だけです。大体、Rにしろpythonにしろ主流ユーザーはWindowsな人ですからね。ツールが充実してるunixならsedで一発変換です。

#! /bin/sh
cat current.csv | \
    sed 's/\(..\)\(..\)\(..\)\(..\)\(.*\)/\1,\2,\3,\4\5/' | \
    head -4

今回みたいに、パターンマッチが複数ある場合、一応emacsでスクリプトを書きます。 touch con.sh; chmod 755 con.sh; emacs con.sh します。コードを書いたら、C-c C-x して、出来栄え確認。よかったら、headを止めて、ファイルに書き出すようにします。

sedで面倒なのは、パターンを指定する括弧にエスケープが必要な事。エスケープしながら書いていると集中力が低下する。よって、エスケープ無しで最初に書いてします。後で、( -> \( のように、一気に変換しちゃう。この方が楽です。

さてDLしたzipを覗いてみるか。zipファイルはWindows御用達なんで、用心の為にdirを作って、その中で展開します。そうしておけば、フィアルをぶちまけられても掃除が簡単ですからね。 前回見繕っておいたunzip.pyを使って展開しました。これ感じ悪い漢字ファイル名にも対応してますからね。

.
└── 図(公開ソース)
    ├── 1章
    │   ├── 本文
    │   │   ├── death.csv
    │   │   ├── population_h26.csv
    │   │   ├── room.csv
    │   │   ├── Rスクリプト（本文）.txt
    │   │   └── Salary.csv
    │   └── 練習問題
    :

思った通り漢字が使われていました。GUIな人ならダブルクリックだけで何とかなるんだけど、CUIな人には面倒な奴です。dirも含めて何とかしよう。こういう場合に威力を発揮するのがemacsに備え付けのdired-modeと言うファイラーです。 Emacs dired wdired インストール 設定 使い方

取り合えず実験と言う事で、１章本文用にsec1なんてdirを用意。そこへ *.csvをコピー（mでマークしておいてからC）。Rスクリプト（本文）.txt と言う、どうしようもないやつは、変換したいファイルを選んでおいてからRを実施。renameしたよ。(s1.R)

で、喜びいさんで実行。roomの所で C-c C-c。emacsはsjisファイルでも、ちゃんとデコードして表示してくれる。

#データの入力
room <- read.csv("room.csv")
population_h26 <- read.csv("population_h26.csv")
death <- read.csv("death.csv")
Salary <- read.csv("Salary.csv")

下記のようなエラーが発生したぞ。

> Error in make.names(col.names, unique = TRUE) :
  invalid multibyte string at '<8b>߂<b3>'
> Error in make.names(col.names, unique = TRUE) :
  invalid multibyte string at '<93>s<93><b9><95>{<8c><a7><96><bc>'

csvファイル内に日本語データが混じっているんだな。はて、どうしたものか？

red.csvの所にカーソルを合わせて　C-c C-v でマニュアルを引きます。

Usage:

     read.table(file, header = FALSE, sep = "", quote = "\"'",
                  :
                fileEncoding = "", encoding = "unknown", text, skipNul = FALSE)

このように、山のような省略時引数が表示された。いかにもって言うfileEncodingとただのencodingが有る。どっちだろう。ずっと見てくとfileEncodingが効きそう。

room <- read.csv("room.csv",fileEncoding="cp932")
population_h26 <- read.csv("population_h26.csv",fileEncoding="cp932")
death <- read.csv("death.csv",fileEncoding="cp932")

んな訳で、上記のように修正。これでちゃんとcsvを読み込んでくれた。後はC-c C-cで、次々と評価していけばよい。らくちんである。

#図1.26のグラフ
SalaryofTokyo.ts<-ts(Salary[,3]/10000,start=c(2000,1),frequency=12)
plot(decompose(SalaryofTokyo.ts), xlab="年月")

#図1.27のコレログラム
acf(Salary$Tokyo,main="コレログラム：東京",xlab="月差")

面白いと思ったのは上のやつ。東京のサラリーマンの月収データを、時系列分析してる。

統計と言うと、確率なんたらとか検定なんたらがやたら多いけど、世の中って、身長とかテストの点をこねくり回すのはまれだろう。それより、株価の動きはどう？　とか、台風来てるけと河川が過去に溢れた事ってないのかしらとか、の方が切実。これらは時系列分析にうってつけ。でも、余り取り上げてられていないんだよな。

ウロウロしてたら、 Tokyo.R とか 勉強会発表内容一覧 なんてのに出会った。軽く読めて、暇つぶしにはうってつけ。

大体、R言語の言語としての特徴は何だろう。扱う分野が統計って事で、俯瞰する事を忘れていたな。そんな場合は、これを見ておけ。 R言語(by Wikipedia)

ベクトル処理言語ですって。一瞬スパコンの事かと思ったぞ。numpyより楽。octaveやjuliaと同類だな。

sum(runif(s <- 100000)^2 + runif(s)^2 <= 1) * 4 / s

piを求めるシュミレータをがこんな風に書けるとな。

1:100 -> n -> Ans
Ans[n %% 3 == 0 -> FizzSet] <- "Fizz"
Ans[n %% 5 == 0 -> BuzzSet] <- "Buzz"
Ans[FizzSet & BuzzSet] <- "FizzBuzz"
cat(Ans)

forとかは嫌いです。。。 そして、Rの巨大な掲示板のリンクが紹介されてた。RjpWiki 困った時は、ここを見ろってのは、昔の2CHに通じるものが有るな。

Rには膨大な実習用サンプルデータが付属している。例えばアマチュア無線家なら気になるデータが用意されていたりする。

> sunspot.month
       Jan   Feb   Mar   Apr   May   Jun   Jul   Aug   Sep   Oct   Nov   Dec
1749  58.0  62.6  70.0  55.7  85.0  83.5  94.8  66.3  75.9  75.5 158.6  85.2
1750  73.3  75.9  89.2  88.3  90.0 100.0  85.4 103.0  91.2  65.7  63.3  75.4
 :

太陽黒点の観測データだ。これを使って、今後どうなるか占うなんて、楽しい試行だ。残念ながら最新のデータは欠落してるけど。データの始点が1749年になってるけど、そんな昔から黒点に注目してたんか？　ちょいと信じ難い。こういう時はman、もとえ help(sunspot.month)

Source:

     WDC-SILSO, Solar Influences Data Analysis Center (SIDC), Royal
     Observatory of Belgium, Av. Circulaire, 3, B-1180 BRUSSELS
     Currently at <URL: http://www.sidc.be/silso/datafiles>

他には、このデータを使っての遊び方まで指南されてたぞ。

> bld <- read.csv("bld.csv")
> tail(bld)
     yy mm dd hh  hi lo pl
6289 20  8 29  5 149 72 49
6290 20  8 29 21 133 67 59
6291 20  8 30  5 138 76 47
6292 20  8 30 21 133 72 58
6293 20  8 31  5 142 76 48
6294 20  8 31 21 138 71 59

> table(bld$hh)

   2    3    4    5    6   20   21   22   23
   2  356 1821  923   47    5 3137    2    1

> table(bld$yy, bld$hh)

       2   3   4   5   6  20  21  22  23
  12   0  90 225  40   5   2 358   0   0
  13   1  72 214  67   4   2 357   1   0
  14   0  59 248  55   0   1 359   1   0
  15   0  82 237  45   1   0 364   0   1
  16   1  30 247  88   0   0 363   0   0
  17   0  15 176 168   6   0 364   0   0
  18   0   3 149 182  31   0 365   0   0
  19   0   4 185 175   0   0 364   0   0
  20   0   1 140 103   0   0 243   0   0

> subset(bld, hh >= 22)
     yy mm dd hh  hi lo pl
1344 13 11 14 22 119 65 69
2085 14 11 23 22 119 66 66
2469 15  6  3 23 120 75 64

> am <- subset(bld, hh < 12)
> dim(am)
[1] 3149    7

> tail(am, n= 4)
     yy mm dd hh  hi lo pl
6287 20  8 28  5 135 72 48
6289 20  8 29  5 149 72 49
6291 20  8 30  5 138 76 47
6293 20  8 31  5 142 76 48

+ aggregate(x=am[c("hi","lo","pl")], by=list(yy=am$yy,mm=am$mm), mean)
+
>     yy mm       hi       lo       pl
1   12  1 125.0645 79.32258 59.29032
2   13  1 124.0323 77.70968 58.32258
:
103 18 12 131.2258 72.54839 55.90323
104 19 12 138.8710 74.93548 53.70968

+ aggregate(hi ~ yy + mm, data=am, FUN=mean)
+
>     yy mm       hi
1   12  1 125.0645
2   13  1 124.0323

 + qqq <- aggregate(hi ~ yy + mm, data=am, function(ps) signif(mean(ps), 4))
!> head(qqq, n=2)
   yy mm    hi
 1 12  1 125.1
 2 13  1 124.0

+ qqq[order(qqq$yy,qqq$mm),]
+
>     yy mm    hi
1   12  1 125.1
10  12  2 123.3
19  12  3 121.0
:

 > qqq$yy <- qqq$yy + 2000
!> head(qqq,n=3)
     yy mm    hi
 1 2012  1 125.1
 2 2013  1 124.0
 3 2014  1 122.0