6514.　enormpとは何をする関数ですか？　プログラマーA　2005/04/30 (土) 20:47
├6520.　Re: enormpとは何をする関数ですか？　青木繁伸　2005/05/01 (日) 10:59
├6516.　Re: enormpとは何をする関数ですか？　青木繁伸　2005/04/30 (土) 22:31
│└6517.　Re^2: enormpとは何をする関数ですか？　青木繁伸　2005/04/30 (土) 23:41
│　└6518.　Re^3: enormpとは何をする関数ですか？　プログラマーA　2005/05/01 (日) 01:06
│　　└6519.　Re^4: enormpとは何をする関数ですか？　プログラマーA　2005/05/01 (日) 01:31
│　　　└6521.　Re^5: enormpとは何をする関数ですか？　プログラマーA　2005/05/01 (日) 13:58
│　　　　└6523.　Re^6: enormpとは何をする関数ですか？　青木繁伸　2005/05/01 (日) 15:45
│　　　　　└6525.　Re^7: enormpとは何をする関数ですか？　プログラマーA　2005/05/01 (日) 19:15
│　　　　　　└6528.　Re^8: enormpとは何をする関数ですか？　青木繁伸　2005/05/01 (日) 21:17
│　　　　　　　└6529.　Re^9: enormpとは何をする関数ですか？　プログラマーA　2005/05/02 (月) 01:15
└6515.　Re: enormpとは何をする関数ですか？　プログラマーA　2005/04/30 (土) 20:52

6514.　enormpとは何をする関数ですか？　プログラマーA　　2005/04/30 (土) 20:47 以下のプログラムは，Paul Johnsonという人が書いたプログラムで，Rで言うとGRASSやfBasicのパッケージのサブルーチンとして使われているパッケージ内部のそのまた内部の関数です。 質問1）何をどうする関数ですか？ 質問2）行番号20のところ以降のx4のところがバグではありません？ FUNCTION enormp(x) REAL x DOUBLE PRECISION x1,x2,x3,x4,yy1,yy2 xp1 = 0.771058495001320D-04 xp2 = -0.00133733772997339D0 xp3 = 0.0323076579225834D0 xp4 = 0.0479137145607681D0 xp5 = 0.128379167095513D0 xq1 = 0.00301048631703895D0 xq2 = 0.0538971687740286D0 xq3 = 0.375795757275549D0 xr1 = -1.36864857382717D-07 xr2 = 0.564195517478974D0 xr3 = 7.21175825088309D0 xr4 = 43.1622272220567D0 xr5 = 152.989285046940D0 xr6 = 339.320816734344D0 xr7 = 451.918953711873D0 xr8 = 300.459261020162D0 xs1 = 1.0D0 xs2 = 12.7827273196294D0 xs3 = 77.0001529352295D0 xs4 = 277.585444743988D0 xs5 = 638.980264465631D0 xs6 = 931.354094850610D0 xs7 = 790.950925327898D0 xs8 = 300.459260956983D0 xt1 = 2.10144126479064D0 xt2 = 26.2370141675169D0 xt3 = 21.3688200555087D0 xt4 = 4.65807828718470D0 xt5 = 0.282094791773523D0 xu1 = 94.1537750555460D0 xu2 = 187.114811799590D0 xu3 = 99.0191814623914D0 xu4 = 18.0124575948747D0 x3 = 0.564189583547756D0 x1 = abs(x) IF (x1.GT.0.5D0) GO TO 10 x4 = x*x yy1 = ((((xp1*x4+xp2)*x4+xp3)*x4+xp4)*x4+xp5) +1.0D0 yy2 = (((xq1*x4+xq2)*x4+xq3)*x4) + 1.0D0 enormp = x* (yy1/yy2) RETURN 10 IF (x1.GT.4.0D0) GO TO 20 yy1 = ((((((xr1*x1+xr2)*x1+xr3)*x1+xr4)*x1+xr5)*x1 + + xr6)*x1+xr7)*x1+xr8 yy2 = ((((((xs1*x1+xs2)*x1+xs3)*x1+xs4)*x1+xs5)*x1 + + xs6)*x1+xs7)*x1+xs8 enormp = 1.0D0-exp(-x*x)*yy1/yy2 IF (x .LT. 0.0D0) enormp = -enormp RETURN 20 x2 = x*x x4=1.0D0*x4 yy1 = ((((xt1*x4+xt2)*x4+xt3)*x4+xt4)*x4) + xt5 yy2 = ((((xu1*x4+xu2)*x4+xu3)*x4+xu4)*x4) + 1.0D0 enormp = (x3/x1)-(yy1*x1)/(x2*yy2) enormp = 1.0D0-exp(-x2)*enormp IF (x .LT. 0.0D0) enormp = -enormp RETURN END

