2022年10月1日土曜日

ズンドコを判定する関数を作った

 先日アップした記事でズンドコキヨシをデータステップで作成した。その時は関数化まではしていなかったが、やはり元ネタに則り関数化までせねばならないと使命感に駆られてしまった。このクソ忙しい時に何をやっているのだろうか。実は私は暇だったのかもしれない。本当か?

お題が「関数化」なのであえてマクロではなくfcmpで関数を定義している。fcmp自体触ったことが無かったので随分と苦戦した。特にfcmpでのarrayはデータステップでのarrayとそのまま同じではないことを初めて知った。まだ二つの違いを理解できていないので難儀なもんである。sasも奥が深い。その違いによるものかは分からないが、乱数で0をズン、1をドコに割り当てたかったが既に数字として定義されているのでarrayの中で文字変数を作れなかったため、やむなく0/1のまま処理をしている。00001がズンズンズンズンドコにあたるので、その場合にキヨシが発生したとみなしている。

ややうまくいかないところがあったが、何はともあれfcmpを使った関数化には成功したので思う存分ズンドコキヨシを判定することができる。多分動くからリリースしようぜとは昔の偉い人もよく言ったものである。関数の出力はズンドコした回数だけにしているのでlogに過度に出力されることもなく、プログラムにこっそり混ぜていてもログが大惨事になることは無い。さあletsズンドコ。この記事が弊社にバレないことを祈る。

/*warning回避用*/

options cmplib = work.piyo ;


proc fcmp outlib = work.hoge.func ;

    function kiyoshi() $40 ;


        __n = 0 ;

        do until (cats(of temp1 - temp5) = "00001") ;


            __n + 1 ;

            /*任意の範囲の乱数を発生させるにはint((上限-下限+1)*乱数+下限)*/

            /*今回は0か1の2値が欲しいので上限=1,下限=0*/

            /*01の乱数をズンとドコに割り当てたかったがうまくいかないので01のままにしている*/

            array ch[5] temp1 - temp5 ;

            do i = 1 to dim(ch) ;

                ch[i] = int( 2*rand("uniform") ) ;

            end ;   


        end ;


        length V_RET $40 ;

        V_RET = cats(__n , "回ズンドコしました") ;


        return(V_RET) ;


    endsub ;

run ;



options cmplib = work.hoge ;

data _null_ ;

    HOGE = kiyoshi() ;

    putlog HOGE ;

run ;


/*実行例*/


880  options cmplib = work.hoge ;

881  data _null_ ;

882      HOGE = kiyoshi() ;

883      putlog HOGE ;

884  run ;


19回ズンドコしました

NOTE: DATAステートメント処理(合計処理時間):

      処理時間           0.04 秒

      CPU時間            0.04 秒


2022年9月1日木曜日

sasで幾何平均を出すときの注意事項

幾何平均は掛け算して根を計算するので値に0や負の数が含まれていると計算できないです.sasのプロシジャで求めるときに値に0が含まれていると,プロシジャによって結果が異なることを紹介します.

幾何平均を出すプロシジャのunivariate,surveymeans,ttestでそれぞれ幾何平均を求めます.使うのは0を含んだ3obsのデータです.私としては0を含むと幾何平均は0になるので,そもそも計算しないことが正しいと思いますが,sasだと0で出すか計算しないか0を除いて計算するかの3つなので,どのように数字を扱いたいかを事前に確認するのがよいと思います.
確認してないですが負の値が含まれていると欠損で値が返ってくるのですかね

/*使うデータ*/
data hoge ;
    a = 5 ; output ;    a = 0 ; output ;  a = 4 ; output ;
run;

/*  univariは0を含む幾何平均を0で出力する  */
proc univariate data=hoge noprint;
 var a;
 output out=out1 geomean=geomean;
run;

NOTE: データセットWORK.OUT1は1オブザベーション、1変数です。
NOTE: PROCEDURE UNIVARIATE処理(合計処理時間):
処理時間 0.00 秒
CPU時間 0.01 秒

data _null_ ;
    set out1 ;
    putlog geomean= ;
run ;

geomean=0
NOTE: データセットWORK.OUT1から1オブザベーションを読み込みました。
NOTE: DATAステートメント処理(合計処理時間):
処理時間 0.00 秒
CPU時間 0.00 秒

/*  surveymeansは0を含むと幾何平均を計算しない  */
proc surveymeans data=hoge geomean ;
  var a ;
  ods output GeometricMeans=out2;
run;

ERROR: A variable must be positive when geometric mean is requested.
NOTE: エラーが発生したため、このステップの処理を中止しました。
NOTE: PROCEDURE SURVEYMEANS処理(合計処理時間):
処理時間 0.00 秒
CPU時間 0.00 秒
WARNING: 出力'GeometricMeans'は作成されていません。出力オブジェクト名、ラベル、
パスが正しく記述されているかを確認してください。また、
要求した出力オブジェクトを作成するために、適切なプロシジャオプションが使われ
ているかも確認してください。たとえば、
NOPRINTオプションが使われていないことを確認してください。

