2019年7月25日木曜日

[R] 正則化パラメータの選択

Lasso推定値は正則化パラメータλの値に依存するため,λの値が異なると当然Lasso推定値も異なってきます.
特に,Lasso推定値の違いは変数選択の結果の違いに直接影響してくるため,正則化パラメータの値の選択は重要です.

前回の投稿(アメリカの犯罪例にLassoを適用)ではアメリカの犯罪データにLassoを適用しました.
その続きとして,用意した正則化パラメータの値ごとに10分割交差検証法を用いて得られた結果を図示すると下図のようになります.
上図において,y軸は誤差2乗和(交差検証法の値)で,x軸は正則化パラメータの対数値になります.赤色の点から伸びている帯は,10分割交差検証法で計算された10個の検証誤差の標準誤差を表します.また,縦の破線は,選択された正則化パラメータの位置を表します.

CVの値が最小になる正則化パラメータの値(λ = 3)に対応するモデルを最適なモデルとして選択しています.

上図を表示するには,前回の投稿(こちら)に引き続いて以下のように入力します.

まずはCVを計算します.
> res.cv <- cv.glmnet(x=X, y=y)

続いて,CV値の推移をプロットします.
> plot(res.cv, xlab="Logarithmic value of regularization parameter", ylab="Regression coefficient")

CV値が最小となる正則化パラメータの値を出力します.
> res.cv$lambda.min
[1] 15.15854

1標準誤差ルールにより選択された正則化パラメータの値を出力します.

> res.cv$lambda.1se

[1] 155.1523

2019年7月24日水曜日

[R] アメリカの犯罪例に Lasso を適用

スパース推定法による統計モデリング (統計学One Point)の 1.3 にある例を実行してみた結果のメモです.

アメリカ合衆国の都市における犯罪データに関して Lasso を適用する.

まずは,Lassoの代表的なパッケージである glmnet を読み込みます.glmnet パッケージのインストール方法はこちら
> ### Lasso
> library(glmnet)
Loading required package: Matrix
Loading required package: foreach
Loaded glmnet 2.0-18

基となるアメリカ合衆国の都市における犯罪データ(crime.txt *)には,以下の6つの項目を50の都市から取得した結果が記載されています.なお,crime.txtはこちらからダウンロードできるデータです.

X1 = total overall reported crime rate per 1 million residents
         人口100万人あたりの犯罪率(目的変数 Y とする)
X2 = reported violent crime rate per 100,000 residents
         人口10万人当たりの暴力的な犯罪率 ← 今回の例では使用しない
X3 = annual police funding in $/resident
         警察の年間資金($ /住民)(説明変数 X1 とする)
X4 = % of people 25 years+ with 4 yrs. of high school
         25歳以上で,高校を卒業した人の割合(説明変数 X2 とする)
X5 = % of 16 to 19 year-olds not in highschool and not highschool graduates.
         16〜19歳のうち,高校に通っていない,もしくは卒業していない人の割合(説明変数 X3 とする)
X6 = % of 18 to 24 year-olds in college
         18〜24歳のうち,大学生の割合(説明変数 X4 とする)

X7 = % of people 25 years+ with at least 4 years of college
         25歳以上で、4年制大学を卒業した人の割合(説明変数 X5 とする)

このデータに対して,X1のデータを目的変数(Y)とし,X3 - X7の5つを説明変数(X1 - X5)として回帰分析を行います.

まずは,データ(crime.txt)を読み込みます.
> crime <- read.table("crime.txt")  # Load crime data

"as.matrix" を用いて行列に変換します.
> crime <- as.matrix(crime)

どのデータを使用するか,また,目的変数,説明変数の設定を行います.
> X <- crime[, 3:7]  # Explanatory variable
> y <- crime[, 1]    # Objective variable

説明変数の標準化,目的変数の中心化を行います.
> X <- scale(X)      # Standardize explanatory variables
> y <- y - mean(y)   # Centralize objective variables
ただし,glmnet パッケージでは,説明変数を標準化しなくても,初期設定で標準化され,係数の推定値は元の尺度に戻して出力されるとのことです.

Lasso推定を行います.
> # Lasso estimation
> res <- glmnet(x=X, y=y)

解パス図の描画を行います.
> # Drawing Solution path
> plot(res, xvar="lambda", label=TRUE, 
+      xlab="Logarithmic value of regularization parameter", 
+      ylab="Regression coefficient", col="black", lwd=2.5) 