　　　　　[このページのトップへ]

6520.　Re: enormpとは何をする関数ですか？　青木繁伸　　2005/05/01 (日) 10:59 なあんだ。 誤差関数 error function ですね。 |x| > 4 のときは密度関数が0になるので文番号20以下でも問題はありません（プログラム的には，だめだめプログラムですが）。 R 的に名前を付けて有れば(perfとか)分かりやすかったのにね。

　　　　　[このページのトップへ]

6516.　Re: enormpとは何をする関数ですか？　青木繁伸　　2005/04/30 (土) 22:31 > 質問1）何をどうする関数ですか？ 何なのでしょうね。 > 質問2）行番号20のところ以降のx4のところがバグではありません？ FORTRANには行番号というのはなく，あれは文番号ですね。 というのは，さておき，X4 のところって，具体的にはどれ？ コンパイラーエラーが起きるのですか，それとも，実行時のエラーですか。 バグというのとエラーは大きな違いがありますね。 バグは，一応動くけど正しく動かないというものでしょう。 私がやってみたところでは（Mac OS 10.4）ちゃんとコンパイルできて，呼んだら値は返すようですが？ あなたは，どのような環境でこれをコンパイルしたのでしょうか？

　　　　　[このページのトップへ]

6517.　Re^2: enormpとは何をする関数ですか？　青木繁伸　　2005/04/30 (土) 23:41 x4=1.0D0*x4 ですか？ 意味ないじゃないかということでしょうが， internet で検索して出てくる C言語版で， x4 *= 1.; /* huh, what? See original FORTRAN */ などとコメントされているのではあるが。 間違いではないというか， x4=1.0D0+x4 とかじゃないかといっても，この関数が何を計算しているのかが分からないと確認のしようがないですね。

　　　　　[このページのトップへ]

6518.　Re^3: enormpとは何をする関数ですか？　プログラマーA　　2005/05/01 (日) 01:06 おっしゃる通りです。お答えありがとうございます。感謝します。 そうです。それを移殖した人が，文番号20に飛んだ際に，x4を前もって定義する数値がないので，0.5というふうにそのコメントをしてある後で再定義しているのです。 試 しに，これと同じものをグラフにー5から＋5のxの値に対してプロットしてみると，どうやら−1から1の間を動く関数で，このプログラムのままでは，−4 以下のところが1に戻ってしまっているのですが。私の知識では，この関数がやっていることは，謎なのですが。そして，ご指摘のコメントをしている移殖プロ グラムもx4を0.5としているところも。

　　　　　[このページのトップへ]

6519.　Re^4: enormpとは何をする関数ですか？　プログラマーA　　2005/05/01 (日) 01:31 失礼，このプログラムでは，グラフは1に戻ることはなく，単調にー1から1までを滑らかに増加する関数になっているようです。グラフに問題はありません。 それを変形して，データxからその平均を引いて標準偏差で割ったものを使って，ルート2および2で割ったもの全体に0.5を足したものを0から1までを動くFx（おそらく標準正規分布か？）のような計算をしているのですが。 該当するそのサブルーチンを使ったプログラム箇所は，XBARという平均およびSDXという標準偏差が与えられているとき， X(I)=(X(I)-XBAR)/SDX FX(I)=.5+(ENORMP(X(I)/SQRT(2.0))/2.0) となります。このFXとは標準正規分布関数なのでしょうか？それとも何か他のものなのでしょうか？具体的には，そういう質問です。

　　　　　[このページのトップへ]