/* ttestは0を含む幾何平均は,0を除いた値で計算してwarningを出す */
proc ttest data=hoge dist=lognormal;
  var a;
  ods output Statistics=out3;
run;

WARNING: 1 observations with invalid response values have been deleted. The response was less than or equal to zero for the
lognormal distribution.
NOTE: データセットWORK.OUT3は1オブザベーション、12変数です。
NOTE: PROCEDURE TTEST処理(合計処理時間):
処理時間 0.31 秒
CPU時間 0.17 秒
data _null_ ;
    set out3 ;
    putlog geommean= ;
run ;

GeomMean=4.4721
NOTE: データセットWORK.OUT3から1オブザベーションを読み込みました。
NOTE: DATAステートメント処理(合計処理時間):
処理時間 0.00 秒
CPU時間 0.00 秒

長々と書きましたが要はログを読んでくれ

2022年8月4日木曜日

データステップでズンドコキヨシ

 このクソ忙しいときに何をやっているんだろうか.いやでも煽られたならなぁ!しょうがねえよなぁ!?Twitterで「どうぞ」なんて言われた日にはやってみざるを得ない.しょうがないんよね.元のお題はこちら.ズンドコキヨシを自作関数で表現したとの話.2016年と昔の話ではあるが…一度見かけてまあだれかsasでもやるやろとスルーしていたのは内緒だ.

元のお題は「関数を作成」ですが今回は普通のデータステップでの記載としました.マクロ化や関数化はあきらめました.FCMPよく分からんし.以下にプログラムと実行結果を示しますが,このプログラムにオサレポイントは特にないです.5個発生させた0か1の2値をとる乱数にそれぞれ0なら「ズン」1なら「ドコ」を割り当て,5変数を結合した文字列が既定の並びになったらループを終了させているだけです.arrayで生成する変数は指定がなければ数値変数になることにお気を付けください.数値変数にはズンとドコの文字列の割り当てができないので終了条件が達成されずに無限ループになります(1敗

data _null_ ;
  /*初期値を与えてもいい*/
  /*call streaminit(123) ;*/

 do until (V_RET="キヨシ!!") ;

    /*任意の範囲の乱数を発生させるにはint((上限-下限+1)*乱数+下限)*/
    /*今回は0か1の2値が欲しいので上限=1,下限=0*/
    length HOGE1 - HOGE5 $4 ;
    array ch HOGE1 - HOGE5 ;
      do i = 1 to dim(ch) ;
        ch(i) = ifc( int( 2*rand("uniform") ) = 0 , "ズン" , "ドコ") ;
      end ;

    if cats(HOGE1 , HOGE2 , HOGE3 , HOGE4 , HOGE5)="ズンズンズンズンドコ" then V_RET="キヨシ!!" ;

    /*logに出力*/
    putlog HOGE1 HOGE2 HOGE3 HOGE4 HOGE5 V_RET ;

 end ;
run ;

/*実行結果 特に上の初期値を使ってないのでこの結果は再現できない*/
/*sas studioじゃ日本語が化けるのでローマ字でやってください*/

ズン ズン ドコ ドコ ドコ ドコ ズン ズン ズン ズン ドコ ズン ズン ドコ ドコ ズン ズン ドコ ズン ズン ドコ ドコ ズン ドコ ズン ズン ズン ドコ ズン ズン ズン ズン ドコ ズン ズン ズン ズン ズン ズン ズン ドコ ズン ズン ズン ズン ドコ ドコ ズン ズン ズン ドコ ドコ ドコ ズン ズン ズン ズン ドコ ドコ ズン ドコ ズン ドコ ドコ ドコ ドコ ドコ ズン ドコ ズン
ズン ズン ズン ズン ドコ キヨシ!! NOTE: DATAステートメント処理(合計処理時間):
処理時間 0.00 秒
CPU時間 0.00 秒


2022年4月27日水曜日

変更後のエスケープ文字を確認する話

 sasはエスケープ文字をデフォルトから変えることができます.例えばODS ESCAPECHAR="^" ;などと指定すると,[^]になりますね.自分で直近に作ったプログラムならエスケープ文字を何にしたかは大体見当がつきますが,他の人のを引き継いだり,随分前に作ったプログラムだとエスケープ文字が何かわからないことがあります.だいたい使いそうな文字をエスケープ文字に指定することは無い(はず)ですが、適当に変えるのも怖いので事前に確認はしたいところです

確認するには頑張ってエスケープ文字を指定している箇所を見に行くか,現在のエスケープ文字を格納している自動マクロ変数を見に行くかのどっちかですね.この自動マクロ変数は最近知りました.まだ便利だという状況に出会ったことが無いですがエスケープ文字の指定している箇所を探しに行ったことは何度かあるので,折があれば便利に使えるかなと楽しみです.

4     ODS ESCAPECHAR="^" ;

5    %put &=SYSODSESCAPECHAR ;

SYSODSESCAPECHAR=^