正則化パラーメータの値を 20 に固定して,係数の推定値を求めます.
> # Fix value of regularization parameter to 20
> res1 <- glmnet(x=X, y=y, lambda=20)
> res1$beta  # Estimated value of coefficient
5 x 1 sparse Matrix of class "dgCMatrix"
          s0
V3 133.50551
V4 -25.22804
V5  19.45576
V6   .      
V7   .      

上記の例のJupyter Notebookファイルは,GitHubStatical Modeling via Sparse Estimation.ipynbの"Apply Lasso on crime data"で見ることができます.

..-. --- --- - -.--- - . * crime.txt の中身
478 184 40 74 11 31 20
494 213 32 72 11 43 18
643 347 57 70 18 16 16
341 565 31 71 11 25 19
773 327 67 72 9 29 24
603 260 25 68 8 32 15
484 325 34 68 12 24 14
546 102 33 62 13 28 11
424 38 36 69 7 25 12
548 226 31 66 9 58 15
506 137 35 60 13 21 9
819 369 30 81 4 77 36
541 109 44 66 9 37 12
491 809 32 67 11 37 16
514 29 30 65 12 35 11
371 245 16 64 10 42 14
457 118 29 64 12 21 10
437 148 36 62 7 81 27
570 387 30 59 15 31 16
432 98 23 56 15 50 15
619 608 33 46 22 24 8
357 218 35 54 14 27 13
623 254 38 54 20 22 11
547 697 44 45 26 18 8
792 827 28 57 12 23 11
799 693 35 57 9 60 18
439 448 31 61 19 14 12
867 942 39 52 17 31 10
912 1017 27 44 21 24 9
462 216 36 43 18 23 8
859 673 38 48 19 22 10
805 989 46 57 14 25 12
652 630 29 47 19 25 9
776 404 32 50 19 21 9
919 692 39 48 16 32 11
732 1517 44 49 13 31 14
657 879 33 72 13 13 22
1419 631 43 59 14 21 13
989 1375 22 49 9 46 13
821 1139 30 54 13 27 12
1740 3545 86 62 22 18 15
815 706 30 47 17 39 11
760 451 32 45 34 15 10
936 433 43 48 26 23 12
863 601 20 69 23 7 12
783 1024 55 42 23 23 11
715 457 44 49 18 30 12
1504 1441 37 57 15 35 13
1324 1022 82 72 22 15 16

940 1244 66 67 26 18 16

Crime data

Crime data for 50 U.S cities.
The data (X1, X2, X3, X4, X5, X6, X7) are for each city.
X1 = total overall reported crime rate per 1 million residents
X2 = reported violent crime rate per 100,000 residents
X3 = annual police funding in $/resident
X4 = % of people 25 years+ with 4 yrs. of high school
X5 = % of 16 to 19 year-olds not in highschool and not highschool graduates.
X6 = % of 18 to 24 year-olds in college

X7 = % of people 25 years+ with at least 4 years of college


2019年7月23日火曜日

[R] glmnetパッケージをインストール

ターミナルからRを立ち上げます.
$ r

R version 3.6.0 (2019-04-26) -- "Planting of a Tree"
Copyright (C) 2019 The R Foundation for Statistical Computing
Platform: x86_64-apple-darwin17.7.0 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

During startup - Warning messages:
1: Setting LC_COLLATE failed, using "C" 
2: Setting LC_TIME failed, using "C" 
3: Setting LC_MESSAGES failed, using "C" 
4: Setting LC_MONETARY failed, using "C" 

glmnetがインストールされてない状態でglmnetを読み込もうとしても,以下のようにそのようなパッケージが存在しないと返ってきます.
> library(glmnet)
Error in library(glmnet) : there is no package called ‘glmnet’

インストールは以下のようにして行います.
> install.packages("glmnet")

すると,以下のようにCRANミラーを選択するよう支持されます.
Installing package into ‘/usr/local/lib/R/3.6/site-library’
(as ‘lib’ is unspecified)
--- Please select a CRAN mirror for use in this session ---
Secure CRAN mirrors 

 1: 0-Cloud [https]                   2: Algeria [https]                
 3: Australia (Canberra) [https]      4: Australia (Melbourne 1) [https]
 5: Australia (Melbourne 2) [https]   6: Australia (Perth) [https]      
 7: Austria [https]                   8: Belgium (Ghent) [https]        
 9: Brazil (PR) [https]              10: Brazil (RJ) [https]            
