UNIVARIATEプロシジャ

例4.25 折り重ねられた正規曲線の追加表示

この例は、HISTOGRAMステートメントでサポートされていない当てはめた曲線の表示方法を示しています。多くの製造済み組み立て品の取り付け位置のオフセット(mm)が測定され、測定値(Offset)はAssemblyという名前のデータセットに保存されます。次のステートメントはデータセットAssemblyを作成します。

data Assembly;
   label Offset = 'Offset (in mm)';
   input Offset @@;
   datalines;
11.11 13.07 11.42  3.92 11.08  5.40 11.22 14.69  6.27  9.76
 9.18  5.07  3.51 16.65 14.10  9.69 16.61  5.67  2.89  8.13
 9.97  3.28 13.03 13.78  3.13  9.53  4.58  7.94 13.51 11.43
11.98  3.90  7.67  4.32 12.69  6.17 11.48  2.82 20.42  1.01
 3.18  6.02  6.63  1.72  2.42 11.32 16.49  1.22  9.13  3.34
 1.29  1.70  0.65  2.62  2.04 11.08 18.85 11.94  8.34  2.07
 0.31  8.91 13.62 14.94  4.83 16.84  7.09  3.37  0.49 15.19
 5.16  4.14  1.92 12.70  1.97  2.10  9.38  3.18  4.18  7.22
15.84 10.85  2.35  1.93  9.19  1.39 11.40 12.20 16.07  9.23
 0.05  2.15  1.95  4.39  0.48 10.16  4.81  8.28  5.68 22.81
 0.23  0.38 12.71  0.06 10.11 18.38  5.53  9.36  9.32  3.63
12.93 10.39  2.05 15.49  8.12  9.52  7.77 10.70  6.37  1.91
 8.60 22.22  1.74  5.84 12.90 13.06  5.08  2.09  6.41  1.40
15.60  2.36  3.97  6.17  0.62  8.56  9.36 10.19  7.16  2.37
12.91  0.95  0.89  3.82  7.86  5.33 12.92  2.64  7.92 14.06
;

折り重ねられた正規分布をオフセットの測定に当てはめることが決定されます。変数 $\text{[math]}$ は $\text{[math]}$ の場合に折り重ねられた正規分布になります。ここで、 $\text{[math]}$ は $\text{[math]}$ のような分布です。当てはめた密度は次のとおりです。

$\text{[math]}$

ここで、 $\text{[math]}$ となります。

SAS/IMLを使用して、Elandt (1961)が定義した積率法に基づく $\text{[math]}$ および $\text{[math]}$ の予備的な推定値を計算できます。これらの推定値はElandt (1961)の(19)を解くことによって計算されます。この式は次のように定義されます。

$\text{[math]}$

ここで、 $\text{[math]}$ は標準正規分布関数です。

$\text{[math]}$

次に、 $\text{[math]}$ および $\text{[math]}$ の推定値は次のように求められます。

$\text{[math]}$

まず、PROC MEANSによる1番目および2番目の積率の計算と、次のDATAステップによる定数 $\text{[math]}$ の計算を行います。

proc means data = Assembly noprint;
   var Offset;
   output out=stat mean=m1 var=var n=n min = min;
run;

* Compute constant A from equation (19) of Elandt (1961);
data stat;
   keep m2 a min;
   set stat;
   a  = (m1*m1);
   m2 = ((n-1)/n)*var + a;
   a  = a/m2;
run;

次に、SAS/IMLサブルーチンNLPDDを使用して、 $\text{[math]}$ を最小化することにより式(19)を解き、 $\text{[math]}$ および $\text{[math]}$ を計算します。

proc iml;
   use stat;
   read all var {m2}  into m2;
   read all var {a}   into a;
   read all var {min} into min;

   * f(t) is the function in equation (19) of Elandt (1961);
   start f(t) global(a);
     y = .39894*exp(-0.5*t*t);
     y = (2*y-(t*(1-2*probnorm(t))))**2/(1+t*t);
     y = (y-a)**2;
     return(y);
   finish;

   * Minimize (f(t)-A)**2 and estimate mu and sigma;
   if ( min < 0 ) then do;
      print "Warning: Observations are not all nonnegative.";
      print "     The folded normal is inappropriate.";
      stop;
      end;
   if ( a < 0.637 ) then do;
      print "Warning: the folded normal may be inappropriate";
      end;
   opt = { 0 0 };
   con = { 1e-6 };
   x0  = { 2.0 };
   tc  = { . . . . . 1e-8 . . . . . . .};
   call nlpdd(rc,etheta0,"f",x0,opt,con,tc);
   esig0 = sqrt(m2/(1+etheta0*etheta0));
   emu0  = etheta0*esig0;

   create prelim var {emu0 esig0 etheta0};
   append;
   close prelim;

   * Define the log likelihood of the folded normal;
   start g(p) global(x);
      y = 0.0;
      do i = 1 to nrow(x);
         z = exp( (-0.5/p[2])*(x[i]-p[1])*(x[i]-p[1]) );
         z = z + exp( (-0.5/p[2])*(x[i]+p[1])*(x[i]+p[1]) );
         y = y + log(z);
      end;
      y = y - nrow(x)*log( sqrt( p[2] ) );
      return(y);
   finish;
   * Maximize the log likelihood with subroutine NLPDD;
   use assembly;
   read all var {offset} into x;
   esig0sq = esig0*esig0;
   x0      = emu0||esig0sq;
   opt     = { 1 0 };
   con     = { . 0.0, .  .  };
   call nlpdd(rc,xr,"g",x0,opt,con);
   emu     = xr[1];
   esig    = sqrt(xr[2]);
   etheta  = emu/esig;
   create parmest var{emu esig etheta};
   append;
   close parmest;
