FCMP.SOLVEで方程式の実数解をガウス・ニュートン法で求める

 まずSASの公式に乗ってる例から

y = 1 / √x

でy=20となるxの値を求めてみます

グラフにすると

なので，まあ　だいたい値は0.0023～26あたりかなぁって感じですが

proc fcmp outlib=work.common.solve ;

/*解きたい式を定義*/

   function f1(x);

     y = 1/sqrt(x);

      return(y);

   endsub;

 

/*yを指定して、xを求めて戻すルーチン*/

  subroutine solver(y,x);

   outargs x  ;

   x=solve("f1",{.}, y , .);

  endsub ;

run;

options cmplib=work.common;

data wk1;

 call missing(x);

 call solver(20,x);

run;

proc print data=wk1;

run;

👆
と，このように，あらかじめ解きたい式をFCMPで指定し，
FCMP.SOLVE関数の第一引数に　解きたい関数名，第三引数に目標値（ここでいうところの左辺，y）を指定します，
第二引数には計算のための5つのパラメータを配列で与えるのですが，{.}とするとデフォルト設定になります
最後の引数は，第一引数にとった関数に追加引数をここで与えることができます

つぎに3次方程式

を解いてみます

グラフにすると

👆どうも実数解は3つありそうです．ぱっとみ，交点は1,2,3らへん

この場合，初期値の設定によって，収束する解が異なるので，コールルーチンをちょっと改造して，パラメータを細かく指定できるようにして実行してみます

proc fcmp outlib=work.common.solve ;

/*解きたい式を定義*/

   function f1(x);

      return(x**3 -6*x**2 + 11*x -6);

   endsub;

 

/*yを指定して、xを求めて戻すルーチン*/

  subroutine solver(y,x,initial_value,absolute_criterion,relative_criterion,maximum_iterations,solve_status);

  /*設定パラメータ */

   array opts[5] initial_value absolute_criterion relative_criterion maximum_iterations solve_status;

   outargs x ,solve_status,y ;

   x=solve("f1"

               ,opts

               , 0

               , .);

   put x= solve_status=;

  endsub ;

run;

options cmplib=work.common;

data wk1;

 call missing(of x);

 absolute_criterion=1.0e-12;

 relative_criterion=1.0e-6;

 maximum_iterations=100;

 solve_status=.;

 y=0;

 initial_value=0.001;

 call solver(20,x,initial_value,absolute_criterion,relative_criterion,maximum_iterations,solve_status);

 output;

 initial_value=1.9;

 call solver(20,x,initial_value,absolute_criterion,relative_criterion,maximum_iterations,solve_status);

 output;

 initial_value=2.9;

 call solver(20,x,initial_value,absolute_criterion,relative_criterion,maximum_iterations,solve_status);

 output;

run;

proc print data=wk1;

run;

グラフから，おそらくｘ＝1,2,3あたりに解がありそうなので
初期値(initial_value)を，0.001，1.9, 2.9として回してみると

👆
やはり ｘ=1,2,3が解でした．solve_status=0は収束していることを示します

UpSet plot (アップセット プロット)の実装例

 UpSet-Plotというプロットをよくみかけます．アイテムの組み合わせ集計を，1個での集計，2個のアイテムセットでの集計，3個でのアイテムセットの集計...... n個のアイテムセットの集計　って，バスケット分析のような集計をしたときにそれをグラフ化します．
星とり表の部分がアイテムセットを示して，上のヒストグラムはそのカウント数です．
ここで，たとえば以下のグラフはある仮想の疾患に罹患した際に，患者に発現した症状のプロットです．Aさんが発熱と倦怠感の症状があったい場合，「発熱＆倦怠感」のアイテムセットとして集計されます

👆
ここで注意は
上図のヒストグラム(上側)の一番左は「咳嗽」で，●が一つですが，これは，「咳嗽」だけを症状として，もった人がカウントされます．「咳嗽」と「発熱」がある人は「咳嗽」にはカウントされずに「咳嗽＆発熱」でカウントされます．

では全体として，「咳嗽」を少なくとも1症状として持つ人数は？？となった場合は横に延びてるヒストグラムがそれにあたるという仕組みです

SASコードは以下です（冒頭のグラフ描いたやつです）．一応，データが変わっても，ID , item, itemnameをあわせてもらえれば可変で動くとは思いますが，多少調整いるかもです．

