2015年5月27日水曜日

小数点以下の微妙な誤差について

みなさん、SASやExcel等で出した結果が小数点以下の深い桁で微妙に合わない
「なぜ?」という経験はありませんか?


Excelだと、こんな感じの現象。













SASだと、こんな感じの現象。





今回はこの現象について、なるべく簡単にお伝えできたらと思います。

深く知りたい方は、コンピュータの仕組や世界的に準拠されている
技術標準 IEEE(アイ・トリプリ・イー)辺りを調べていただくと
よいかと思います。


では、改めてこの現象が起こる理由についてですが
それはズバリ!!

コンピュータでの計算は2進数で行われているということ
そして『10進数の小数を2進数へ変換 及び 2進数から10進数へ戻す際の
丸め誤差によるズレ』 これが、理由です。


...???


と、言葉で言うとこれだけですが
まったくわかりませんよね...

なので、少しだけ深堀したいと思います。



『有効桁数』

コンピュータ上で表現できる、数値の桁数のことを「有効桁数」とよんでいます。

コンピュータ上で計算する場合、人力でやるには無理があるような桁数でも
計算が可能ではありますが、普段使用する中でそこまで大きい桁数を必要とするわけでもなく
また、たとえコンピュータとはいえ、無限に桁数を増やせるわけではありません。

そこで、有効と思われる桁数を決めて、その範囲内で計算するようになっており
これを「有効桁数」と読んでいます。

「有効桁数」は一般的には「IEEE」が定めている基準を準拠しており
小数点以下の桁数は15桁となっています。
それに伴い「有効桁数」を超えた場合、丸め処理が行われています。



『循環小数』

ある桁から先で同じ数字の列が無限に繰り返される少数のことを「循環小数」とよびます。
例えば、「0.3333333333・・・・」のことです。

この「循環小数」は桁数を決めておかなければ、、際限なく桁数を増やしてしまいます。
その桁数を決めているのが「有効桁数」になります。



『計算は2進数』

コンピュータ上では私たちが普段行う10進数での計算は行われません。
通常は2進数を用いて計算が行われます。

なぜ、2進数で行われているかについては、ここでは省略しますが
整数や小数であっても2進数に変換され、計算後10進数に再変換され
画面上に表示されています。



『誤差発生の仕組み』

コンピュータ上において「計算は2進数」で行われるため、小数であっても10進数から2進数に変換しますが、ほとんどの場合において「循環小数」になってしまい正確に表現できないという
問題が存在します。
そのため「有効桁数」にあわせるため、丸め処理が行われます。
丸め処理が行われた数値を10進数に戻すため、誤差が生まれます。


例えば、10進数で「0.1」を2進数に変換した場合、以下の図のように
「000110011001100110011....」と「循環小数」なります。
ここで、「有効桁数」を超えた部分で丸め処理を行います。











この丸め後の2進数を10進数に変換すると、「0.10000000000000000055511 …」となり
赤色の数字の部分が誤差になるというわけです。



まとめ

ザックリとではありますが、わかりましたでしょうか。
結論を言いますと、『コンピュータ上で計算する限り、計算誤差は必ずある』です。

こう言いますと、どうにもならないようにも思えますが対策することは可能です。
すぐに実践可能な対策に以下のものがあります。

・あらかじめ必要な桁数を決めておき、四捨五入や切り捨てなどを使用する。
    あらかじめ、使用する桁数を決めておけば、誤差が出る深い部分を
    使用しないで済みます。
    ほとんどの場合において小数点以下第2位、多くても第5位くらいではないでしょうか。

    また、比較する際においても、決めておいた桁数のみで比較を行い、
    それ以下の桁は許容範囲内として処理する。
    但し、今回の現象とは無関係の部分での誤差まで許容してしまわないように
    注意が必要です。


いかがでしたでしょうか、深く知っていなくても
「こんな現象が起きてるんだな~」
っていうのだけでも、覚えておいてもらえれば
実際に直面した時にスムーズに対処できるのではないでしょうか。


2015年5月11日月曜日

事務所移転のお知らせ

ネイチャーインサイトは、2015年5月11日に事業拡大のため事務所を下記の通り移転をしました。

今後も社員一同、より一層のサービス向上に励み、お客様のご期待添いますよう努めてまいる所存です。
お引き立てを賜りますよう、よろしくお願い申し上げます。