11: Brazil (SP 1) [https]            12: Brazil (SP 2) [https]          
13: Bulgaria [https]                 14: Chile [https]                  
15: China (Hong Kong) [https]        16: China (Guangzhou) [https]      
17: China (Lanzhou) [https]          18: China (Shanghai) [https]       
19: Colombia (Cali) [https]          20: Czech Republic [https]         
21: Denmark [https]                  22: Ecuador (Cuenca) [https]       
23: Ecuador (Quito) [https]          24: Estonia [https]                
25: France (Lyon 1) [https]          26: France (Lyon 2) [https]        
27: France (Marseille) [https]       28: France (Montpellier) [https]   
29: France (Paris 2) [https]         30: Germany (Erlangen) [https]     
31: Germany (Göttingen) [https]      32: Germany (Münster) [https]      
33: Germany (Regensburg) [https]     34: Greece [https]                 
35: Hungary [https]                  36: Iceland [https]                
37: Indonesia (Jakarta) [https]      38: Ireland [https]                
39: Italy (Padua) [https]            40: Japan (Tokyo) [https]          
41: Japan (Yonezawa) [https]         42: Korea (Busan) [https]          
43: Korea (Gyeongsan-si) [https]     44: Korea (Seoul 1) [https]        
45: Korea (Ulsan) [https]            46: Malaysia [https]               
47: Mexico (Mexico City) [https]     48: Norway [https]                 
49: Philippines [https]              50: Serbia [https]                 
51: Spain (A Coruña) [https]         52: Spain (Madrid) [https]         
53: Sweden [https]                   54: Switzerland [https]            
55: Turkey (Denizli) [https]         56: Turkey (Mersin) [https]        
57: UK (Bristol) [https]             58: UK (London 1) [https]          
59: USA (CA 1) [https]               60: USA (IA) [https]               
61: USA (KS) [https]                 62: USA (MI 1) [https]             
63: USA (MI 2) [https]               64: USA (OR) [https]               
65: USA (TN) [https]                 66: USA (TX 1) [https]             
67: Uruguay [https]                  68: (other mirrors)                

ここでは,40番のJapan(Tokyo)を選択します.
Selection: 40

するとインストール*が始まります.

・・-・ --- --- - -・ --- - ・
* インストール時には以下のように表示されます.
also installing the dependencies ‘iterators’, ‘foreach’

trying URL 'https://cran.ism.ac.jp/src/contrib/iterators_1.0.10.tar.gz'
Content type 'application/x-gzip' length 290575 bytes (283 KB)
==================================================
downloaded 283 KB

trying URL 'https://cran.ism.ac.jp/src/contrib/foreach_1.4.4.tar.gz'
Content type 'application/x-gzip' length 360705 bytes (352 KB)
==================================================
downloaded 352 KB

trying URL 'https://cran.ism.ac.jp/src/contrib/glmnet_2.0-18.tar.gz'
Content type 'application/x-gzip' length 3862714 bytes (3.7 MB)
==================================================
downloaded 3.7 MB

