データステップ中に一時的に使って、自動的にドロップしたい変数の話

データステップ中に一時的に作成するけど、完成したデータセットからは落としたい変数の定義法について困っている方がいたので考えてみました。

ひとつは、一律、変数の接頭にXX_とかつけて、最後にdrop XX_: とかで一括ドロップしてもいいと思います。

もうひとつは、以下のようにテンポラリー配列で定義した配列に値を出し入れしたりする方法とかどうでしょう?

以下のコードは、意味のない処理ですが、色々やっても結局データセットにはxとyだけが残ります。

いかがでしょう?

data A1;
array temp{3} _temporary_;

x=1;

temp{1}=10;
temp{2}=15;
temp{3}=x;

x=temp{1}+1;
temp{3}=temp{3}+x;
y=temp{3}*temp{1}+x;

run;

毎月第2火曜日は卵が100円の日だからマスタ更新に流し込むデータセット作るときに 気をつけてねみたいなこと言われたらの話

昔、SASを覚え始めた頃に何年何月の第何週の何曜日を起点としてナントカみたいに何月何週目の何曜日という条件で日付を使って処理する複雑なコードを書いていて気が狂いそうになったことがありました。
病院にデータ入力の督促だか、なんか回収の管理だかそんな系で。

バッチでまわすので、日付を動的に判断する必要がありました。

どんな処理だったかもうあんまり覚えてないですが今ならnwkdom関数を使っていたと思います。

nwkdom関数は、以下のコードと結果を見てもらえたらすぐわかるように何年何月の第何週の何曜日というものを全てパラメータで与えてやると該当の日付を返してくれるという関数です。

data A1;
 /* 
  第一引数:第何週か
  第二引数:曜日指定
 1=日曜
   2=月曜   
 3=火曜
 4=水曜
 5=木曜
 6=金曜
 7=土曜
第三引数 :何月か
第四引数 :何年か
*/

/*第1週~5週の月曜日の日付*/
  do week= 1 to 5;
  x=nwkdom(week,2,2,2016);output;
  end;

/*第1週~5週の火曜日の日付*/
  do week= 1 to 5;
  x=nwkdom(week,3,2,2016);output;
  end;

  format x nldatew30.;
run;



















臨床試験分野でコード書いていると、ふつう投与日等からのDayXXで考えることが多いので便利さがピンとこないですが、他の分野だと、例えば第何週の何曜日はセールだから売り上げ予測がナントカとかいった処理があったりして結構役立つんですね。

注意点としては、ちょっと癖のある仕様で、例えば、先のコードにも入れましたが2016年2月の第5週の火曜日を指定した場合、2016年2月には第5週月曜日までしかないので、代わりに第4週の火曜日の日付が返ります。

つまり第一引数に5を指定した場合、5週目の指定曜日がなかった場合は4週の値、イコールその月の最後の指定曜日の日付が返るというわけです。

よく月の最後の何曜日にナントカの処理があるとかって場合は、5を指定しておけば必ず指定曜日の最終日付がとれるわけです。

4と5に同じ日付が戻ることにを意識の片隅においておかないと、冷凍食品4割引をさらに4割引したりして投売りみたいなことになりかねないので気をつけてくださいね

強さとはなにか。イロレーティングの話

最近教えて貰ったんですが、ドラゴンボールのフリーザって、最近の映画で復活して、体が金色になって戦闘力が1垓(京の上、100000000000000000000)になったそうです。
僕が知っているのは「私の戦闘力は53万です」ですげーっ!!て時代だったんで、そんなジンバブエドルみたいな単位になってるとは知りませんでした

しかし、スカウターみたいな便利なものがあればいいですが(よく爆発するけど)、強さといった複合的要因でなりたっているものをスコア化するのはなかなか難しいです。
また、普段からスカウター慣れしていないと算出された戦闘力の差が、どの程度ならどのくらい勝てるかいうことが実感できないです。
(正確な値だしているのに、ちっ壊れてやがる!グシャっていうシーンあった気もするけど、慣れすぎているが故の弊害でしょうね)


そういう実際の能力を数値化する方法と違うアプローチとして、測定対象の対戦結果データがある程度存在する場合、強さ(勝率と捉える)が正規分布に従うと仮定し、数値化した値の差から勝率が逆算できるように作られたイロレーティングというものがあります。

原理など詳しくはwikiや参考ページなど検索してみてください。

【wikipedia】
https://ja.wikipedia.org/wiki/%E3%82%A4%E3%83%AD%E3%83%AC%E3%83%BC%E3%83%86%E3%82%A3%E3%83%B3%E3%82%B0