【移転先】
〒101-0052
東京都千代田区神田小川町3-3
神田小川町トーセイビルⅡ 5階
TEL:03-3518-6061
FAX:03-3518-6062


(アクセス)
都営新宿線「小川町駅」より徒歩6分
東京メトロ半蔵門線・都営三田線・都営新宿線「神保町駅」より徒歩6分
東京メトロ丸ノ内線「淡路町駅」より徒歩7分
東京メトロ千代田線「新御茶ノ水駅」より徒歩7分
JR総武線・中央線「御茶ノ水駅」御茶ノ水橋口より徒歩8分



                


新しいオフィスのエントランスです。
お近くにお越しの際は、是非お立ち寄り下さい!


2015年5月1日金曜日

【求人情報】 2015年7月入社 SASエンジニア / 生物統計解析業務



                                                                           



明日からゴールデンウィークですね。
皆さんは、どこへお出かけの予定でしょうか?

私はとくに予定がないので、運動でもして体力アップを図りたいと思います!






さて、ネイチャーインサイトでは、2015年第一弾の求人募集をスタートしました。
ゴールデンウィーク明けには事務所を移転し、新オフィスでお迎えいたします。

当社では、新たなチャレンジで一緒に成長できる仲間を求めております。
DODA求人広告(下記のリンク)より、是非ご応募をお待ちしております。


未経験からはじめる SASエンジニア


臨床試験のDM・統計解析業務(新規部署立ち上げメンバー)


皆さまのご応募をお待ちしております。

説明力のある回帰分析モデルは使えない?

・株価の動きを説明する回帰モデルを作れば株で儲けられるのでは?

今回も回帰分析の話をします。例として、日経平均株価について回帰分析をして株価予測モデルを作ってみましょう。日経平均株価は「日経平均プロフィル」の以下のページから csv ファイルをダウンロードしました
http://indexes.nikkei.co.jp/nkave/index
この csv ファイルをSASに読み込みます。ただし、以下のプログラムでは、読み込んだデータセットを半分に分けています。これは、前半のデータでモデルを推定し、後半のデータで、推定したモデルを使って擬似的に株価予測をするためです。
/* パスは適宜変更してください */
%let csvpath=C:\Users\HOGEHOGE\Downloads\nikkei_stock_average_daily_jp.csv;
libname tmp "C:\Users\HOGEHOGE\Documents";
data work.stock;
    attrib  date  length=8 label='経過日数'          
    end   length=8 label='終値'
    start length=8 label='始値'
    high  length=8 label='高値'
    low   length=8 label='安値'
    _date_ length=8 label='日付' format=yymmdds10.        ;
    infile "&csvpath." dsd firstobs=2;
    input date: yymmdd10. end start high low;
    _date_=date;
    date = date- '03Jan2012'd;
run;
data work.stock_1 work.stock_2;
    set work.stock;
    /* 2013年8月14日で分割 */
    if _N_ <= 400 then output work.stock_1;
    else                output work.stock_2; run; 
日毎の株価終値の折れ線グラフが以下のようになります。
日々の上がり下がりはあるものの、全体として右上がりの傾向があります。そこで、経過日数を説明変数にして回帰してみましょう。このようなモデルです。
左辺の yt は株価終値を、右辺の t は経過日数を表し、 εt
は誤差を表します。先ほど分割したデータのうち、前半のデータを使い、reg プロシジャで計算します。
proc reg data=work.stock_1;
    NIKKEI:model end = date;
quit; 
結果が以下です。
係数が有意で、決定係数の値も 0.70 とかなり高めです。株価終値の実際の値と回帰直線をグラフにすると、こうなります (グラフはプロシジャ実行時に自動で作成されます)。
しかし、線形回帰モデルはつまり、データの特徴を直線で表現しようというものです。直線からはずれた大きな上下の波を説明することはできません。直線と実際の株価とで1000円単位の誤差があるのでは、実用に問題がありますね。

・多項式モデルなら株価の動きをかなり正確に説明できる

