k-db から hesonogoma で大騒ぎ

　２０１７年末を持ってk-dbが閉鎖ということになりました。お世話になりました。代替サイトですが、hesonogomaさんの株価一覧表を使うこととしました。証券会社で一括データの提供くらいしてくれてもよさそうなものですが、調べた限りでは見つかりませんでした。で、いまだ、対応途上ですが、目鼻が付きましたので、それに伴う様々なトラブルをメモしておきます。

1. データのダウンロード・・・k-dbではurlの指定でcsvがダウンロードできましたが、hesonogomaさんではフラッシュを使用しているようで、そのようなわけにはいかず、手作業でファイルをダウンロードすることとしました。
2. データファイルが４つ必要・・・全銘柄・ETF等・REIT・ファンドと４種類のファイルが必要です。
3. 文字コードの問題１・・・文字コードの問題は相当重症です。ダウンロードしたファイルの文字コードは「UTF-16LE BOM付き」でした。UTF-16をRで読み込むとファイルが途中で切れます。当初原因がわからず、試行錯誤しましたが、StackOverflowに中国人の方が、同じ現象に会って質問をしていました。特定のコード以降読み込まないという趣旨だろうと思います。回答はありませんでした。で、nkf.exeを使っていったんUTF-8に変換することにしました。
4. 文字コードの問題２・・・最初「UTF-8 BOMなし」に変換したのですが、EXCELに読み込むと文字化けする。で、BOM付きに変換。そのせいか、そのファイルをRで読み込むと列数が多いといわれて読み込めない。

ということで現在のところの文字コード部分のコードは次の通り。

system(paste(nkf,"-w8 --overwrite",fnams[1]))
system(paste(nkf,"-w8 --overwrite",fnams[2]))
system(paste(nkf,"-w8 --overwrite",fnams[3]))
system(paste(nkf,"-w8 --overwrite",fnams[4]))
d1=read.csv(fnams[1],sep="\t",na.strings = "-",fileEncoding="UTF-8-BOM",stringsAsFactors=F)
d2=read.csv(fnams[2],sep="\t",na.strings = "-",fileEncoding="UTF-8-BOM",stringsAsFactors=F)
d3=read.csv(fnams[3],sep="\t",na.strings = "-",fileEncoding="UTF-8-BOM",stringsAsFactors=F)
d4=read.csv(fnams[4],sep="\t",na.strings = "-",fileEncoding="UTF-8-BOM",stringsAsFactors=F)

なお、「UTF-8 BOM付き」をEXCELで読み込むと、tabで区切ってくれないのでもうひと手間かかります。

また、fileEncodingに「UTF-8-BOM」はありますが「UTF-16-BOM」はありませんでした。「UTF-16」はありますが上記の問題に会います。

さらに、４つのファイルの項目名が違う。全銘柄ファイルでは業種、それ以外のファイルでは種別となっている。

その部分を乗り越え、一つのファイルにするのが次のコード

names(d2)[which(names(d2)=="種別")]="業種" #names(d2)=gsub("種別","業種",names(d2))のほうがおしゃれかな？
names(d3)[which(names(d3)=="種別")]="業種"
names(d4)[which(names(d4)=="種別")]="業種"
orinames=c("SC","名称","市場","業種","始値","高値","安値","株価","出来高","売買代金.千円.")
d=rbind(d1[,orinames],d2[,orinames],d3[,orinames],d4[,orinames])

とりあえず、一つのデータフレームにできました。次の問題は

1. 銘柄名が長い。・・・新しい銘柄名に変えようか、迷っている。
2. 市場名が違う。・・・これは古いものに統一しようと思っている。
3. 銘柄コードが違う。・・・k-dbでは「1301-T」のような形式でしたが新しいほうは「1301」のような形式です。新しい形式にしたときに問題になりそうなことは、
• データ参照の時にyt[["1301-T"]]という参照方法とyt[[1]]という参照方法がありましたが、「1301」形式では文字列として扱われるのか整数として扱われるのか、悩ましい。Rが勝手に変換してくれるのもよいのやら。
• k-dbでは同一コード同一銘柄で別市場のものがありましたが、新しいほうでは同一コードが見当たりません。東証と別市場がある場合は東証のみになっているのでしょうか？