【ヘキサドライブ(ゲーム会社)公式ブログ】
http://hexadrive.sblo.jp/article/67433861.html


イロレーティングの肝となるのは以下三点。

・ゲームの結果は一方の勝ち、一方の負けのみとし、引き分けは考慮しない(0.5勝0.5敗と扱うものとする)。
・200点のレート差がある対局者間では、レートの高い側が約76パーセントの確率で勝利する。
・平均的な対局者のレートを1500とする。

特に2点目、レーティングの差から、仮に相手を見たこともなくても勝率が推定できるという
のがいいんです。

もともとチェスで運用されていたもので、今はネット将棋、囲碁など主たるサイトはこれを基本にそれぞれカスタマイズして使ってユーザーの実力を測り、対戦マッチング等に利用しています。

僕もネット将棋やるのですが、これが結構よくできてるなぁって実感します。
実際、後で自分の対戦成績と対戦相手のレーティングで集計してみると200差の相手には
3割勝ててないですし、レーティングで100あいてると、やっててもちょっと自分より強いなって感覚を持ちます

算出法はいたって簡単、何局か後でまとめて計算するのと都度計算するので多少ロジック代わりますが、

プレイヤー1のレーティングがrate1、プレーヤー2がrate2とすると
プレイヤー1が勝利する確率e1は以下で算出できます
e1 = 1/( 1 + 10 ** ((rate2-rate1)/400) );

そして、対局後のレーティングは,
勝利の場合S=1 敗北はS=0 引き分けはS=0.5として
対局後のレーティング = 現在のレーティング + 16 * ( S  - e1);

16は定数で、プロなら16、アマなら32にすると安定するといわれています

例えば以下のように現在のレーティングデータと
data Rating;
length name $10.; 
input name rate;
cards;
A 1500
B 1500
C 1500
D 1500
;
run;












時系列に並んだ対局結果のデータセットがあった場合
(1はプレイヤー1が勝者、2は2が勝者、3は引き分け)
data hosi;
length player1 player2 $10.;
input player1 player2 win;
winner=choosec(win,player1,player2,'引き分け');
cards;
A B 1
B A 2
D A 2
B C 1
C B 2
B D 1
D B 1
A C 1
A D 1
C D 3
D C 1
C A 2
D C 1
C D 2
B C 1
;
run;




























コードは

data cal;
if _N_=0 then set rating;
if _N_=1 then do;
declare hash h1(dataset:'rating',ordered:'Y');
h1.definekey('name');
h1.definedata('name','rate');
h1.definedone();
end;
set hosi end=eof;

if h1.check(key:player1) ne 0 then do;
h1.add(key:player1,data:player1,data:1500);
end;
if h1.check(key:player2) ne 0 then do;
h1.add(key:player2,data:player2,data:1500);
end;

h1.find(key:player1);
rate1=rate;
h1.find(key:player2);
rate2=rate;

e1 = 1/( 1 + 10 ** ((rate2-rate1)/400) );
e2 = 1/( 1 + 10 ** ((rate1-rate2)/400) );
rate1 = rate1 + 16*(ifn(win=1,1,ifn(win=3,0.5,0)) - e1);
rate2 = rate2 + 16*(ifn(win=2,1,ifn(win=3,0.5,0)) - e2);

h1.replace(key:player1,data:player1,data:rate1);
h1.replace(key:player2,data:player2,data:rate2);

h1.output(dataset:cats('rating',_N_));

run;

とすれば、1対局ごとのレーティングの状態をデータセット化できます
(わかりやすくするために1局ごとに作ってますが、普通はif eofで最後の
時だけratingのデータセットを吐けばいいです)

1局目、A対BでAが勝利した後のデータセットをみてみると、
このようにAが少しあがり、Bが少し下がります。












最終的に全対局後をみると以下のようになり











全勝のAと全敗のCでは100ほど開いています。
100だと勝率64%ぐらいです

ある程度データがたまらないと安定したレーティングは得られませんが
まあ、AとCが戦ったら6割5分Aが勝ちそうだというのは感覚的に納得できます

ちなみに式をみればわかりやすいですが、レーティング差が大きい相手との
対局ほど、移動するレーティングが大きくなります。
雑魚に負けた方はレーティングが大きく下がり、下克上した方は大きく上がります。
わかりやすいですね。

イロレーティングの改良法については色々研究されているみたいなので興味のある方はぜひ
勉強してみてください

2地点間の緯度経度から直線距離でもっとも近いデータを取得する話 geodist関数

親戚の子供の宿題を手伝っていたんですが、以下の問題に頭抱えました。