もうちょっと，洗練されてどこかで出そうと思いますが，第一ドラフトということで．
barchartparamで上側のヒストグラム，scatterとhighlowで星取表部分，横のヒストグラムはhighlowで描いてます．

/*テストデータ*/

data test;

call streaminit(123);

do ID = 1 to 100;

  iend=rand("integer",1,3);

  do i = 1 to iend;

    item = rand("integer",1,5);

    itemname = choosec(item,"鼻汁","倦怠感","咽頭痛","咳嗽","発熱");

    output;

  end;

end;

drop i iend;

run;

/*アイテム数*/

proc sql noprint;

 select count(distinct item) into: item_n trimmed from test;

quit;

%put &=item_n;

/*アイテム*/

proc sort data=test(keep=item itemname) out=itemlist nodupkey;

  by  item;

run;

/*アイテムformat*/

data itemformat;

 set itemlist;

 FMTNAME="itemformat";

 START=item;

 LABEL=itemname;

run;

proc format lib=work cntlin=itemformat;

run;

/*個人内での重複除去*/

proc sort data=test out=test1 nodupkey;

  by ID item;

run;

/*アイテム単純集計*/

proc summary data=test1 nway ;

 class item itemname;

 output out=_single_count(drop=_TYPE_ rename=(_FREQ_=single_count) );

run;

data single_count;

 set _single_count;

 single_start=0;

 run;

/*個人内での組み合わせ生成のため転置*/

proc transpose data=test1 out=out1 prefix=item_;

 by ID;

 id item;

 var item;

run;

/*組み合わせ生成*/

data out2;

length comb $200.;

 set out1;

 array ar item_:;

 do over ar;

  if ^missing(ar) then comb=catx("/",comb,ar);

 end;

  count=1;

run;

/*組み合わせ集計*/

proc summary data=out2 nway missing;

 class comb: item_:;

 output out=_comb_count(drop=_TYPE_ rename=(_FREQ_=comb_count) );

run;

 proc sort data=_comb_count;

 by descending comb_count;

run;

/*組み合わせ数と最大カウント*/

proc sql noprint;

 select count(*) into:comb_n from comb_count;

 select max(comb_count) into:comb_max from comb_count;

quit;

/*星取り表のつなぎ線用に最小・最大値*/

data comb_count;

set _comb_count;

 array dummy_{&item_n};

 do i = 1 to &item_n;

    dummy_{i}=i;

 end;

min = min(of item_:);

max = max(of item_:);

x = _N_;

drop i;

run;

/*グラフ用データ完成*/

data graph_data;

 set comb_count single_count;

run;

/*星取表*/

%macro scatterplot();

 %do i = 1 %to &item_n;

    scatterplot x=x y=item_&i /yaxis=y2 markerattrs=(size=12 symbol=circlefilled);

    scatterplot x=x y=dummy_&i /yaxis=y2 datatransparency=0.9 markerattrs=(size=12 symbol=circlefilled);

 %end;

%mend;