ということで、大晦日から元旦にかけての作業でした。

R ディレクトリのみ表示

ディレクトリのみ表示する。

t=dir()
t[file.info(t)$isdir]

１文にするなら
dir()[file.info(dir())$isdir]

こんな命令があった。
list.dirs(".",recursive = F)

参考https://stackoverflow.com/questions/4749783/how-to-obtain-a-list-of-directories-within-a-directory-like-list-files-but-i

R 関数を文字列で呼びだす

do.call("mean",list(1:10,na.rm=T))
第２引数は必ずリストで。これが関数への引数リストとなる。

eval((parse(text= "mean(1:10,na.rm=T)")))
関数に限らず文字列を実行したいとき。

R 練習問題

Rの練習問題を作ってみました。今回の対象は次のようなデータフレームです。試験の点数処理を題材としてみました。データを作るスクリプトです。乱数でとも思いましたが、同じ結果の方が見やすいでしょうから。

データはこんな風です。３組それぞれ５人の国・数・英の点数です。

問題と出力です。
#問１　dに個人ごとの計を付加

#問２　計順にソートして表示

#問３　dに計で学年順位を付加

#問４　組ごと計順に表示

#問５　教科間の相関係数

#問６　科目別箱髭図

#問７　国の組別箱髭図

#問８　教科集計

#問９　計の度数分布

#問１０　科目別平均

#問１１　計の組別平均

#問１２　科目別組別平均

#問１３　国の組別平均・最高・最低

#問１４　国50未満を抽出

#問１５　国50未満を抽出、組番国のみ表示

#問１６　数40以上英50以上を抽出

#問１７　数30以上を計順に表示

#問１８　科目別、成績順、組番

#問１９　計の２０点刻み度数分布図

#問２０　科目別５点刻み度数分布表


解答


終わりに
この練習問題では科目は横に展開し、組は縦に展開しているので、科目ごとと組ごとの扱いが違って来る。
apply系の命令をうまく使うことは難しい。この練習問題では、結構面白い使い方ができたと自負している。

簡単にポイントにふれておく。

問１　①新たな列の付けくわえ方②d[3:5]はデータフレームdの３列目から５列目に限定している。データフレームはベクトルのリストであるので、d[3:5]で３列目から５列目が取り出される。これがmatrixならば、d[,3:5]とカンマとつけなければならない。もちろんデータフレームでもこのように指定してもよいが。③行集計にrowSumsを用いた。apply(d[3:5],1,sum)としてもよい。

問２　①ソートにはorderを用いるのが良い。単一のベクトルをソートするにはsortを用いればよいが、キーでソートするようなときはorderがよい。orderは昇順にソートするので、高得点の方が先に来るように、マイナスした。desc=Tを用いてもよいが、おしゃれじゃない。②行番号で指定するので、d[order(-d$計),]と、コンマが必要となる。問１で述べたように、ここでコンマがないと、列の指定となってしまう。

問３　①順位はrank。②引数のt=は ties.method = c("average", "first", "random", "max", "min")である。引数の変数名は誤解のない範囲で省略できる。

問４　orderでソート項目は並べればよい。

問５　特にない

問６　これだけで科目ごとの箱髭図を書いてくれるなんて、なんて便利なんだろう。

問７　組ごとの箱髭図はこう指定する。問６との違いは、先に述べたように、科目は横に展開し、組は縦に展開していること。

問８　特にない

問９　特にない

問１０　d[3:5]がデータフレームなので、sapplyはリストの要素の各ベクトルすなわち各列に対して動作する。

問１１　問１０との違いは、前述のとおり

問１２　applyの中でtapplyを使ってみた。うまくいった。sapplyでもいいのだね。

問１３　いろいろと試行錯誤し、もっともよさげなものをあげた。

問１４～１７　抽出