During startup - Warning messages:
1: Setting LC_TIME failed, using "C" 
2: Setting LC_MESSAGES failed, using "C" 
3: Setting LC_MONETARY failed, using "C" 
* installing *source* package ‘iterators’ ...
** package ‘iterators’ successfully unpacked and MD5 sums checked
** using staged installation
** R
** inst
** byte-compile and prepare package for lazy loading
During startup - Warning messages:
1: Setting LC_TIME failed, using "C" 
2: Setting LC_MESSAGES failed, using "C" 
3: Setting LC_MONETARY failed, using "C" 
** help
*** installing help indices
** building package indices
** installing vignettes
** testing if installed package can be loaded from temporary location
During startup - Warning messages:
1: Setting LC_TIME failed, using "C" 
2: Setting LC_MESSAGES failed, using "C" 
3: Setting LC_MONETARY failed, using "C" 
** testing if installed package can be loaded from final location
During startup - Warning messages:
1: Setting LC_TIME failed, using "C" 
2: Setting LC_MESSAGES failed, using "C" 
3: Setting LC_MONETARY failed, using "C" 
** testing if installed package keeps a record of temporary installation path
* DONE (iterators)
During startup - Warning messages:
1: Setting LC_TIME failed, using "C" 
2: Setting LC_MESSAGES failed, using "C" 
3: Setting LC_MONETARY failed, using "C" 
* installing *source* package ‘foreach’ ...
** package ‘foreach’ successfully unpacked and MD5 sums checked
** using staged installation
** R
** demo
** inst
** byte-compile and prepare package for lazy loading
During startup - Warning messages:
1: Setting LC_TIME failed, using "C" 
2: Setting LC_MESSAGES failed, using "C" 
3: Setting LC_MONETARY failed, using "C" 
** help
*** installing help indices
** building package indices
** installing vignettes
** testing if installed package can be loaded from temporary location
During startup - Warning messages:
1: Setting LC_TIME failed, using "C" 
2: Setting LC_MESSAGES failed, using "C" 
3: Setting LC_MONETARY failed, using "C" 
** testing if installed package can be loaded from final location
During startup - Warning messages:
1: Setting LC_TIME failed, using "C" 
2: Setting LC_MESSAGES failed, using "C" 
3: Setting LC_MONETARY failed, using "C" 
** testing if installed package keeps a record of temporary installation path
* DONE (foreach)
During startup - Warning messages:
1: Setting LC_TIME failed, using "C" 
2: Setting LC_MESSAGES failed, using "C" 
3: Setting LC_MONETARY failed, using "C" 
* installing *source* package ‘glmnet’ ...
** package ‘glmnet’ successfully unpacked and MD5 sums checked
** using staged installation
** libs
gfortran  -fPIC  -g -O2  -c glmnet5dpclean.f -o glmnet5dpclean.o
clang -I"/usr/local/Cellar/r/3.6.0_3/lib/R/include" -DNDEBUG   -I/usr/local/opt/gettext/include -I/usr/local/opt/readline/include -I/usr/local/include  -fPIC  -g -O2  -c glmnet_init.c -o glmnet_init.o
clang -dynamiclib -Wl,-headerpad_max_install_names -undefined dynamic_lookup -single_module -multiply_defined suppress -L/usr/local/Cellar/r/3.6.0_3/lib/R/lib -L/usr/local/opt/gettext/lib -L/usr/local/opt/readline/lib -L/usr/local/lib -o glmnet.so glmnet5dpclean.o glmnet_init.o -L/usr/local/opt/gcc/lib/gcc/9/gcc/x86_64-apple-darwin17/9.1.0 -L/usr/local/opt/gcc/lib/gcc/9 -lgfortran -lquadmath -lm -L/usr/local/Cellar/r/3.6.0_3/lib/R/lib -lR -lintl -Wl,-framework -Wl,CoreFoundation
ld: warning: text-based stub file /System/Library/Frameworks//CoreFoundation.framework/CoreFoundation.tbd and library file /System/Library/Frameworks//CoreFoundation.framework/CoreFoundation are out of sync. Falling back to library file for linking.
installing to /usr/local/lib/R/3.6/site-library/00LOCK-glmnet/00new/glmnet/libs
** R
** data
** inst
** byte-compile and prepare package for lazy loading
During startup - Warning messages:
1: Setting LC_TIME failed, using "C" 
2: Setting LC_MESSAGES failed, using "C" 
3: Setting LC_MONETARY failed, using "C" 
** help
*** installing help indices
** building package indices
** installing vignettes
** testing if installed package can be loaded from temporary location
During startup - Warning messages:
1: Setting LC_TIME failed, using "C" 
2: Setting LC_MESSAGES failed, using "C" 
3: Setting LC_MONETARY failed, using "C" 
** checking absolute paths in shared objects and dynamic libraries
** testing if installed package can be loaded from final location
During startup - Warning messages:
1: Setting LC_TIME failed, using "C" 
2: Setting LC_MESSAGES failed, using "C" 
3: Setting LC_MONETARY failed, using "C" 
** testing if installed package keeps a record of temporary installation path
* DONE (glmnet)

The downloaded source packages are in
‘/private/var/folders/sd/hzc6sg0j5tj5p340n4h992_r0000gn/T/RtmphxoZOU/downloaded_packages’

2019年7月1日月曜日

[gnuplot] gnuplotのインストール (with Aquaterm)

