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=^

2021年7月16日金曜日

warterfall plotの棒が少ない時に間隔を狭くする話

 オプションの一発で解決する内容なのですが,機会があったので話のネタにしようと.…身バレしそうやけどこんな場末ブログ誰も見とらんやろ

以下のような3本しかbarが出てないので間隔が広いのがどうにかならんかと言われたとします.


こういう時は私はxaxisステートメントのoffsetmin/maxで調整しています.ほかにいい方法があるかもしれませんが.今回だとoffsetmin = 0.3    offsetmax = 0.3 くらいを指定すると以下の図になり多少は棒の間隔がマシになります.今回は0.4まで広げると個人的に見た目がイマイチだったので0.3か0.35くらいですかね.正直間隔が広くてもええやろって思うのですがなかなか許してもらえない


この状態で棒の幅をvbarステートメントでbarwidth = 1.0を指定すると,以下の感じに棒の
間隔そのものがなくなります.offsetmin/maxである程度見た目整えて,微調整をbarwidthでやるようにしています.初めて棒くっつけろ言われた時はどうしようか悩んだのが懐かしいですね.頼むから余計なリクエスト言わんでくれ…


2021年6月10日木曜日

特定文字のunicodeが知りたい話

 βとかαなどギリシア文字に半角用フォントを充てるときや,≧や③といった文字を出力するときにunicodeを直接指定するときがあると思います.βなら^{unicode 2265}だったりする4文字のやつですね.毎度unicodeを検索しなくてもunicodec関数を使用すればsasである程度は表示できるので,その関数の紹介です.

例えば以下の実行結果ですが,上で取り上げたβと≧のunicodeをlogに出しています.

131    data _null_ ;
132      str1 = unicodec("β" ) ;
133      str2 = unicodec("≧" ) ;
134      putlog "β:" str1 "≧:" str2 ;
135    run ;

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

β:\u03B2 ≧:\u2267

他にも③だと以下の通りです.

136    data _null_ ;
137      str1 = unicodec("③" ) ;
138      putlog "③:" str1 ;
139    run ;

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

③:\u2462

分からなくて毎回調べるのが大変だったのでユニコードを出力できるこのunicodec関数は結構助かっています.sasのhelpを見ると「This function reads characters that are in the current SAS session encoding and converts them to Unicode encoding.」とあるので,現在のsasのセッションエンコーディングの文字を読み取って,unicodeに変換する関数のようです.なのでsjisの環境だと上付きの数字などのunicodeを出せないのが残念でならないです…上付きの数字は結構使うので惜しいところです.unicode環境で上付き文字を指定するとunicodeは出せると思います.試してないですが.

ちなみにunocodec関数には第二引数があり出力時のエスケープ文字を変えれるのですが,今回のようなユニコードを調べるだけならデフォルトのままで不都合がないのでそのままにしています.例えば第二引数に"paren"と指定すると<03B2>のようにカッコくくりになります

2021年6月9日水曜日

sas発表会に参加した

 前回延期したのが2020年1月末,実にそこから一年強空いたがsasの発表会に参加した.昨今の情勢もありオンライン開催となったが,会場を提供してくださった方に多大な感謝を.

今回の私の発表はswimmerplotについて…こいついつもグラフ描いてんな,なんてことは言ってはいけない.最近は闇の勢力によってexcelでKM曲線出せだの,excelでforestplot出せだの摩訶不思議なworldに足を踏み入れている.地獄か?最近はやりのspotfire君は早く私を助けてくれ.excelでKM曲線書くときのヒゲをどうするかと言う偉大な問題は人力で解決したが,設定で出せるやろVS設定考えてる暇あったら人力でやってしまえ,どうせこんな話二度とこない.の葛藤の末のことだった.何とか良い感じに処理できるように工夫したいところ.

今回の発表会も奇々怪々,いや失敬,多様な内容だった.複数のグラフを1枚のwordに並べるときによくわからない事象が起きるのはどこも共通のようだった.あと野球のシミュレーションの話がなかなか面白そう.ああいうの見るとちょっとやりたくなりますなります.

次回開催は未定.ユーザー会が秋なので冬にするか,いっそ夏にするか,夏だと発表の準備が間に合わんから冬だったらうれしいなって思う.

2021年のsasユーザー会のアブストの締め切りが7月にあるが,正直気づいてなかったし今年の私の予定がやべーので多分間に合わない.あきらめの境地である.オンライン開催なら有給で潜り込んでもバレへんかな…などと考えているが,やっぱり今年も申請して参加するのだろう.予定が合えばだが.

最後に共有フォルダの設定がポンコツなのは修正しますので少しお待ちください…すみません…