問１７　抽出と整列を同時にできないかというチャレンジ。一時変数を使わないでやりたかったが、最終的には用いてしまった。敗北感あり。使わない版もいくつか作ったが、同じ計算を２度やらせていて、許せない。いつの日か思いつくかもしれない。

問１８　科目別のベストテンなどの表はよく作るだろう。それである。文字列ベクトルを作っておいて、その並び順を添え字で指定するというこのパターンは初めての経験かもしれない。

問１９　breaksについては問２０参照

問２０　科目別度数分布表。苦労したのは、区間の指定。たとえば、breaks=c(10,20,30,40)とすると、この４つの数値によって、３つの区間に分けられる。10≦x≦20,20＜x≦30,30＜x≦40。最初の10の不等号に＝が入っている

サンプルデータでは、回帰係数や回帰分析や、主成分分析、クラスタ分析も面白くなかろうと思ってやめた。

移動平均について

移動平均について考察してみる。なお、単純な移動平均についてである。
本日の終値を C0、前日の終値をC1、前々日の終値をC2、…と書く。
騰落率は C0/C1-1である。
５日移動平均MA5=(C0+C1+C2+C3+C4)/5である。

移動平均からの乖離率について

今、１日の騰落率が５日間ともrであったとしよう。すると、
C0 = C4×(1+r)^4 ≒ C4×(1+4r)
ｒが十分0に近いとき、(1+r)^n ≒ 1+nr であることを用いている。 今後これをどしどし使う。

同様に、C1≒C4×(1+3r)、C2≒C4×(1+2r)、C3≒C4×(1+r)、であるから、
MA5 ≒ C4×{(1+4r)+(1+3r)+(1+2r)+(1+r)+1}/5 ≒ C4×(1+2r) となる。
C0のMA5からの乖離率は、
C0/MA5-1 ≒ (1+4r)/(1+2r)-1
　　　　　　　≒ (1+4r)(1-2r)-1
　　　　　　　≒ 2r
すなわち、騰落率の２倍の値となる。逆に言えば、乖離率の２分の１が１日当たりの騰落率ということになる。
同様に、（２５日移動平均からの乖離率）÷１２ ＝ ２５日間の平均騰落率
　　　　　（７５日移動平均からの乖離率）÷３７ ＝ ７５日間の平均騰落率
と考えることができる。

移動平均曲線の傾き

移動平均が右肩上がりならば、上昇トレンドにあると考えられる。しかし、チャートを比べて上昇の度合いを考えたいと思っても、チャートの縦の比率が一定ではないので、単純には比べられない。
移動平均曲線の傾きはどう求めたらよいだろうか。前述のMA5で調べてみる。

最終のMA5の傾きは、前日MA5と当日MA5の差である。
(MA5の傾き)＝(当日MA5)-(前日MA5)
　　　　　　　　＝ (C0+C1+C2+C3+C4)/5-(C1+C2+C3+C4+C5)/5
　　　　　　　　＝ (C0-C5)/5

 なんと、５日前との終値の差の５分の１である。途中の値は関係ないのだ。
C0-C5はモメンタムというらしい。
同様に、２５日移動平均の最終日の傾きは、２５日前との終値の差の２５分の１となる。
ただの終値の差に、こんな意味があったとは！！

さて、実際には、株価が高いところでは当然その傾きも大きくなるのだから、比較のためには、
C5/C0-1の値を使うことになるであろう。この値はROC(rate of change)というらしい。