quit;

出力4.25.1に示されているように、予備的な推定値がデータセットPrelimに保存されます。

出力4.25.1 $\text{[math]}$ 、 $\text{[math]}$ および $\text{[math]}$ の予備的な推定値

Obs	EMU0	ESIG0	ETHETA0
1	6.51735	6.54953	0.99509

ここで、 $\text{[math]}$ および $\text{[math]}$ を初期推定値として使用して、NLPDDサブルーチンを呼び出し、折り重ねられた正規分布の対数尤度 $\text{[math]}$ を(ここでは定数になるまで)最大化します。

$\text{[math]}$

* Define the log likelihood of the folded normal;
start g(p) global(x);
   y = 0.0;
   do i = 1 to nrow(x);
      z = exp( (-0.5/p[2])*(x[i]-p[1])*(x[i]-p[1]) );
      z = z + exp( (-0.5/p[2])*(x[i]+p[1])*(x[i]+p[1]) );
      y = y + log(z);
   end;
   y = y - nrow(x)*log( sqrt( p[2] ) );
   return(y);
finish;
* Maximize the log likelihood with subroutine NLPDD;
use assembly;
read all var {offset} into x;
esig0sq = esig0*esig0;
x0      = emu0||esig0sq;
opt     = { 1 0 };
con     = { . 0.0, .  .  };
call nlpdd(rc,xr,"g",x0,opt,con);
emu     = xr[1];
esig    = sqrt(xr[2]);
etheta  = emu/esig;
create parmest var{emu esig etheta};
append;
close parmest;
quit;

出力4.25.2に示されているように、データセットParmEstに最尤推定値 $\text{[math]}$ 、 $\text{[math]}$ (および $\text{[math]}$ )が格納されます。

出力4.25.2 $\text{[math]}$ 、 $\text{[math]}$ および $\text{[math]}$ の最終推定値

Obs	EMU	ESIG	ETHETA
1	6.66761	6.39650	1.04239

曲線をヒストグラムに追加するため、まずヒストグラム間隔の幅と終点を計算します。次のステートメントは、これらの値をOutCalcというデータセットに保存します。なお、プロットはこの時点では作成されません。

proc univariate data = Assembly noprint;
   histogram Offset / outhistogram = out normal(noprint) noplot;
run;

data OutCalc (drop = _MIDPT_);
   set out (keep = _MIDPT_) end = eof;
   retain _MIDPT1_ _WIDTH_;
   if _N_ = 1 then _MIDPT1_ = _MIDPT_;
   if eof then do;
      _MIDPTN_ = _MIDPT_;
      _WIDTH_ = (_MIDPTN_ - _MIDPT1_) / (_N_ - 1);
      output;
   end;
run;

出力4.25.3はデータセットOutCalcのリストを示しています。ヒストグラムの棒の幅は、変数_WIDTH_の値として保存されます。最初と最後のヒストグラムの棒の中間点は、変数_MIDPT1_および_MIDPTN_の値として保存されます。

出力4.25.3 データセットOutCalc

Obs	_MIDPT1_	_WIDTH_	_MIDPTN_
1	1.5	3	22.5

次のステートメントはAnnoという名前の注釈データセットを作成し、そのデータセットに当てはめた曲線の座標を格納します。

data Anno;
   merge ParmEst OutCalc;
   length function color $ 8;
   function = 'point';
   color    = 'black';
   size     =  2;
   xsys     = '2';
   ysys     = '2';
   when     = 'a';
   constant = 39.894*_width_;;
   left     =  _midpt1_ - .5*_width_;
   right    =  _midptn_ + .5*_width_;
   inc      = (right-left)/100;
   do x = left to right by inc;
      z1 = (x-emu)/esig;
      z2 = (x+emu)/esig;
      y  = (constant/esig)*(exp(-0.5*z1*z1)+exp(-0.5*z2*z2));
      output;
      function = 'draw';
   end;
run;

次のステートメントはANNOTATE=データセットを読み込み、ヒストグラムと当てはめた曲線を表示します。

title 'Folded Normal Distribution';
ods graphics off;
proc univariate data=assembly noprint;
   histogram Offset / annotate = anno;
run;

出力4.25.4はヒストグラムと当てはめた曲線を表示しています。

出力4.25.4 折り重ねられた正規曲線を追加表示したヒストグラム

この例のサンプルプログラムuniex15.sasは、Base SASソフトウェアのSASサンプルライブラリに含まれています。