「県庁所在地で、直線距離がもっとも近いのは何県と何県?
地図帳とものさしを使って考えてみよう」

うわぁ~、凄いめんどい。

「名古屋と岐阜であってるよね」って言われても、そんなのわかんない。
大阪-神戸だと思ったけど自信なし

あたりをつけて、さし当ててみるけど、結構難しいぞこれ。
わかります??

困った時のSAS On Demand
ネットさえつながればどこでも自由にSAS使えるので、ちょっとPC借りて
以下のコードをパチパチ。

data q1;
length no 8. shi $10. ido keido 8.;
input no shi ido keido;
cards;
1 札幌市 43.063968 141.347899
2 青森市 40.824623 140.740593
3 盛岡市 39.703531 141.152667
4 仙台市 38.268839 140.872103
5 秋田市 39.7186 140.102334
6 山形市 38.240437 140.363634
7 福島市 37.750299 140.467521
8 水戸市 36.341813 140.446793
9 宇都宮市 36.565725 139.883565
10 前橋市 36.391208 139.060156
11 さいたま市 35.857428 139.648933
12 千葉市 35.605058 140.123308
13 新宿区 35.689521 139.691704
14 横浜市 35.447753 139.642514
15 新潟市 37.902418 139.023221
16 富山市 36.69529 137.211338
17 金沢市 36.594682 136.625573
18 福井市 36.065219 136.221642
19 甲府市 35.664158 138.568449
20 長野市 36.651289 138.181224
21 岐阜市 35.391227 136.722291
22 静岡市 34.976978 138.383054
23 名古屋市 35.180188 136.906565
24 津市 34.730283 136.508591
25 大津市 35.004531 135.86859
26 京都市 35.021004 135.755608
27 大阪市 34.686316 135.519711
28 神戸市 34.691279 135.183025
29 奈良市 34.685333 135.832744
30 和歌山市 34.226034 135.167506
31 鳥取市 35.503869 134.237672
32 松江市 35.472297 133.050499
33 岡山市 34.661772 133.934675
34 広島市 34.39656 132.459622
35 山口市 34.186121 131.4705
36 徳島市 34.06577 134.559303
37 高松市 34.340149 134.043444
38 松山市 33.84166 132.765362
39 高知市 33.559705 133.53108
40 福岡市 33.606785 130.418314
41 佐賀市 33.249367 130.298822
42 長崎市 32.744839 129.873756
43 熊本市 32.789828 130.741667
44 大分市 33.238194 131.612591
45 宮崎市 31.91109 131.423855
46 鹿児島市 31.560148 130.557981
47 那覇市 26.212401 127.680932
;
run;

県庁所在地と緯度経度をネットから適当に拾ってきて
geodist関数にあてます。
geodist関数は(地点1の緯度,地点1の経度,地点2の緯度,地点2の経度)で、距離を返します(デフォルトなら単位はキロメートル、オプションでマイルにもできます)

ちなみに緯度経度には日本測地系と世界測地系があって、
微妙に違うので、できれば世界に変換して同じ測地系で比べましょう。

変換法はPHPだけど式は同じなので以下のページで紹介されているコードなんかがいいと
思います
http://d.hatena.ne.jp/nakamura001/20080501/1209660263


で、書いたコードは以下のとおり、これぐらいのオブザベーション数なら
直積作ってもよかったですけどハッシュでもわりとあっさり書けました。

data a1;
length _shi $10. _ido _keido 8.;
if _N_=1 then do;
call missing(_shi,_ido,_keido);
declare hash h1(dataset:'q1(rename=(shi=_shi ido=_ido keido=_keido))');
h1.definekey('no');
h1.definedata('_shi','_ido','_keido');
h1.definedone();
end;
set q1;
do no= 1 to 47;
rc=h1.find();
dist=geodist(ido,keido,_ido,_keido);
if dist^= 0 and (dist<mindist or mindist=.) then do;
mindist=dist;
moyori=_shi;
end;
end;
keep shi moyori mindist;
run;

proc sort data=A1;
by mindist;
run;
proc print;
run;

結果は



















































































滋賀の大津市と京都の京都市か~、確かに近いな!
10キロちょっとか

神戸にとっての最寄は大阪だけど大阪の最寄は奈良か。
へ~、そういう感じだったんですね。
勉強になりました。

詰めSAS データステップで新しく作成された変数か、元からデータセットに あった変数かを判定する方法を考えてみる。vinarray関数とか

結構、よく質問されるのが、データステップ中で新しく追加された変数か、
読み込んだデータセットに元からあった変数どうかを判別したいみたいな話です。