quantmodでは momentum(x, n = 1, na.pad = TRUE)
ROC(x, n = 1, type = c("continuous", "discrete", na.pad = TRUE)
となっている。上の計算と合わせるなら"discrete"である。

数独を解く２

１０月２８日付けブログの数独を解くプログラムで、数独を解いてみた。連続的に解くために若干の修正をしたが、本質的には変わりない。

サンプルとして挙げた「超極・激」には、対角線指定や、変形３×４などが含まれており、それらを除いた通常のもの１３０題を解いたものである。

出力の冒頭を以下に示す。



問題と解答と所要時間を示した。

解答時間の平均は０．６９２秒　最大３．４４秒であった。使用マシンのCPUはi7である。解答時間の分布は



詳細を検討してみると、checkだけで解けているものも多数あり、solveを用いて試行錯誤的にアプローチしているものも多い。solveを用いた場合、再帰の深さは結構深いのが目につく。

前回の解法アルゴリズムで、１つのグループ内に例えば５の可能性のあるマスが１つだけならばそこが５である、というアルゴリズムを最後に付け加えた。

これを拡張して、１つのグループ内に可能性が５と８だけであるマスが２つあれば、そのグループのほかのマスからは５と８の可能性が削除できる。

また、１つのグループ内に　可能性５と８を含むますが２つだけならば、その２つのマスの５と８以外は削除できる。

さらに拡張すると、１つのグループ内に例えば可能性が５と８と９だけであるマスが３つあれば、そのグループのほかのマスからは５と８と９の可能性が削除できる。この場合、この３つのマスのうちの１つは「５と８」でもよい。など。

これらについてはまだ、実装してない。

数独を解く (sudoku ナンプレとも)

Rで数独を解いてみた。その中で、R独特の機能が結構有効であったので、解説してみる。サンプルとして、Take-Shobo「ナンプレ超極・激」の３０番を使用した。次のような問題である。

結果は

整形してないので見にくいが、９文字ごとに改行を入れると答えになる。もちろん先頭の[1]などは除いて。

Rのプログラムというのだろうか、ソースは以下のようなものである。

１行目、qが問題である。８１文字の数字で、０は空欄を表す。q　すなわち　question である。

２行目、qを１文字ずつに分解し、更に数値に変換している。０～９の数字８１個のベクトルである。名前bはボードの意である。
プログラムの中で、bは重要な位置を占める。
数独の問題は、視覚的には平面であるが、プログラム内では１次元ベクトルとして表現している。実は、この方がよいのである。そのことはだんだんと分かってくる。
数独のマスに順に１～８１までの番号を付け、この番号を「座標」と呼ぶことにする。座標はiであらわす。

３行目から１０行めはグループの指定である。
数独では各１行の中に１～９のすべての数字が入る。各１列の中にも同様である。さらに、３×３に分けた９個の各正方形の中においても同様である。この、行・列・３×３正方形の１つずつを「グループ」と呼ぶことにする。gsは各グループの座標を持つ行列である。gsはグループスの意味である。内容は次のようになる。

[1,]は１行目を表す。[10,]は１列目を表す。[19,]は左上の３×３正方形を表す。
gsを用いて、各グループに同じ数字がないことを確認するのである。実は思いがけない利用法がある、対角線上に１～９があることを条件としたり、変な形の領域に１～９があることを条件にする問題があるが、ここを変えればよいのである。
グループは27個あり、その番号はkで指定する。また、一つのグループはベクトルgで表す。

１１行めから４４行めの関数checkが肝である。
引数が複雑だが、外から使うときはcheck(b)のように使えばよい。このように使うと、１２行めで行列sbを作る。
このsbがこのプログラムのもう一つの重要点である。サブボードのつもりでsbという名前を付けた。
ボードの各マスには１～９のいずれかの数字が入るのだが、その候補を記録するのがsbである。
sbは８１行９列の行列で、各行が、各座標に入る数字の候補を記録している。
sb[i,n]は、座標iに数字nが入る可能性があれば1、なければ0をとる。
例えば１行目が 1,1,0,0,0,0,1,1,1 となっているとすると、座標１のところには、1,2,7,8,9のいずれかが入る。逆に言えば3,4,5,6は入らないことを意味する。
sbの候補をどんどん減らしていって答えを見つけるのである。
最初はすべて１である。すなわち、すべての数字の入る可能性がある。

どのように候補を絞っていくかが、そのあとに続く。

１５・１６行めは、ボードに数字が入っている場合、すなわち０でない場合は、その座標のsbは、その数字のところを１とし、それ以外を０とする。

b[i]は座標iの位置にある数字である。「sb[i,-b[i]]=0」はsbのi行めのうちの、「b[i]すなわち座標iの数字」以外のところを０にする。普通の言語だとループでやるところだ。こんなに簡単にできる。（簡単というのは表現の文字数が少ないという意味だ。表現の理解はそう簡単ではない）

１８行めから２６行め。for文を使ってしまった。applyか何かでやりたいところだが、今後の課題。
ここでは、ボード上で数字が入っていたら、その座標を含むグループのほかの座標ではその数字は使えないので、sbのその位置を０にする。もっと具体的に言うと、座標iの位置にある数字がnのとき、座標iを含むグループの他の座標jについて、sb[j,n]を0にするのだ。
プログラムに沿って解説する。
kはグループ番号で、１～２７まで変わる。
k番目のグループをgとする。gは長さ９のベクトルで、k番目のグループの座標を持っている。
pはg内をポイントし、１～９変わる。
g[p]がグループgのp番目の座標なので、その位置の数字はn=b[g[p]]となる。nが0の時は数字はないので何もしない。0でないときは座標g[p]に数字nがあるので、そのグループ内のほかの座標のsbのnのところを0にしなければならない。なんとそれが、
sb[g[-p],n]=0でできてしまう。g[-p]はグループgのポイントしている以外のすべての座標からなるベクトルだ。それらの行のn番目を0にするのだ。

ここまではbの情報をもとにsbを変更した。

次の２７行めから３０行めは、sbをもとに、bを変更する部分である。

iはボード上の座標である。１～８１の値をとる。
sb[i,]は座標iのサブボード情報である。すなわち、各数字が可能性があるかどうかを示す９つの数（０か１）からなる。
sum(sb[i,])が0だということは、その座標に入る数字がないことになる。すなわち、解（甲斐）なしである。NULLを返すこととした。
sum(sb[i,])が１だということは、座標iに入る数字が一つに決まったということである。その場合はボードに数字を書き込む。

１５行めから３０行目までを繰り返すと、易しい問題は解ける。繰り返しは必要で、sbからbに変更が加わるとそのbからまたsbに変更が生し、・・・ということがあるからである。１３行めから４２行目までがrepeatで囲まれているのはそのためである。
繰り返して、変更がなくなったらrepeatを抜けるのだが、その判断は、sbのフラグの総数の変化で調べている。sbのフラグ（１）の総数は単調減少だ。

３１行めから４０行めの説明は後回しにする。

４３行め、返す値は、ボードbとサブボードsbのリストである。

checkだけで解ける問題は易しい問題である。checkだけでは解けない問題を試行錯誤で解くようにしたのがsolveである。
solveは最初の空欄に１～９までの数字を当てはめてcheckを繰り返す。１段階の当てはめでうまくいかない場合は、再帰的に多段階で挑戦させる。

solveの引数はl(エル)である。checkで返すリストを引数にする。リストなのでlだ。l[[1]]はボード、l[[2]]はサブボードだ。サブボードは渡さなくても計算できるが、せっかくそれまで計算したものをはじめから計算するのもどうかと思って、引数として渡すことにした。

さて、最後に残った３１行めから４０行めだ。
当初、これはつけないつもりでいた。試行錯誤させればよいと考えていたからだ。ところがサンプルの問題では１時間ほど動かしても答えが出ない。空欄が多い問題ではこうなるのは必然だ。

ここでやっていることは、１つのグループ内のsbのフラグをすべて見て、例えば５のフラグが１個しか立っていないならば、そこが５だ
ということを調べている。

このルーチンを付け加えたら、サンプル問題は一瞬で解けた。試行錯誤なしで、checkのみで解けた。「超極・激」というには情けない。まあ、３０番だから。このあと難しいのがあるのだと思う。

このように実際の数独を解くときに使う戦略を組み込むと速度は劇的に改善する。今回Rで作ってみたが、以前に一度python3でプログラムしたことがある。その時の経験に基づいた発言である。

以上でこの稿は終わる。
しばらく、更新してなかったので、更新しなければという強迫観念で書いた。
きっと読みにくいだろうと思う。
面白いと思ってくれる人がいるのだろうか？