macOSにHomebrew経由でgnuplotをインストールすると,terminal type が 'qt' のみとなってしまい,非常に困ってしまいます.
以前は,Homebrewのインストールで
$ brew install gnuplot --with-aquaterm
とすれば良かったのですが,現段階では何故かエラーとなってしまい,AquaTermが一緒にインストールすることができません.

AquaTermがないと.epsファイルへの書き出しができずに困ってしまうので,以下の方法でインストールを行いました.

まずは,Homebrew経由でインストールしたgnuplotを削除します.

Homebrewでインストールしたものを検索してみます).
$ brew list
gnuplot
(実際には,gnuplot以外にも色々と出てくると思うが,上記の例では,gnuplotのみを表示しています.

gnuplotがあることを確認して,以下のコマンドで削除します.
$ brew uninstall gnuplot

続いて,gnuplotのホームページからソースコードを入手します(この段階では ver. 5.2 を選択したので,ダウンロードしたファイルは gnuplot-5.2.7.tar です).
gnuplot-5.2.7.tar を展開して,適当なディレクトリに移動します(以下の例では xxx というディレクトリに移動したとして説明します).

gnuplot-5.2.7.tar を展開して移動したディレクトリ(xxx)に移ります
$ cd /Users/xxx/gnuplot-5.2.7 

以下のコマンドを入力します.
$ ./configure --with-readline=builtin --with-aquaterm

このコマンドを入力すると,以下のように相当ズラズラと流れます(以下の例では途中を省略しています).
checking for a BSD-compatible install... /usr/bin/install -c
checking whether build environment is sane... yes
checking for a thread-safe mkdir -p... ./install-sh -c -d
checking for gawk... no
checking for mawk... no
checking for nawk... no
checking for awk... awk
... <- Omission
gnuplot will install the following additional materials:

  cfg file for epslatex terminal: yes
  TeX *.sty for lua/tikz terminal: yes
  TeX files will be installed in /usr/local/texlive/texmf-local/tex/latex/gnuplot
                               (use --with-texdir=DIR to change)
  Help file: ${datarootdir}/gnuplot/5.2/gnuplot.gih
  PostScript prologue files: ${datarootdir}/gnuplot/5.2/PostScript/

続いて make します.
$ make

make においても,以下のように相当ズラズラと流れます(以下の例では途中を省略しています).
/Applications/Xcode.app/Contents/Developer/usr/bin/make  all-recursive
Making all in config
... <- Omission
make[3]: Nothing to be done for `all'.

cp -p ./Gnuplot.app-defaults Gnuplot

make[2]: Nothing to be done for `all-am'.

続いて,以下のコマンドを実行します.この段階ではパスワードの入力を求められるので,administrator としてログインする際のパスワードを入力します.
mini:gnuplot-5.2.7 hide$ sudo make install
Password:

この際にも以下のように,多少ズラズラと流れてインストールが完了します(例では途中を省略しています).
aking install in config
make[2]: Nothing to be done for `install-exec-am'.
make[2]: Nothing to be done for `install-data-am'.
... <- Omission
make[3]: Nothing to be done for `install-exec-am'.
 .././install-sh -c -d '/usr/local/share/gnuplot/5.2/app-defaults'
 /usr/bin/install -c -m 644 Gnuplot '/usr/local/share/gnuplot/5.2/app-defaults'
 .././install-sh -c -d '/usr/local/share/gnuplot/5.2'
 /usr/bin/install -c -m 644 colors_default.gp colors_podo.gp colors_mono.gp gnuplotrc '/usr/local/share/gnuplot/5.2'
make[2]: Nothing to be done for `install-exec-am'.
make[2]: Nothing to be done for `install-data-am'.

確認のためにgnuplot を立ち上げてみます.
$ gnuplot

G N U P L O T
Version 5.2 patchlevel 7    last modified 2019-05-29 

Copyright (C) 1986-1993, 1998, 2004, 2007-2018
Thomas Williams, Colin Kelley and many others

gnuplot home:     http://www.gnuplot.info
faq, bugs, etc:   type "help FAQ"
immediate help:   type "help"  (plot window: hit 'h')

Terminal type is now 'aqua'
gnuplot> 
 
上記のように
Terminal type is now 'aqua'
と表示されていれば成功です.

macOS High Sierra は上記の方法でインストールに成功しましたが,Majaveでは駄目でした...