Rでペンローズタイルを描く
この記事は R Advent Calendar 2014 の25日目の記事です。
ペンローズ・タイルって皆様ご存知でしょうか。
ペンローズ・タイル - Wikipedia
簡単にまとめると以下のような特徴を持つ不思議なタイルとその平面充填方法です。
- 非周期的であり、同じ充填方法が現れない。
- kiteとdirtと呼ばれる二つの四角形から構成される。
- それぞれの四角形は同じ形のロビンソンの三角形という二等辺三角形を2つ組み合わせたもので出来ている。
- 三角形の長辺との長さは黄金比になっている。
- 再帰的に分割することができ、フラクタル図形の一種と言える。
このペンローズタイルと同じ性質を持つ図形はイスラム建築の模様や、電動シェーバーの網刃、準結晶等様々なところに見られます。
とっても不思議な図形ですね。まとまった解説記事がないのが残念です。
ということで、このペンローズタイルをRで描いてみました。せっかくなのでクリスマスツリー状にしています。
PHI<<-(1+sqrt(5))/2
L<<-2
decompose_kite<-function(t1,k,rev=FALSE){
if(rev){
x1<-c((t1$coord$x[2]/PHI+t1$coord$x[1])/PHI,t1$coord$x[3],t1$coord$x[1])
y1<-c((t1$coord$y[2]/PHI+t1$coord$y[1])/PHI,t1$coord$y[3],t1$coord$y[1])
x2<-c((t1$coord$x[2]/PHI+t1$coord$x[1])/PHI,t1$coord$x[3],(t1$coord$x[2]+t1$coord$x[3]/PHI)/PHI)
y2<-c((t1$coord$y[2]/PHI+t1$coord$y[1])/PHI,t1$coord$y[3],(t1$coord$y[2]+t1$coord$y[3]/PHI)/PHI)
x3<-c(t1$coord$x[2],(t1$coord$x[2]+t1$coord$x[3]/PHI)/PHI,(t1$coord$x[2]/PHI+t1$coord$x[1])/PHI)
y3<-c(t1$coord$y[2],(t1$coord$y[2]+t1$coord$y[3]/PHI)/PHI,(t1$coord$y[2]/PHI+t1$coord$y[1])/PHI)
}else{
x1<-c(t1$coord$x[3],t1$coord$x[1],(t1$coord$x[2]/PHI+t1$coord$x[3])/PHI)
y1<-c(t1$coord$y[3],t1$coord$y[1],(t1$coord$y[2]/PHI+t1$coord$y[3])/PHI)
x2<-c((t1$coord$x[2]/PHI+t1$coord$x[3])/PHI,t1$coord$x[1],(t1$coord$x[2]+t1$coord$x[1]/PHI)/PHI)
y2<-c((t1$coord$y[2]/PHI+t1$coord$y[3])/PHI,t1$coord$y[1],(t1$coord$y[2]+t1$coord$y[1]/PHI)/PHI)
x3<-c(t1$coord$x[2],(t1$coord$x[2]+t1$coord$x[1]/PHI)/PHI,(t1$coord$x[2]/PHI+t1$coord$x[3])/PHI)
y3<-c(t1$coord$y[2],(t1$coord$y[2]+t1$coord$y[1]/PHI)/PHI,(t1$coord$y[2]/PHI+t1$coord$y[3])/PHI)
}
small_coord <-list(x=x1,y=y1)
small_coord2<-list(x=x2,y=y2)
big_coord <-list(x=x3,y=y3)
polygon(small_coord$x, small_coord$y, col = "green")
polygon(small_coord2$x, small_coord2$y, col = "green")
polygon(big_coord$x, big_coord$y, col = "red")
polygon(-small_coord$x, small_coord$y, col = "green")
polygon(-small_coord2$x, small_coord2$y, col = "green")
polygon(-big_coord$x, big_coord$y, col = "red")
if(k<L && rev==FALSE){
decompose_dart( list(coord=big_coord), k=k+1, rev)
decompose_kite(list(coord=small_coord), k=k+1,!rev)
decompose_kite(list(coord=small_coord2), k=k+1, rev)
}
if(k<L && rev==TRUE){
decompose_dart( list(coord=big_coord), k=k+1,rev)
decompose_kite(list(coord=small_coord), k=k+1,rev)
decompose_kite(list(coord=small_coord2), k=k+1,rev)
}
return(list(small=list(coord=small_coord), big=list(coord=big_coord)))
}
plot_Tri_kite<-function(a,b,r,theta,rev=FALSE){
x<-c(0,-1*cos(pi*3/5),-2*cos(pi*3/5))
y<-c(0,sin(pi*3/5),0)
x1<-(r*x*cos(theta)-r*y*sin(theta))+a
y1<-(r*x*sin(theta)+r*y*cos(theta))+b
coord<-list(x=x1,y=y1)
decompose_kite(list(coord=coord),1,rev)
return(list(coord=coord))
}
decompose_dart<-function(t1,k,rev=FALSE){
if(rev){
x1<-c((t1$coord$x[3]+t1$coord$x[1]/PHI)/PHI,t1$coord$x[1],t1$coord$x[2])
y1<-c((t1$coord$y[3]+t1$coord$y[1]/PHI)/PHI,t1$coord$y[1],t1$coord$y[2])
x2<-c(t1$coord$x[3],(t1$coord$x[3]+t1$coord$x[1]/PHI)/PHI,t1$coord$x[2])
y2<-c(t1$coord$y[3],(t1$coord$y[3]+t1$coord$y[1]/PHI)/PHI,t1$coord$y[2])
}else{
x1<-c(t1$coord$x[2],t1$coord$x[3],(t1$coord$x[3]/PHI+t1$coord$x[1])/PHI)
y1<-c(t1$coord$y[2],t1$coord$y[3],(t1$coord$y[3]/PHI+t1$coord$y[1])/PHI)
x2<-c(t1$coord$x[2],(t1$coord$x[3]/PHI+t1$coord$x[1])/PHI,t1$coord$x[1])
y2<-c(t1$coord$y[2],(t1$coord$y[3]/PHI+t1$coord$y[1])/PHI,t1$coord$y[1])
}
small_coord<-list(x=x1,y=y1)
big_coord<-list(x=x2,y=y2)
polygon(small_coord$x, small_coord$y, col = "green")
polygon(big_coord$x, big_coord$y, col = "red")
polygon(-small_coord$x, small_coord$y, col = "green")
polygon(-big_coord$x, big_coord$y, col = "red")
if(k<L){
decompose_dart(list(coord=big_coord), k=k+1)
decompose_kite(list(coord=small_coord), k=k+1)
}
return(list(small=list(coord=small_coord), big=list(coord=big_coord)))
}
plot(0,0, xlim=c(-10,10), ylim=c(-10,10))
t1<-plot_Tri_dart(0,0,5,0)
plot_Tri_dart<-function(a,b,r,theta,rev=FALSE){
x<-c(0,cos(pi*1/5),2*cos(pi*1/5))
y<-c(0,sin(pi*1/5),0)
x1<-(r*x*cos(theta)-r*y*sin(theta))+a
y1<-(r*x*sin(theta)+r*y*cos(theta))+b
if(rev){coord<-list(x=x1[3:1],y=y1[3:1])}
else{coord<-list(x=x1,y=y1)}
decompose_dart(list(coord=coord),1)
return(list(coord=coord))
}
plot(0,0, xlim=c(-5,5), ylim=c(-1,9))
t1<-plot_Tri_kite(0,0,1,pi/10)
t1<-plot_Tri_dart(t1$coord$x[2],t1$coord$y[2],1/PHI,-3*pi/10)
t2<-t1
t1<-plot_Tri_kite(t1$coord$x[2],t1$coord$y[2],1,-5*pi/10)
t2<-plot_Tri_kite(t2$coord$x[1],t2$coord$y[1]+1/PHI,1,-5*pi/10)
t3<-t2
t2<-plot_Tri_dart(t2$coord$x[2],t2$coord$y[2],1/PHI,-9*pi/10)
t2<-plot_Tri_kite(t2$coord$x[1],t2$coord$y[1],1,-7*pi/10)
t3<-plot_Tri_dart(t3$coord$x[1],t3$coord$y[1]+1,1/PHI,-5*pi/10)
t4<-t3
t3<-plot_Tri_kite(t3$coord$x[2],t3$coord$y[2],1,-7*pi/10)
t3<-plot_Tri_dart(t3$coord$x[1],t3$coord$y[1],1/PHI,-3*pi/10)
t3<-plot_Tri_dart(t3$coord$x[2]+cos(3*pi/10)/PHI,t3$coord$y[2]+sin(3*pi/10)/PHI,1/PHI,-9*pi/10)
t4<-plot_Tri_dart(t4$coord$x[1],t4$coord$y[1],1/PHI,-1*pi/10)
t5<-t4
t4<-plot_Tri_dart(t4$coord$x[3],t4$coord$y[3],1/PHI,5*pi/10)
t4<-plot_Tri_dart(t4$coord$x[1],t4$coord$y[1],1/PHI,9*pi/10)
t5<-plot_Tri_kite(t5$coord$x[1],t5$coord$y[1],1,pi/10)
t5<-plot_Tri_kite(t5$coord$x[3],t5$coord$y[3],1,3*pi/10)
t5<-plot_Tri_dart(t5$coord$x[2],t5$coord$y[2],1/PHI,-1*pi/10)
t6<-t5
t5<-plot_Tri_dart(t5$coord$x[3],t5$coord$y[3],1/PHI,5*pi/10)
t5<-plot_Tri_kite(t5$coord$x[1],t5$coord$y[1],1,pi/10)
t5<-plot_Tri_kite(t5$coord$x[3],t5$coord$y[3],1,3*pi/10)
t5<-t6
t5<-plot_Tri_kite(t5$coord$x[1],t5$coord$y[1],1,pi/10)
t5<-plot_Tri_kite(t5$coord$x[3],t5$coord$y[3],1,3*pi/10)
t5<-plot_Tri_dart(t5$coord$x[2],t5$coord$y[2],1/PHI,-1*pi/10)
t6<-t5
t5<-plot_Tri_dart(t5$coord$x[3],t5$coord$y[3],1/PHI,5*pi/10)
t5<-plot_Tri_kite(t5$coord$x[1],t5$coord$y[1],1,pi/10)
t5<-plot_Tri_kite(t5$coord$x[3],t5$coord$y[3],1,3*pi/10)
t6<-plot_Tri_kite(t6$coord$x[1],t6$coord$y[1],1,pi/10)
t6<-plot_Tri_kite(t6$coord$x[3],t6$coord$y[3],1,3*pi/10)
t6<-plot_Tri_dart(t6$coord$x[2],t6$coord$y[2],1/PHI,-1*pi/10)
t6<-plot_Tri_dart(t6$coord$x[1],t6$coord$y[1],1/PHI,3*pi/10)
t7<-plot_Tri_dart(t6$coord$x[3],t6$coord$y[3],1/PHI,-7*pi/10)
t6<-plot_Tri_kite(t6$coord$x[2],t6$coord$y[2],1,pi/10)
実行結果が以下の通り、
ところどころロビンソンの三角形がむき出しになっていますね。。。
失敗したようなので来年の宿題としてリベンジしたいところです。
Numcraw for R User
Tomorrow, Japan.R will be held.
It is one of the biggest R user meetup.
If you don't know about it, I recommend you to google it.
Unfortunately, I can not join Japan.R because of a business trip.
So following slide is to celebrate Japan.R from U.S.
The function Numacraw() in the slide is here.
numacraw<-function(get2=FALSE){
if(require(png) && require(RCurl)){
set.seed(Sys.time())
path_back<-paste0(tempdir(),"\\back.png")
download.file("http://pbs.twimg.com/media/B4AEO4LCEAMjpx7.png"
, destfile = path_back,mode="wb")
path_front<-paste0(tempdir(),"\\front.png")
download.file("http://pbs.twimg.com/media/B4AB8oDCIAA_04r.png"
, destfile = path_front,mode="wb")
vec<-par()$usr
SIZE<-runif(min=0.4,max=1, 1)
X<-runif(min=0, max=1, 1)*(vec[2]-vec[1])
Y<-runif(min=0, max=1, 1)*(vec[4]-vec[3])
D<-ifelse(runif(min = 0 ,1)>0.3, path_front, path_back)
if(!get2){
rasterImage(readPNG(D)
,X-SIZE*(vec[2]-vec[1])/2
,Y-SIZE*(vec[4]-vec[3])/2
,X+SIZE*(vec[2]-vec[1])/2
,Y+SIZE*(vec[4]-vec[3])/2
)
}else{
rasterImage(readPNG(D)
,X-SIZE*(vec[2]-vec[1])/2
,Y-SIZE*(vec[4]-vec[3])/2
,X+SIZE*(vec[2]-vec[1])/4
,Y+SIZE*(vec[4]-vec[3])/2
)
rasterImage(readPNG(D)
,X-SIZE*(vec[2]-vec[1])/4
,Y-SIZE*(vec[4]-vec[3])/2
,X+SIZE*(vec[2]-vec[1])/2
,Y+SIZE*(vec[4]-vec[3])/2
)
}
}else{
stop("Package png and RCurl are required")
}
}
財務情報を分析した〜試しに財務情報のトレンドでグルーピング
ネット上で企業の決算情報を公開している面白いサイトを見つけたので、前からやりたかった財務情報の分析を三連休最後の日にしてみた。
なお、本記事は本記事やその誤りによって発生した問題や不利益について一切の責任は負いません。念のため。
入手できたデータで出来そうな分析を考えた結果、財務情報のトレンドで企業をグルーピングすることにした。
データ
当該ページより各企業の決算情報を大体2008年くらいから2014年までExcelでダウンロード出来る。
とりあえずなるべく色々な財務情報から企業のデータを分析できないかと考えたが、欠損値が多いので、欠損値の少ない売上高と経常利益を用いることにした。
このデータから「個別」の「日本基準」の情報の中で2008年期から2014年期の7年間分のデータが存在する企業のみを抽出、7年分の売上高と経常利益を取りだし、入力データとした。
これをさらに、各企業の売上高と経常利益それぞれにおいて標準化した。
つまり、企業規模によらない時間によるトレンドを分析の対象としており、似たようなトレンドの企業がグルーピングされる。
入力データのイメージは以下の図のようになっており、Xは標準化された売上高をYは経常利益を表わしている。
分析手法
上記のデータのうち、データクレンジングの結果、利用できる2855企業に関する7年×2種類（売上高・経常利益）の計14データを
単一の個体としてSOMに入力、SOMを用いた理由はあまり面倒なことを考えたくなかったからである。
SOMのセル数は4×5の20にした。RのSOMパッケージを用いて、特にパラメータチューニングは行ってない。
結果
結果のプロットは以下の通り、各セルは類似する特徴をもつ企業グループがnで示された数含まれている。
セルの企業グループの特徴を表わすのが折れ線グラフである。
折れ線の左から7つは2008年から2014年の売上高、残りの7つは同年度別の経常利益を表わしている。
また、隣り合うセル同士には類似した特徴をもつ企業が含まれやすい。
これより、例えばセル0は年度別の変動の少なかった企業があり、セル15には2008年以降から成長した企業が含まれていると推論できる。
セル4は2008年以降業績が悪くなっているので、リーマンショックの影響を受けた企業だろう。
セル19はリーマンショック以降、売上は大きく下がったが、徐々に経常利益を回復している企業である。
各セルに属する企業を見てみた結果、セル番号ごとに以下のような業種の企業が多かった。
（もちろん他の企業もたくさんあった。）
各セルに特徴的な業種
セル0=変動が少ない企業：電力会社
セル4=リーマンショック以降業績が悪くなった企業：鉄鋼、印刷、メガバンク
セル15=順調に成長している：Webサービス、製薬、ガス
セル19=売上高は横ばい、経常利益は伸びている企業:地方銀行、精密機器
一番興味深いのはセル19の企業だと思うんだけど、これって何なんだろうか。
今度、そのあたりに明るい人に聞いてみよう。
ちなみに自動車はこの4つを除く、真ん中のセルに多かった。
結論
たった、2種類のデータでも意外と面白かった。
より細やかな分析を行って、今回は公開を避けた個別の企業名を含めた分析をしてみたい。
色々な財務諸表を使いたいし、時系列データを時系列のデータとして反映した分析を行ってみたい。
Stanで人類最強の男を決定する(1)~モデルを作ってみた。
お久しぶりです久しぶりのblog更新になります。
去る7月12日、第3回BUGS/Stan勉強会が開催されました。
近日中にレビューblogがアップされると思いますが、今回の特徴としては、BUGSに関連する発表がなかったことが挙げられます。
初期のBUGS/Stan勉強会はStanユーザは少なかったことを考えると、この1年くらいでずいぶん状況が変わってきたことはとても感慨深いですね。
今回は総合格闘技の試合データから、世界最強最強の選手を決めるモデルを構築してみた（けど、収束しなかったので、アドバイスあればください）という発表をしました。格闘技好きの人間がいつも議論している話題ですね。
内容に関しては下記のSlideShareの資料をご覧になってください。
とりあえず、収束していないことに関しては、パラメータ数がデータ数に対して多いため、パラメータ数列のバラつきが少ないことが原因かなと思っています。ただし、結果は再現性が高く、ほぼ安定しているので正しい結果が返ってきているようです。
今後はモデルに年齢による選手のキャリアピークの調子を反映させれるように調整したいと思っています。
が、先はパラメータの収束ですね。
追記：Andrwe Gelmanがワールドカップの結果から似たようなモデルを作ってました。参考にします。
http://andrewgelman.com/2014/07/13/stan-analyzes-world-cup-data/?utm_medium=twitter&utm_source=twitterfeed
第2回BUGS/Stan勉強会を開催しました
新年、あけましておめでとうございます。今年もよろしくお願い致します。
さて、昨年末のことで遅くなってしまい大変恐縮ですが、第2回BUGS/Stan勉強会を開催致しました。
BUGS/Stanって何？
前回のBlogの繰り返しではありますが、BUGS/Stanとは何ぞやということを説明しておくとベイズ推定による柔軟なモデリングを実現するためのDSL言語です。
これでモデルを記述すると、非エンジニアなデータ解析屋のみなさんでも（比較的）簡単にデータ分析が出来ます！
このあたりの技術は日本語で枯れていないのでこの勉強会を通して色々と技術を枯らして行きたいと思います。
特にベイズモデリングは恣意性が高く、細かいテクニックと経験が必要なのでそういった悩みのあたりも共有できたらと思っています。
Togetterまとめ
xiangze750さんがTogetterで今回の勉強会に伴うつぶやきを纏めてくださっています。ありがとうございます。
１．「入門セッション(Stanチュートリアル2.0)」（発表者 @TeitoNakagawa）
- インストールについて・・・話題のStan2.0に関してですが、私の環境（Windows7、64bit）ではスムーズにインストール出来ました。普通にrstan1.Xをremove.packageで削除後、通常の手順で入れました。なお、聞いたところによるとMacでのPyStanのビルドに躓いている方が多いようです。こちらのITO Hriokiさんのblog記事が参考になると思います。
- 型変換について・・・質疑応答にありましたがStanではIntegerとRealの型変換ができません。この辺りがStanの使いづらさです。離散的な値などを多用するモデルではJAGS等の方がよいのかもしれません。
こちらのスライドは発表したものより古いです。Stan2.0やPyStanへの対応も含めて随時アップデートをかけていく予定です。
２．「とある病んだ院生の体内時計（サーカディアンリズム）」（発表者 @berobero11）
発表者コメント：睡眠時刻・起床時刻を扱った循環する変数の統計モデリングの話でした。
とある病んだ院生さんの593日間の記録を可視化した後にRStanの練習に使っています。
StanのHello Worldとも言える内容からvon Mises分布、そして累積時刻に変換して潜在変数のARモデルなどバリエーションに富む内容でした。
各変数の範囲が異なる場合に範囲を指定する方法が参考になると思います。
なお、データやStanコードはこちら（http://heartruptcy.blog.fc2.com/blog-entry-91.html）で公開されています。
berobero11さんの発表でした。とある大学院生の睡眠時間記録のデータに関してのモデリング例を4つほど、紹介して頂きました。
私にとっては珍しい（？）フォン・ミーゼス分布を使った例も出てきますが、後半の仮説検証的なモデリングになるにつれデータ前処理部そのものでのモデリングが増えるため、Stanでのモデリングテクニック⇒データ前処理テクニックということを教えられたような気がします。
最初のフォン・ミーゼス分布を使った時刻モデルから最後の時刻と時間のARを使った累積時刻のモデルへと昇華していく過程はデータ解析の仕事に就いた人が華々しい統計モデリングに挫折し、データ前処理の仕事に終始するようになり、みじめな気分を味わいながらも強く前向きに生きていく壮大な叙事詩、データ解析版のレ・ミゼラブルを見ているような気分にさせられます。
実際にはベイズ推定使おうが他の手法使おうが、前処理でほぼ決まります。このスライドはその泥臭い過程を初めて歩む方にとってのいいロールモデルになるのではないでしょうか。明るく前向きに前処理しましょう。
３．「PyStanで自然言語処理へ向けて」（発表者 @xianze750）
Pystan for nlp from xiangze
発表者コメント：Pystanで自然言語処理を行うにあたっての基礎的事実の紹介、特にLDAとノンパラメトリックなモデルに関する説明。
xianze750さんの発表でした。すみません、私の自然言語処理の知識が足らず途中までしか理解できませんでした。
発表の趣旨としては自然言語処理の階層Dirichletモデルに対してPyStan、JAGSで実装を試みるというものでした。
自然言語処理特有の離散的なデータの扱いはStanとは相性が悪いようです。Stanの今後の機能の改善が期待されます。
４．「『予測にいかす統計モデリングの基本』の売上データの分析をトレースしてみた」（発表者 @berobero11）
発表者コメント：この世の中でありそうな売上の時系列データに対して、
いかにドメイン知識や経験を統計モデリングに組み入れるか、そして各パラメータがいかに簡単にMCMCサンプリングで推定できるかを示した発表になります。
Stanは時系列データに対して相性がよさそうとのことです。詳細（データ含む）はこちら（http://heartruptcy.blog.fc2.com/blog-entry-90.html）で公開されています。
またまた、berobero11さんの発表です。『予測にいかす統計モデリングの基本』の本に出てくる居酒屋の売上予測モデルをStanで実装してみた例です。
データおよびモデル共にberobero11さんの執念を感じるトレースっぷりです。
『予測にいかす統計モデリングの基本』は良書ですが、本にはコードが出てきません。
これを参考に本を読みながら、手元のデータをモデリングをするとマーケティング系の方には勉強になると思います。
５．「尤度フリーMCMCの紹介」（発表者 @kos59215）
発表者コメント：Stan を使えば尤度を勝手に処理してくれますが，複雑な系では尤度を自分で計算してやる必要があります。
今回は尤度が不要な推定手法の ABC を紹介しました。尤度がいらない手法なのに， Stan で使うために尤度の話ばかりになって混乱させてすみません。
kos59215さんの発表でした。ABCという手法により尤度関数の形が分からない状況でもベイズ推定（？）が出来るというお話でした。
ABCはイメージ的には尤度関数の代わりに、恣意的にモデリングしたシュミレーションと重みづけの調整可能なデータの類似度を用いたウェイトの2層をかませることによりパラメータを推定するベイズ推定の手法のようです。
最初の感想としては「これさらに恣意的だけどやって大丈夫なのか？」という感じがしました。
ただ、普段データの自明に持つ関係性がモデルに組み込まれていないことがあり、その関係性だけでも組み込めればと思うことがあります。
そういった状況において、うまくささる状況ではすごく活躍する可能性を秘めた手法なんだと思いました。
こういった手法を生態学という分野の方から教えて頂けるのはありがたいです。
xianze750さんの発表にもありましたが「生態学すごい」です。
補足ですが、Stan2.0では対数尤度の記載の仕方が大きく変わっており、スライドにもある通りlp__からvoid型のincrement_log_probという関数に変わっています。Stan1.Xから使用している方は一度、Stan2.0のマニュアルの23.2章を参照されると良いかと思います。
感想
会場は今回もドリコムさんの100人くらいは入る広い会場をお借りしました。ドリコムさんに感謝！！
参加者は全員で20名程度で、満員でした。主催者側が勉強会の内容に集中できるよう参加者は20人に限定して実施しています。
参加できなかった皆様、本当にすみません。
今回はアカデミック寄りの生態学系の研究者の方とマーケティング系の方が増えたかなという印象を持っています。
裾野が広がって参入障壁が減るといいのですが…
発表内容は今回も素人を寄せつけない、かなりディープな内容でした。
内容はかなりコアに見えますが、例えばberobero11さんの発表内容に関してベイズ推定を実務で応用していこうという初心者の方々には大変勉強となる内容です。
流し読みではなく、しっかり理解できるまでじっくり読むといいと思います。
学生の参加者の方もいらっしゃいましたが、上記の資料を踏み台（？）にして若い方々へ活用の輪が広がっていくと主催者としては感無量です。
次回予定
次回の予定ですが、恐縮ですが主催者側の都合により、スケジュールの決定が2月中旬くらいになりそうです。大変申し訳ありません。
3月か4月あたりの開催を目標にしていますのでよろしければみなさんのご参加、お待ちしております。
発表者募集しています
今回も、ほぼ前回と同じメンバーによる発表であったため、各方面からの発表を期待しています。
次回の発表を希望される方はTwitterで@TeitoNakagawaまで、ご連絡ください。
以上です。関係者の皆様、本当にありがとうございました。
QuantumGISに同心円をらくらくで描く方法
仕事柄、地図を記載するときはスケールバーでスケールを表すことよりも中心点からの同心円で表す事が多いです。
QuantumGISで同心円を描こうとするとき、通常は「ベクタ」→「空間演算ツール」→「バッファ」で
同心円を記述しますが、これだと単一の地点に円が一つしか描けないため、同心円でスケールを表す事ができません。
同心円の個数だけバッファの処理を行ってshpファイルを保存するのも面倒です。これを回避するために以下のような手法を思いついたのでメモ。
用意するのはQuantumGISとPostGIS（理屈的にはSpatial Liteでもいいはず）、最後にQGISプラグインのDB Managerをインストールしてセットアップしておきます。
今回は東京タワーに5km圏のバッファ円を4つプロットしてみます。
DB ManagerでSQLを使って5km圏の同心円のポリンゴンデータを読み込みます。
①「データベース」から「DB Manager」を選択します。
②「PostGIS」からデータベースを選択して、「SQL Window」ボタンを押下します。
③以下のSQLを「SQLクエリ」に記述して、データをロードする。。
select id, ST_Buffer(ST_GeomFromText('POINT(139.745447 35.65861)'), id*0.049456) geom from unnest(ARRAY[1,2,3,4]) id;0.049456は5kmを表す度、POINT(139.745447 35.65861)は東京タワーのことです。
ポイントはunnest句を用いてARRAYから仮のテーブルを作成しているところです。
これをST_Bufferのバッファ半径に掛けることで同心円が複数作成されます。
SQLを記述して「実行」ボタンを押下して、「新レイヤとしてロードする」にチェックを入れ「整数のユニークな値で構成されるカラム」に「id」を、「ジオメトリカラム」に「geom」を指定します。
レイヤ名は適当に名づけて「今ロードする！」ボタンを押すと同心円のポリゴンがロードされます。
④書式設定を変更する。
書式設定を変更して、塗りつぶしをなくして、枠線だけにすると以下のような同心円バッファが記載されます（CRSの設定で楕円になってますが。）。
やったね！
deldirでコンビニ位置データのボロノイ分割を行う。
この記事はRAdventCalendar2013に参加しています。
この記事は4日目になります。
これがはじめてというわけではないですがRのdeldirパッケージを使ってみました。
パッケージdeldirはドロネー（Delaunay）分割とディリクレ(Dirichlet)分割（ボロノイ分割、ティーセン分割）に関連する処理を含むパッケージです。
パッケージdeldirの胆になるのが関数delirです。分割に当たってはLee-Schaterアルゴリズムを使用しているようです。
deldir関数の解説
deldir(x, y, dpl=NULL, rw=NULL, eps=1e-09, sort=TRUE, plotit=FALSE,
digits=6, z=NULL, zdum=NULL, suppressMsge=FALSE, ...)
引数
x,y:分割したい座標の入力点
dpl:リスト
rw:分割を行う境界領域を指定する。デフォルトだと入力点の端に対してプラス10％になる。
sort:TRUEがデフォルト。TRUEにセットするとデータが分割前にソートされ計算が効率的になる。
FALSEだと追加点が最後に付加されるので、デバッグ用に設定するための変数である。
plotit:plotを行うか否かのフラグである。
digists:数値の丸め桁数、デフォルトは6桁である。
z:各入力点に対応する重みの値である。
zdum：各ダミー入力点に対応する重みの値である。
返り値
delsgs:ドロネー三角形エッジリストを表すデータフレームである。エッジの座標点（x1,y1,x2,y2）が入力され、各ドロネー点のIDがそれに続く。
dirsgs:ディリクレ領域（ボロノイ領域）エッジリストを表すデータフレームである。エッジの座標点（x1,y1,x2,y2）が入力され、各ドロネー点のIDがそれに続く。
7番目の点はエッジの両端が境界領域と接しているかどうかを表す。
summary:分割されるデータ点とダミー点の情報を含むデータフレーム
-x,y(,z):点とその重みの位置情報
-n.tri:当該点におけるドロネー三角形の数
-del.area:当該点におけるドロネー三角形の面積の3分の1
-del.wts:del.areaを比率に変化させたもの。
-n.tside:ディリクレタイルの他のディリクレタイルとの接触数
-nbpt:ディリクレタイルの境界領域との接触数
-dir.area:ディリクレタイルの面積
-dir.wts:dir.areaの面積を比率に変化させたもの。
n.data:分割されたデータ点、重複は除く。
n.dum:分割されたデータ点、重複点は除く。
del.area:三角に分割された点の凸包の面積、del.areaの合計と一致。
dir.area:境界領域によって囲まれている面積、dir.areaの合計と一致。
rw:正方形領域の枠。
パッケージdeldirその他の関数
tile.list():deldirを受け取って、ボロノイ分割されたタイルを返す。
plot.deldir():deldirを受け取って、プロットする。
以下では、前回の記事でも用いたコンビニデータから東京のコンビニ6009店*1をボロノイ分割してポリゴンをシェープファイルに出力するコードを記載します。
コード
library(deldir)
library(maptools)
#データは事前に準備している。
#> head(df)
#storename storetype x y
#1 セブンイレブン 世田谷中央病院店 セブンイレブン 139.6511 35.64166
#2 ファミリーマート 池袋二丁目店 ファミリーマート 139.7090 35.73495
#3 みんなのイチバ 足立江北２丁目店 その他 139.7637 35.76921
#ボロノイ分割を行う。rwオプションを指定し外枠を拡張させない。
dd<-deldir(df$x, df$y, rw=c(min(df$x),max(df$x),min(df$y),max(df$y)))
#タイルを取得する。
w<-tile.list(dd)
#ポリゴンベクターの初期化
polys<-vector(mode='list', length=length(w))
for (i in seq(along=polys)) {
#ボロノイ点を取得している
pcrds<-cbind(w[[i]]$x, w[[i]]$y)
#閉路化してポリゴン座標に変更する。
pcrds<-rbind(pcrds, pcrds[1,])
#ポリゴンを作成する。
polys[[i]]<-Polygons(list(Polygon(pcrds)), ID=as.character(i))
}
#スパシャルポリゴンオブジェクトの作成
SP<-SpatialPolygons(polys)
#ボロノイ分割したデータをSpatialDataFrameクラスに変更する。
voronoi<-SpatialPolygonsDataFrame(SP,data=data.frame(x=df$x,y=df$y,storename=df$storename, storetype=df$storetype,
row.names=sapply(slot(SP, 'polygons'),function(x) slot(x, 'ID'))))
#shpファイルにして出力する。
writeSpatialShape(voronoi, fn="vornoi-combini")
結果
以下は結果のシェープファイルをQGISで表示させたものです（赤がセブンイレブン、青がローソン、緑がファミリーマート、オレンジがサークルK・サンクス、その他が紫）。
広域図
詳細図
2013年も一月を切りました。皆様よいクリスマス、年越しを。
*1:データ数は6182だったがジオコードの結果、同一座標となった店舗からはランダムに選出した。