proc template ;

  define statgraph upset;

      begingraph ;

       /*2*2レイアウト*/

          layout lattice/ rowdatarange = union

                             columnweights = (0.25 0.75)

                             rowweights = (0.6 0.4)

                             columns = 2 rows = 2 

                             pad=(top=0 bottom=0 left=0 right=0) 

                             ; 

              /*余白 */

              layout overlay / pad=(top=0 bottom=0 left=0 right=0)

                                     outerpad=(top=0 bottom=0 left=0 right=0) 

                                     xaxisopts  = (display =none);

              endlayout;

             

              /*ヒストグラム（barchartparm）*/

              layout overlay /pad=(top=0 bottom=0 left=0 right=0) 

                                   outerpad=(top=0 bottom=0 left=0 right=0)

                                    walldisplay=none

                                    xaxisopts  = (display =none linearopts = (viewmin=1 viewmax=&item_n tickvaluesequence = (start=1 end=&item_n increment=1)))

                                    yaxisopts  = (display= (line tickvalues) linearopts = (tickvaluesequence = (start=0 end=&comb_max. increment=5))) ;

                                   barchartparm category=x response=comb_count/datalabel=comb_count;

              endlayout;

              /*横ヒストグラム(highlowplot)*/

              layout overlay /pad=(top=0 bottom=0 left=0 right=0)

                                    outerpad=(top=0 bottom=0 left=0 right=0)

                                    walldisplay=none

                                    xaxisopts  = (display= none reverse=true)

                                    yaxisopts = (display =none linearopts = (viewmin=1 viewmax=&item_n tickvaluesequence = (start=1 end=&item_n increment=1)))

                                    y2axisopts = (display =(tickvalues) tickvalueattrs=(size=9) linearopts = (viewmin=1 viewmax=&item_n tickvaluesequence = (start=1 end=&item_n increment=1)));

                                    highlowplot y=item high=single_count low=single_start 

                                    /yaxis=y2 group=item type=bar barwidth=0.5 fillattrs=(color=VLIGB) outlineattrs=(color=VLIGB) 

                                     highlabel=single_count labelattrs=(color=black)

                                    ;

              endlayout;

              /*星取表(scatterplot とhighlowplot)*/

              layout overlay  /pad=(top=0 bottom=0 left=0 right=0)

                                     outerpad=(top=0 bottom=0 left=0 right=0)

                                     xaxisopts  = (display =none linearopts = (tickvaluesequence = (start=1 end=&comb_n. increment=1)))

                                     yaxisopts  = (display= (line ) linearopts = (viewmin=1 viewmax=&item_n  ))

                                     y2axisopts  = (display=none linearopts = (viewmin=1 viewmax=&item_n  ))

                                     walldisplay=none;

                                     %scatterplot();

                                     highlowplot x=x low=min high=max/yaxis=y2;

              endlayout;

           endlayout;

      endgraph;

  end;

run;

* グラフ作成実行 ;

proc sgrender data=graph_data template=upset;

format item itemformat.;

run;

IMLを使わずFCMPで行列計算⑦ コレスキー分解

proc fcmp;

   array mat1[3,3] 2 2 3 2 4 2 3 2 6;

   array mat2[3,3] / nosym;

   call chol(mat1, mat2, 0);

   rc = write_array("MAT1", mat1 ) ;

   rc = write_array("MAT2", mat2 ) ;

run;

proc print data=MAT1;

run;

proc print data=MAT2;

run;

mat2の転置行列を作って，mat2とmat2tの行列積をつくればmat1に戻ることが確認できる
つまり，コレスキー分解は　下三角行列とその転置の積に分解する操作

 

IMLを使わずFCMPで行列計算⑥ 行列式の算出

 proc fcmp;

   array mat1[2,2] (2 1 3 4);

   call det(mat1, result);

   put result=;

quit;

result = 5

2次正方行列なので

２ * 4  - 1 * 3 = 5

IMLを使わずFCMPで行列計算⑤ 単位行列の生成

proc fcmp;

  array mat1[3,3] /nosym;

   call identity(mat1);

   rc = write_array("RESULT", mat1 ) ;

quit;

 

IMLを使わずFCMPで行列計算④ 転置行列

 proc fcmp;

   array mat1[3,2] (3,5,0,9,1,2);

   array trans_mat1[2,3] /nosym;

  /*call transposeでmat1の転置行列をtrans_mat1に格納 */

   call transpose(mat1,trans_mat1);

   rc = write_array("MAT1", mat1 ) ;

   rc = write_array("TRANS_MAT1", trans_mat1 ) ;

quit;

proc print data=MAT1;

run;

proc print data=TRANS_MAT1;

run;

IMLを使わずFCMPで行列計算③ 逆行列

 proc fcmp;

   array mat1[3,3] (3,5,0,9,1,2,8,7,9);

   array inv_mat1[3,3] /nosym;

   array result[3,3] /nosym;

  /*call invでmat1の逆行列をinv_mat1に格納 */

   call inv(mat1,inv_mat1);

  /*mat1と逆行列の行列積をresultに格納*/

   call mult(mat1, inv_mat1, result);

   rc = write_array("MAT1", mat1 ) ;

   rc = write_array("INV_MAT1", inv_mat1 ) ;

   rc = write_array("RESULT", result ) ;

quit;

proc print data=MAT1;

run;

proc print data=INV_MAT1;

run;

proc print data=RESULT;

run;

浮動小数点誤差で-0とかってなってますが，逆行列かけたら　単位行列になるのが確認できる