いちいち、データセットの仕様書つくるの面倒くさいから
先にプログラム書いちゃって、その実行ついでに新規変数がどうやって生成されているかを
取得して、仕様書を自動生成しちゃえば楽だよなぁみたいな、おバカなこと考えたことありません?
え、いや、僕はないですけどね。

ところが、SASの場合、これが結構ハードで厳しい。
手抜きするつもりが、真面目にやる1000倍くらい時間がかかったりするのもよくある話です

しかも、いまだに完全に納得できる方法が思いついてません
誰か教えてください。

とりあえず、前後のデータセットをコンペアしてその結果を読み込むという
面倒すぎる方法は除外して、1ステップ内で、本処理のついでにできる方法に限定します。

僕が考えた案を2つほど。

以下のデータセットがあります。



変数 a b c dがあります。

data Q1;
a=1;b=2;c=1;d='A';
run;

まず一つ目の方法は以下です。

data A1;
set Q1;
array an _numeric_;
array cn _character_;

/*新しい変数 e f g 作る*/
e=1;
f='b';
g=a+b;

/*チェックしてみる*/
flg_a=ifc(vinarray(a),'元','新');
flg_b=ifc(vinarray(b),'元','新');
flg_c=ifc(vinarray(c),'元','新');
flg_d=ifc(vinarray(d),'元','新');
flg_e=ifc(vinarray(e),'元','新');
flg_f=ifc(vinarray(f),'元','新');
flg_g=ifc(vinarray(g),'元','新');

put (flg_a--flg_g) (/=);

run;

判定結果をputしているので、それを見てみると











追加した3つの変数が正しく判定されているのがわかります。

なにをしているかというと、最初に数値型、文字型それぞれ全変数を
配列に定義しちゃいます。

で、その後、割り当てステートメントがはいるわけですが、
ここで新規に割り付けられた変数は、最初に定義している配列には含まれていない
わけです。

そこで、判定の箇所でvinarray関数という、引数の変数が何かしらの配列に属しているか
どうかを1 , 0で返すニッチな関数で判定します。
配列に入ってないということは、新しい変数だろってノリで。

しかしまあ、お気づきの方もいるでしょうが、これは新しく作る変数を配列に定義しているような
プログラムでは正しく判定できない限定的な方法ですね。
e f gが何の配列にもたまたま定義されていないからそれでいけてるわけです。
同じく、新規変数をarrayで生成したりしても判別不能になります。
(他にも弱点はありますがそれは後で)

次に思いついたのはハッシュオブジェクトです。

data A2;
set Q1;

/*この時点での全変数の情報を一旦ハッシュオブジェクトにいれる*/
if _N_=1 then do;
 declare hash h1();
 h1.definekey('vname');
 h1.definedone();
 array an _numeric_;
 array cn _character_;
 do over an;
  vname=vname(an);
h1.add();
 end;
 do over cn;
  vname=vname(cn);
h1.add();
 end;
 drop vname;
end;

/*新しい変数 e f g 作る*/
e=1;
f='b';
g=a+b;

/*チェックしてみる*/
flg_a=ifc(h1.check(key:vname(a)),'新','元');
flg_b=ifc(h1.check(key:vname(b)),'新','元');
flg_c=ifc(h1.check(key:vname(c)),'新','元');
flg_d=ifc(h1.check(key:vname(d)),'新','元');
flg_e=ifc(h1.check(key:vname(e)),'新','元');
flg_f=ifc(h1.check(key:vname(f)),'新','元');
flg_g=ifc(h1.check(key:vname(g)),'新','元');

put (flg_a--flg_g) (/=);

run;

結果は同じなので割愛。

最初のレコードの時点のPDVの変数名情報をハッシュオブジェクトに格納して
ステップの終わりで、checkメソッドで、存在するかどうかを判定する方法
(さっきと逆に、あれば0,なければエラーコードとして0以外の数字が戻ります)
これならば、ステップの途中で新規変数に配列つかったりしてても問題なしです。


ただ、2つの方法とも、致命的に弱点なのがあくまでステップ開始時の
PDVにある変数を、元からある変数だと考えている点です。

つまり length f $100.;といったようにlengthやformatステートメントなど
データを読み込む前に箱を定義しちゃうようなものがはいると、fは実際は元データセットに無いけど
元からある変数ってことになるんですね。
ここが僕の限界。

う~ん、どうするのが正解なんでしょう。
contentsプロシジャとかsqlとかで、ガチで大元のデータセット情報読み込むコードを
作ってステップ内でdosubl関数で回すとか?

あるいは新規データセットのPDVに影響受けずに変数名をとれるのかな?
if _N_=0の状態でvnextルーチン回したりできんのかな?

誰か教えてください