では、直線でなく曲線で回帰してみましょう。つまり、
というふうに経過日数の多項式にすれば上下に波のある株価の動きを表現できるはずです。 t2, t3, ... と続く t のべき乗を、t とは個別の説明変数と考えれば、これは実質的に重回帰と同じなので、最小二乗法で計算できます。まずは3次式
で回帰してみましょう。ここで、毎回コードを書くのは面倒なので、任意の次数の多項式で回帰するためのマクロを用意します。
%macro poly_reg(d=2);
    /* べき乗を計算*/
    data work.stock_1
         work.stock_2;
        set work.stock;
        %do i=1 %to &d.;
            date_&i. = date**&i.;
        %end;
        if _N_ <= 400 then output work.stock_1;
        else                output work.stock_2;
    run;
    /* 最小二乗法で推定*/
    proc reg data=work.stock_1
             outest=work.result_est;
        poly_&d.: model end = date_1-date_&d. ;
        output out=tmp.stock_1_est p=predict;
    quit;
    /* 外挿 */
    proc score data=work.stock_2
         score=work.result_est
         out  =work.temp(rename=(poly_&d.=predict))
         type=parms;
        var date_1-date_&d.;
    run;
    data tmp.stock_1_est_&d.;
        attrib  predict label='終値の予測値';
        set tmp.stock_1_est(keep=_date_ end predict);
    run;
    data tmp.stock_est_&d.;
        attrib  predict label='終値の予測値';
        set tmp.stock_1_est_&d.
            work.temp(keep=_date_ end predict);
    run;
%mend; 
このマクロではまず、次数に応じて経過日数の日付のべき乗の変数を作り、 reg プロシジャで回帰し、さらに score プロシジャを用いて、回帰に使わなかったほうのデータに対しても、推定結果から予測値を計算しています。大雑把ですが score プロシジャの構文の解説をすると,
proc score data=予測値を計算したいデータセット
           score=推定した結果(係数)の入っているデータセット
           out  =予測値の結果出力先データセット
           type =parms;
    var モデルの説明変数名の列挙;
run;
というふうになります。よって今回は data= に期間外のデータ work.stock_2 を、 score= には直前の reg プロシジャの推定結果を与えています。 type= score= で指定した入力データセットが何なのかを指定しています. 今回は係数 (パラメータ) を格納したデータセットなので、type=parms と指定しています。このように、推定したモデルを期間外 ―つまりたいていは分析の後に得られた未来のデータ― に対して適用し、予測値を計算することを外挿または補外 (extrapolation) といいます。 予測値の計算だけをするので、外挿に必要なのは右辺の説明変数だけです。
ちなみに、今回は使いませんが、ロジスティック回帰をする場合なら、 logistic プロシジャ内で scoreステートメントを使うことで外挿ができるので、 score プロシジャは不要です。それ以外の回帰分析では、この score プロシジャを使います。
これで、 %poly_reg(d=); とするだけで、任意の次数の多項式で回帰分析してくれます。
重回帰分析になると、自動で当てはめプロットをしてくれないので自分でグラフを作る必要があります。
決定係数の値が大きくなっているため、「当てはまりの良さ」言い換えるなら説明力が向上しています。プロット図からも、
と、前よりも実際の株価の推移に近づいているように見えます。では、もっと多項式の次数を増やせばさらに近づくのでは無いでしょうか?
%poly_reg(d=5);
%poly_reg(d=8);
として、5次多項式と8次多項式についても回帰してみましょう。
5次多項式
8次多項式
5次、8次と次数を増やしていくほど、曲線が実際の株価の推移に近づいています。当てはまりの良さの指標である決定係数も増加しています。これは株で儲けるチャンスでは?

・過去の株価では取引できない

しかしこれで儲ける前に、予測結果についても見てみましょう。この回帰モデルは過去 (2012年1月4日から2013年8月14日) の株価データを使って作成しました。現在分かってる部分だけでなく将来の株価の動きも説明できなければ儲けることはできません。もちろん未来の株価が実際にどうなるかを知ることができません。そこで回帰分析の計算に使った期間より後のデータ (2013年8月15日から2015年4月)を使って、予測値と実績値の比較をすることで、擬似的に未来の株価を予測できているかを確認します。最初に株価のデータセットを2分割したのはこのためです。
まずは 3次式の予測結果...
当てはまっているようにも見えますが、やはり急激な動きには反応してませんし、後半の時期はだいぶ実績値から離れていますし、期間の端では曲線が逆方向へ伸びています。では、5次、8次多項式ではどうでしょうか。3次多項式より当てはまっていたので、もっとよく予測できるはず……
……ん? おかしいですね……、大きく下へ外れています。たまたま調子が悪かったのかもしれません。8次多項式のほうはどうでしょうか?
予測値が上へ大きく伸びているので、本来の株価の推移が平たく潰れてしまっています。予測としてまったく用をなしません

・木を見て森を見ず

