The Nameless City

何故か製薬やSAS関連のブログ、の予定。

倍精度浮動小数点型の限界はあるよ、という話。

データステップ100万回      SAS新手一生: マクロ変数に数値をいれて戻すと誤差がでちゃう場合がある問題について考える話を見て。


しばしば、お客さんに説明しても理解してもらえないので超絶面倒臭い話の周りになるんですがね。


確かぼんやりとした記憶ですが、浜田先生のSASの本で、対応ありのt検定をするのにunivariateプロシージャを使っていたのですが(多分SAS6ぐらいの世界の話です)、そこでround関数を使っていたかと思います。
ちなみに、SASV9V6ではround関数の実装にバグがあって、%mround使いましょうねーというお話もありました。

倍精度浮動小数点型の限界。

data _null_ ;
	a = 1.23456789012345678 ;
	put a best32. ;
run ;
                1.23456789012345

SASのマニュアルにも記述はありますが、IEEEで定められた倍精度浮動小数点型の場合、表現に限界があります。
有効数字だいたい15桁ぐらいしか表現出来ません。十分という話もあるんですが。


もうひとつは、この倍精度小数点型、内部表現(計算機内部で保持されている形式)と表示が、1:1の対応にならないんですね。
特に精度が問題になってくる時はそうです。
IEEE 754 - Wikipedia



数値でデフォルトで使われるフォーマットがbest12.なのもその辺りが理由かと思います。

故に。

2進数表現だと、完全に一致させる事は出来ます。
だから何だとか言われても困りますが。


もうひとつ。
基準とするような値、例えば0.005などと大小比較する際には、等号を使わないなど考えて下さい。引き算してfuzz通すとかして下さい。
古いプロシージャでは計算機イプシロンとかが1e-8のまんまだったりするのですが、オプションで変えられるものもあります。だいたい1e-12辺りが他のプロシージャの精度とも合います。
また、情報落ち、桁落ちは、基本的に数字が小さい方に出ます。
で、オマケですが、この辺り、SAS9.4のヘルプをチマチマ確認してたのですが、バグが残ってます。

892  data _null_ ;
893      a = 0.5-0.4 ;
894      b = 0.4-0.3 ;
895      put 'a  :' a ' b  :' b ;
896      put 'a  :' a hex16. ' b  :' b hex16. ;
897      if a = b then put 'OK' ;
898      else put 'ERROR' ;
899      a2 = round(a,1e-12) ;
900      b2 = round(b,1e-12) ;
901      put 'a2 :' a2 ' b2 :' b2 ;
902      put 'a2 :' a2 hex16. ' b2 :' b2 hex16. ;
903      if a2 = b2 then put 'OK' ;
904      else put 'ERROR' ;
905  run ;

a  :0.1  b  :0.1
a  :3FB9999999999998 b  :3FB999999999999C
ERROR
a2 :0.1  b2 :0.1
a2 :3FB999999999999A b2 :3FB999999999999A
OK
NOTE: DATAステートメント処理(合計処理時間):
      処理時間           0.00 秒
      CPU時間            0.00 秒

ちなみに。

統計プロシージャでは、内部的に数字をどう扱っているかというと、一部十進計算していたりするようです。私が、実行時の機械語を見た訳ではなく、あくまでも伝聞なのですが。
dataステップで要約統計量を出すと、univariateプロシージャと計算が合わなかったりするのは、それ故に仕方ないです。