optimってあんまり信用できないなぁ、って話
タイトルの通りです。
仕事で使うことがあったのでRのoptimを使って回帰を解いてみたのですが、これが意外に安定しません。変数の数なのか、ダミー変数が含まれるとダメなのか、原因についてはよくわかりませんが想定以上に解がバラついてしまいました。
実行したコードは以下の通りです。まずは以下のようにデータとoptimに渡す目的関数を定義します：
## サンプルデータとしてirisを使用 dat <- iris datMdl <- as.data.frame(model.matrix(~.-1, data = dat)) ## Speciesを含まないデータ dat01 <- datMdl[, 1:4] ## Speceisを含むデータ（virginicaを除外） dat02 <- datMdl[, -7] ## 最小二乗法を以下のように定義 myLM <- function(b, dat) { Y <- as.matrix(dat[, 1]) X <- as.matrix(cbind(1, dat[, -1])) b <- as.matrix(b) prd <- X %*% b e <- Y - prd return(t(e) %*% e) }
以上の条件で、まずはSpeciesを含まないデータを用いてlmとoptimをそれぞれ実行します。
## Speciesを含まないデータでSepal.Lengthをyとして回帰をかける resLM_01 <- lm(Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width, data=dat01) ## optimによる回帰係数の取得 init_coef <- rep(1, 4) resOpt_01 <- optim(par=init_coef, fn=myLM, dat=dat01[, 1:4]) ## 係数の確認。だいたい同じ cbind( summary(resLM_01)$coef[, 1], resOpt_01$par) [,1] [,2] (Intercept) 1.8559975 1.8589403 Sepal.Width 0.6508372 0.6501880 Petal.Length 0.7091320 0.7087105 Petal.Width -0.5564827 -0.5559860
コメントでも書いてありますが、係数は概ね一致します。続けてSpeciesを含めたデータによりlmおよびoptimをそれぞれ再度実行します。
## 次にSpeciesを含むデータで同様に回帰を実行 resLM_02 <- lm(Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width + Speciessetosa + Speciesversicolor, data=dat02) ## 同じく、optimによる回帰係数の取得 init_coef <- rep(1, 6) resOpt_02 <- optim(par=init_coef, fn=myLM, dat=dat02) ## 係数の確認。結構違う cbind( summary(resLM_02)$coef[, 1], resOpt_02$par) [,1] [,2] (Intercept) 1.1477685 0.6123026 Sepal.Width 0.4958889 0.4621860 Petal.Length 0.8292439 0.8993939 Petal.Width -0.3151552 -0.1967130 Speciessetosa 1.0234978 1.5512063 Speciesversicolor 0.2999359 0.4695160
Rの最適化ソルバーはoptim以外にもありますし(optimizeとか)、optimを使うにしてもパラメータや用いるアルゴリズム(Nelder-Mead、BFGSとか)など設定できる箇所は色々とあります。なのでこの結果だけを元にoptimを信用できないとは言えないのですが、このくらいのデータ・モデルでこれほどlmとの差が出てしまうと複雑な問題への適用はちょっとためらってしまいますね。
なお当初の目的としてはMarketing Mix Modelによる広告効果の推定を行う際に、残存効果の減衰率をoptimで推定してみようというものでした。optimの結果があまり芳しくないためこういった検証を行っていたのですが、この辺はいずれ記事にまとめようと思います。それからoptimはオプションでHessian(ヘッセ行列)を出力することにより標準誤差を計算できますので、その辺もいずれ書いてみようと思います。
リッジ回帰の解の特性
SAPのソリューションに「InfiniteInsight」というものがある。旧KXENをSAPが買収したもので、搭載されているアルゴリズムはサポートベクターマシンの開発者（VapniK）によるものらしく、統計解析に対する知識や経験なしに高精度のモデルを作成することができるらしい。
さて、このInfiniteInsightでは「分類・回帰」に類するものに対する分析では「リッジ回帰」が適用されるという話を耳にした。事の真偽はさておき、この「リッジ回帰」に思うことを述べてみる。
まずリッジ回帰とは何かと問われると、「L2ノルムの二乗を正則化項（ペナルティ）として与える罰則付き最小二乗法」である。いわゆる縮小推定法の一種で、学習データへの過適合を防止するために用いられたりするのだが、これはデータに対して過適合を生じたり多重共線性が問題となるような状況では一般にパラメータの絶対値が大きくなりやすいため、パラメータスペースに制約を設けることが可能なためである（ちなみにL1ノルムを罰則として与える方法をLasso回帰と呼ぶ）。
ここで思うのだが、リッジ回帰（Lassoも）を適用しようと考えている人は「汎化誤差を小さくする可能性が高まる代償として、不偏性や有効性を失っている」ことにまで考えが及んでいるのだろうか？
誤差が正規分布にしたがう変数に対して最小二乗法により解を推定した場合、この解はBLUEとなり、不偏性や有効性・一致性を備えることが知られているが、リッジ回帰ではそこに正則化項を加えることでパラメータが極端な値を取ることを防いでいる。もちろんこの解は不偏推定量ではないため、解の期待値が真の解と一致しないことを容認していると言える。そこまで踏まえ、それでもなおリッジ回帰を選択したというならば良いのだが、これまでに見てきたケースでは多くの場合、そこまで考えたものではなかった。交差検定による検証を行った結果、テストデータに対する予測性能が高くなることから「汎化性能が高い」→「良いモデル」という安直な考えによるものであったと思う。
近い話では、昨今の機械学習の流行に対しても似たような思いを持っている。Kaggleのようなデータ解析コンペでは様々なアルゴリズムを組み合わせた予測モデルを作ることがトレンドであるようだが、そのような方法で求められた解には一体どのような意味があるのだろう？
統計学における「モデル」という言葉の意味をもう一度よく考えてみたい。
できないと言ってはいけません
仕事をしていて思ったことを。
統計分析が求められる場面というものはいくつもあるが、その中で最もメジャーと言ってよいものの一つが「将来予測」だと思う。身近なもので言えば株価や天気で、「明日はどうなるか？」という問いについて、モデルは一つの答えを出してくれる。もちろん、その精度は検討した変数や手法によるところが大きいだろうが、単に勘や経験に頼ったものよりもよっぽど良い答えを出してくれるだろう。
そんな便利な統計モデルであるが、弱点もある。例えば、将来予測を行う際にはモデルに含まれる変数がすべて既知であらねばならない、というものだ。Rでpredict(newdata)を行うときや、SPSS Modelerで作成済みのモデルノードに新規データを渡す場面を想定してほしい。この時、新規データにはモデルに含まれる変数がすべて格納されている必要がある。
「何を今さらそんな…」と思われるかもしれない。しかし、データを触り、分析にかけ、予測を行うという行為の経験がなく、実際の作業としてはどのようなものであるかイメージができない場合には意外とこの辺りの感覚が抜け落ちている。つまり、一度分析を行ってしまえば無条件で新規のデータに対しても予測が可能であると思ってしまうらしい。
具体的な例を出そう。例えば携帯電話の顧客情報をもとにセグメンテーションしたいとする。顧客のデモグラ情報（年代や性別）の他、携帯電話の使用状況を変数として用いた結果、非常によいモデルができたとしよう。それではこのモデルを用いて、新規顧客がいずれのクラスターに分類されるかを判断できるだろうか？
答えはノーだ。なぜならば、新規顧客については携帯電話の使用状況に関するデータがなく、モデルを当てはめることができないからだ。クラスタリングでなく回帰モデルで例えれば、複数ある変数の内のいくつかがNULLになっている状況である。
それでは、このようなクラスタリングに意味はないのだろうか？
当然、そんなことはない。もし情報が足りなくてクラスタリングができないのであれば、どういった情報を取得すれば良いのか考えれば良い。あるいは取引が開始するまでの接触の段階で情報収集を重ね、少しずつクラスタリングの精度を向上させていくのも良い。最初は間違っていてもよいから、何らかのアクションが打てるように分析結果を活用する、というのがコンサルに求められる態度であり、データ分析者と異なるところなのだろう。
『「できない」という言葉は何の意味もないから、「どうすればできるか」を考えよう』
以前に上司から教えられた言葉である。データ分析者としてはすっきり納得がいくものではないが、コンサルとしては非常に大事な教えであったなぁと思う今日この頃。
Data Scientist Workbench
こんなサービスがあることを最近知った。
■Data Scientist Workbench
datascientistworkbench.com
IBMが提供しているサービスで、ブラウザ上でRStudioやJupyterが実行できる環境を提供してくれる（無料で！）。分析手法よりもHadoopやSparkなど大規模データの捌き方をメインにしている様子で、用意されているサンプルプログラムもRやPythonでのSparkの実行例などが多い。
Big Data Universityのオンラインコースへのリンクもあり、基礎的な学習とR/Pythonでの実装についてはここで十分に学べると思う。
Seahorseというツールは初めて見たけどAzure MLに似ている。このUIは個人的には使いにくいと思うんだけど、他の人はどうなんだろうか。。。
ちなみにBig Data Universityはこちら
Analytics, Big Data, and Data Science Courses - Big Data University
RでBI
忙しかったので久しぶりの更新です。
最近Rでこんなパッケージを発見して驚いた。その名もrpivotTable。なんとRでTableauのようなBIツールを再現してしまうという夢のようなパッケージだ。パッケージのインストールから実行までは以下のようになる：
install.packages("rpivotTable") library(rpivotTable) data(iris) rpivotTable(iris, rows="Species", "100%", "40pt")
上記を実行するとブラウザが開き、Tableauで見慣れたようなピボットテーブルが現れる。あとは好きなように変数を配置して表計算やグラフを作成すればよい。
もちろんTableauの方がグリグリ感は強いし、rpivotTableではドリルダウンなどもできないのだけれども、可視化と分析をシームレスに実行できるという点ではよっぽど使い勝手は良いと感じた。
出来映えも大変良く、shiny以来の衝撃を感じたので紹介しておく。
遺伝的アルゴリズムで変数選択を行うパッケージ
遺伝的アルゴリズム(GA)をRで実行するためのパッケージを探していたら、GAによる変数選択という非常に興味深いものを発見した。パッケージ名はGALGO。Genetic Algorithms to solve Optimization problemsらしい。早速使ってみたが、インストールに少し手間取ったので備忘として残しておく。
初めに成功したパターンを記載する。基本的にこちらのページに書いてある通りで良い。
"R.oo"や"R.methodsS3"、その他の必要なパッケージを予めインストールしておき、ダウンロードしておいたgalgo_1.2.tar.gzをソースからインストールするだけ。
CRANではなくソースからインストールする場合は、ソースが置いてある場所までsetwdした後に
install.packages("galgo_1.2.tar.gz", repos=NULL, type="source")
とすれば良い。
では何に手間取ったかというと、「galgo r example」などでGoogle検索した際にトップに表示される下のページ：
こちらに記載されている通りに進めると、古いバージョンのgalgoがリンクされていることから.tar.gzの中にNAMESPACEファイルが存在していないため、一度解凍したのちにNAMESPACEファイルを作成して再び圧縮する必要がある。
tar -xvf galgo_1.0-10.tar.gz cd galgo echo 'exportPattern( "." ) > NAMESPACE cd .. tar -zvf galgo_1.0_10mod.tar.gz galgo
こんな感じ。
ここまで進めて、さぁインストールするぞ！と思ってinstall.packages()を実行すると、エラーの嵐でインストールできない。Rのバージョンの問題なのかもしれないが未解決。
その後、再度他のページを探した結果、最初に紹介したページに行き着いてインストールが無事にできた。使用感はいずれ。