6521.　Re^5: enormpとは何をする関数ですか？　プログラマーA　　2005/05/01 (日) 13:58 僭越ながら，その他の掲示板（いちご　掲載版）でも聞きましたところ，お答え願えた方がいまして，やはり， 0.5+ENORMP(X/SQRT(2.0))/2.0 は， 標準正規分布のCDFの計算方法の1つであるということが判明しました。まだ，x4の数値のところは解決しておりませんが，enormpなるサブルーチン が，上の式とあいまって，標準正規分布のCDFであるということがわかれば，別の計算方法もあるので，これで相違なければ一見落着となります。 連休中，貴重なお時間ありがとうございました。

　　　　　[このページのトップへ]

6523.　Re^6: enormpとは何をする関数ですか？　青木繁伸　　2005/05/01 (日) 15:45 > 僭越ながら， 僭越という言葉の使い方を間違えていますよ。。たぶん 標準正規分布の累積密度関数とは違います。

　　　　　[このページのトップへ]

6525.　Re^7: enormpとは何をする関数ですか？　プログラマーA　　2005/05/01 (日) 19:15 お言葉を返すようですが，すでにその方のアドバイスにしたがい， 0.5+ENORMP(X/SQRT(2.0))/2.0 のプロットと，別の方法で標準正規分布のCDFをプロットしたものとは一致していることを確認しました。ためしに，両者の差を独自にプロットしても，10 のマイナス-16〜-15乗以下のオーダーで0に一致することが確認されました。 「たぶん」ではなくて，実際に「標準正規分布の累積密度関数」です。どうも長らくお手数とお時間をおかけしてしまいました。感謝いたします。ありがとうございました。

　　　　　[このページのトップへ]

6528.　Re^8: enormpとは何をする関数ですか？　青木繁伸　　2005/05/01 (日) 21:17 なるほど。 fOptions というライブラリには erf という関数が定義されています。pnorm を下請けに使うように書かれていますが，表現形が違うといえばそういうことになるんでしょう。 > library(fOptions) 要求されたパッケージ fBasics をロード中です Rmetrics, (C) 1999-2005, Diethelm Wuertz, GPL fBasics: Markets, Basic Statistics, Date and Time 要求されたパッケージ fCalendar をロード中です Rmetrics, (C) 1999-2005, Diethelm Wuertz, GPL fCalendar: Markets, Basic Statistics, Date and Time 要求されたパッケージ fSeries をロード中です Rmetrics, (C) 1999-2005, Diethelm Wuertz, GPL fSeries: The Dynamical Process Behind Financial Markets Rmetrics, (C) 1999-2004, Diethelm Wuertz, GPL fOptions: Valuation of Options > erf(1.96/sqrt(2))/2 [1] 0.9750021 > erf function (x) { 2 * pnorm(sqrt(2) * x) } ということですので， > pnorm(1.96) [1] 0.9750021 ですね。 このサイトのいくつかの所に，正規分布の上側確率を与えるアルゴリズムの記述が有ります。 AWKだと function gxp(x, is, i, c, p, y, z, pi2) { pi2 = 0.398942280401432677940 is = -1 y = abs(x) c = y*y p = 0.0 z = exp(-c*0.5)*pi2 if (y < 2.5) { for (i = 20; i > 0; i--) { p = i*c/(i*2+1+is*p) is = -is } p = 0.5-z*y/(1.0-p) } else { for (i = 20; i > 0; i--) { p = i/(y+p) } p = z/(y+p) } return (x < 0.0) ? 1.0-p : p } 近似式を使う方が計算量は少ないが。

　　　　　[このページのトップへ]

6529.　Re^9: enormpとは何をする関数ですか？　プログラマーA　　2005/05/02 (月) 01:15 あまり，Rには詳しくないので，私が全く知らないような新たな情報が得られ感激しております。非礼はお許しください。特に違う関数を示されたこと，そして，AWKのアルゴリズムで違う方法でお示しになったことは，たいへん参考になっております。御有意義な連休を。

　　　　　[このページのトップへ]

6515.　Re: enormpとは何をする関数ですか？　プログラマーA　　2005/04/30 (土) 20:52 失礼しました。このプログラムは，当然Fortranです。 fBasicというのは，正確には，fBasicsの誤りです。Rの いくつかのライブラリーパッケージのコンパイル前のソース です。

　　　　　[このページのトップへ]

● 「統計学関連なんでもあり」の過去ログ--- 033 の目次へジャンプ
● 「統計学関連なんでもあり」の目次へジャンプ
● 直前のページへ戻る