種明かしをすると、これは過剰適合 (overfitting) という現象です。機械学習の分野では過学習 (overlearning) とも呼ばれます。推定に用いるのは限られた期間のデータなので、このデータだけにモデルをぴったり合わせようとすると、この期間にたまたま起こったノイズ的な事象にまで合わせようとしてしまい、結果としてそれ以外の期間の予測精度が悪くなります。長期間ある株価の動きのうちの一部の限られた部分だけを見て、都合のよい解釈をしてしまっている、まさに「木を見て森を見ず」という状態です。
もうひとつ、過剰適合が問題だと直感的にわかる例を挙げます。さいころを6回振って、2の目がたまたま3回出たとします。では、この事実からさいころの2の目が出る確率が1/2だと言えるでしょうか?
ましてや株価の動きは、さいころを振ってどの面が上を向くか、という現象よりはるかに複雑な現象です。このさいころの例のような明らかな誤解も、複雑すぎる問題のなかでは誤解であると気づけないことが多々あります。
今回これらのモデルは、モデルの当てはまりの良さの指標である決定係数がかなり高かったのですが、このように、決定係数が高い=予測性能が高い、ではないことが分かります。一般に、説明変数が増えるほど過剰適合が起こりやすくなります。今回の場合、むやみに次数を増やすことで説明変数が増えていました。当てはまりの良さだけを考えて回帰分析するのは非常に危険です。
では、どうすれば過剰適合を回避できるのでしょうか。いくつか方法があります。
  1. 今回のように、データを係数の推定用と外挿用に分割する。外挿結果があまりに実績値と外れていれば、そのモデルは過剰適合である、ということが分かります。この方法は実際には今回のようにグラフを確認するだけでなく、より厳密で複雑な手続きがあり、「交叉検証法(クロス・バリデーション)」と呼ばれます。
  2. 「情報量規準 (information criterion)」を用いる。これは複数のモデルの候補があって、最終的にどれが一番予測性能があるか (=過剰適合が少ないか) を選ぶときに使用できます。モデルの当てはまりの良さを表す指標に対して、説明変数の数に応じて減点して求めた指標「情報量規準 」を用います。
  3. 「調整済み R2 乗」を使う。 reg プロシジャの結果にある「調整済み R2 乗」というのは、この過剰適合に対応した決定係数で、過剰適合があってもあまり値が大きくならないよう、通常の決定係数とは計算方法を変えています。しかし今回のように、過剰適合があっても値が大きくなることもよくあるので、絶対ではありません。そのため、上記2つに比べると、これはあくまで目安です。
最後に重要なのですが、過剰適合というのは、「モデルを説明しない」変数をむやみに追加することで発生します。つまり、過剰適合のあるモデルとは、適当な数字を説明変数として詰め込んで作られたハリボテのようなものです。というわけで、今回の株価予測モデルは、実際には意味のない説明変数を使って説明したように見せかけただけですので、これで儲けることはできません (身も蓋もないですが、そもそもそんなモデルを発見しても誰にも言いません)。過剰適合を回避するとは、無意味な説明変数を切り捨て、本質的な法則性を表す説明変数だけを残すようにモデルを洗練するということです。

・まとめ

  • 分析に使えるデータは常に、全体の限られたごく一部。
  • 限られたデータだけを説明しようとすると「木を見て森を見ず」状態になる (過剰適合)。
  • 過剰適合を回避する方法がある。
  • 重要なのは、背後に隠されたメカニズムを意識してモデルを作ること

・おまけ

xkcd というサイトに掲載されている、無茶な設定の予測モデルに関する漫画です。今回の話のように、限られたデータだけしか見ないで、背後の因果関係を無視していることを笑いの種にしています。
「グラフを見てわかるように、翌月には夫が40人を超えます。ウェディングケーキをまとめ買いするといいでしょう」
「グラフを見てわかるように、翌月には夫の数が大きなマイナスになります。全員分 (マイナス個) のウェディングケーキを買えば儲かるでしょう」
上記画像は xkcd により、 クリエイティブ・コモンズ 表示 - 非営利 2.5 一般 ライセンスの下に提供されています

ネイチャーインサイト サイトリニューアル&NIBLOGの引っ越し

ご連絡が遅くなりましたが、 ネイチャーインサイトの際とがリニューアルしました。 https://www.n-insight.co.jp/ それに伴い、NIBLOGも引っ越しすることになりました。 https://www.n-insight.co.jp/niblog/ ...