?

Log in

No account? Create an account

Entries by tag: webtobinq

WebTobinQ, サーバーはGoogle Spread Sheetのcsv公開を使うように変更
karino2
Rのweb版のサブセット、WebTobinQという物を作っているのですが、このデータを独自サービスから、Google Spread Sheetのcsv公開に変更しました。
管理用のUIとかが手抜きで管理が面倒だったのと、職場などでデータ集めるついでに増やしていきたい、と思ったので。
その結果独自関数だったread.serverを廃止してRのread.csvを実装する事にしました(まだまだ手抜きですが)
Exampleを一通り更新しておいたのでこちらを見ていただくと早いかと。
http://www47.atwiki.jp/karino2/pages/36.html

公開しているCSVに関してはこちら
http://www47.atwiki.jp/karino2/pages/46.html

WebTobintって何?という人はこちら:
http://www47.atwiki.jp/karino2/pages/35.html
Tags:

substitute
karino2
WebTobinQでhpfilterをサポートしよう、と思いたち、必要な機能を上から順番に実装している。
現在substituteとdeparseの所で、substituteという物の振る舞いをいろいろ調べていると、これが結構面妖な動きをする。

substituteのドキュメントはこちら。
http://wipaed.wiso.uni-goettingen.de/~holdenb1/R/library/base/html/substitute.html

envの事は忘れれば、通常の呼び出しは普通に渡された引数の見解釈な表現、というか構文木を返す感じに見える。

> substitute(1+2)
1 + 2


でも、関数の中でも使える。

> test <- function(a) { substitute(a); }
> test(1+2+3)
1 + 2 + 3

現状のWebTobinQの実装ではコール時に引数評価を行なっているので、ちょっと呼び出しに細工するだけでは不十分な模様。
ちょっと実験してみよう。存在しない変数を使ってみる。

> test2 <- function(a, b, c) { b+c }
> test2(hoge, 1, 2)
3


参照されなければ、存在しない変数を使っても呼び出しはエラーにならない。
うーむ、なるほど。使うまで評価しちゃいけなかったのか。

> test <- function (a) { testb(a); }
> testb <- function(c) { substitute(c); }
> ika <- 2
> test(ika)
a
あれ?ikaじゃないの?うーむ、面妖な。

どうせ現状だとdeparseとsubstituteはセットでラベル表示するだけなので、envに突っ込む時に元のS式もっとく位でいいかなぁ。
あんまし複雑な式のdeparseを頑張る気は今の段階では無いけれど、自前関数でラベル表示したい、ってのは良くあるので、簡単なシンボルくらいはdeparse出来るようにしておくか。

それにしてもこんなレベルで自分が理解していない部分があるとは思わんかった…
Tags:

所得の中央値を求めて
karino2
以下は前エントリ、「日本人にとっての「経済成長の意味」とは?」を作成する過程の話です。

WebTobinQ
のドッグフードがてら何か手ごろなグラフを描こう、と思ったので、所得の中央値を描いてみよう、と思いました。
なお、まだ足りない機能はたくさんあるので、必要になる都度実装していく方針です。

中央値としては厚生労働省の国民生活基礎調査に近年の結果は載っているのですが、あまり昔のはありません。
とりあえずこの値を計算してみるか、と使えそうな元データを探した所、平成21年の国民生活基礎調査のページに、ここ数年(2004~2008)の相対度数があったので、まずはこれを元に計算してみる事にしました。
所得金額別世帯数相対度数分布という名前でuploadしてあります。

read.server("http://webtobins.appspot.com/t/", "所得金額別世帯数相対度数分布",num=100)

medianの計算にはRの組み込みのmedianは使い出がいまいちなので実装しない事にしました。
if文などを実装するのもかったるいなぁ、と思いどうしようかと悩みながらwebを調べていると、青木さんのページに似た用途の関数が。
いかにもRっぽいコードで気に入ったのでこれを動かす事に。
とりあえずユーザー定義関数とcumsumを実装しました。

さて、この関数と同じような計算で中央値は計算出来そうです。
当該年度をベクトルとして渡し、上記関数と同じような関数で結果を返させるのに必要な機能を考えます。
read.serverの戻りはdata.frameなのですが、階級は今の所titleに文字列としてかえってきます。
Rで言うとattributesのnames属性です。階級の値を算術演算する為にas.numericも実装しました。
data.frameの行をベクターにするには普通はどうするのかな?とちょっと調べてみると、
http://tolstoy.newcastle.edu.au/R/help/04/11/7784.html
には、

as.matrix(df)[1,]

とするのが筋らしい。[1,]の配列アクセスは実装してもいいのですが、この為だけにmatrix型を実装するのも少し遠回りだなぁ、と思っていた所、

as.numeric(df[1,])

でも同じ振る舞いをする模様(ただしas.numericはrowが1行の時以外はNGなのでこのケースでしか使えない)
都合がいいのでas.numericの実装をこのケースもサポートするように変更する事でRとの互換性を維持しつつこのケースを乗り越える事にしました。

以上を組み合わせて、相対度数分布からmedianが計算出来るはずです(物価調整はせず)

df <- read.server("http://webtobins.appspot.com/t/", "所得金額別世帯数相対度数分布")
mymedian <- function(dosuu, levels) {
  pos <- length(dosuu[dosuu<50]);
  dosuu_until_50 <- 50 - dosuu[pos];
 dosuu_haba <- dosuu[pos+1]-dosuu[pos]; 
  levels[pos] + (levels[pos+1]-levels[pos])*dosuu_until_50/dosuu_haba
}
mymedian2 <- function(df, i) {
  dosuu <- as.numeric(df[i,])[-1]
  levels <- as.numeric(attributes(df)[["names"]][-1])
  mymedian(dosuu, levels)
}
med <- 0
for(i in 1:(length(df[["Year"]])) { 
   med[i] <- mymedian2(df, i)
}
med

対話的に確認しながらやっていたら関数が二つになってしまいましたが、本来なら半分くらいのコードに出来そうです。
余談ですが、このfor文

for(i in 1:(length(df[["Year"]])) { 
   med[i] <- mymedian2(df, i)
}
はfor文無しで書けませんかね?少し考えたのですが分かりませんでした。
Rに詳しい人の意見を伺ってみたい所です。

何はともあれこれで2000年から2008年までの中央値は出ました。結果は(2008~2000年)

430.35714285714283 449.1525423728813 454.0816326530612 460.5769230769231 465.4545454545455 477.55102040816325 478.8461538461538 485.7142857142857 500.96153846153845

これを公表されている値と比較してみたい。
厚生労働省の所にある2006年の結果では2006年は451万となっている。比べると、454万なので3万円ほどずれています。
これは450万~500万の間である所までしか度数分布では分からず、その後の計算の仕方が線形補間なせいなのかなぁ、と思いつつもオリジナルのデータの計算方法が分からないのでこれ以上は追求出来ず。

厚生労働省のデータ(を抜き出したと思われるこちらのblog)の値とをグラフにして比較してみると
(注: 厚生労働省のデータは調査年が書いてあるので、どの年次の値なのかを調べるには一年戻す必要がある)

y2 <- c(2007, 2006, 2005, 2004, 2003, 2001, 1999, 1997, 1995)
medS <- c(448,451,458,462,476,485,506,536,550)
plot(df[["Year"]], med, type="l", ylim=c(400,560))
lines(y2, medS)


2000年のデータが違うように見えるのは、2000年のデータは緑の方には無いのでそう見えるだけなのに注意。
完全に一致はしてませんが、だいたい一致しています。このデータでだいたいいいようです。

ただ、相対度数のデータは2000年までしか入手出来ません。元の米国版の研究では1975年からな訳だし、やはり80年代も見てみたい。
どうしたらいいでしょうか?続く。

日本人にとっての「経済成長」の意味とは?
karino2
himaginary氏のブログを見ていたら、「米国人にとっての「経済成長」の意味とは?」という記事がありました。
そこには、実質GDPと所得の中央値の1975年を基準とした比率の計算が。

なかなか興味深かったしWebTobin向きなので日本での結果を求めてみよう、と思い、いろいろ作業してみました。
その過程は後述するとしてまずは結果から。



所得は実質調整済みで、1980年を1としています(元の米国版は1975年基準なので注意)

元になるデータは実質GDPは内閣府の国民経済計算から、所得は国税庁申告所得税標本調査結果から、CPIは総務省統計局からのデータを使っています。
以下も参考までに。
http://webtobins.appspot.com/tables

コードはこちら(注: 以下のコードは http://webtobinq.appspot.com/ にコピペしてEval Allをすると実行出来ます)

median_from_kaikyuusetaisuu <- function(df, i) {
  dosuu <- as.numeric(df[i,])[-1]
  cumdosuu <- cumsum(dosuu)
  total <- sum(dosuu)
  pos <- length(cumdosuu[cumdosuu < total/2])
  levels <- as.numeric(attributes(df)[["names"]][-1])
  levels[pos]+(levels[pos+1]-levels[pos])*(total/2 - cumdosuu[pos])/dosuu[pos+1]
}
df3 <- read.server("http://webtobins.appspot.com/t/", "所得階級別納税者数1989~2009", num=100)
med3 <- 0
for(i in 1:(length(df3[["Year"]])) { 
   med3[i] <- median_from_kaikyuusetaisuu(df3, i)
}
df4 <- read.server("http://webtobins.appspot.com/t/", "所得階級別納税者数1975~1988", num=100)
med4 <- 0
for(i in 1:(length(df4[["Year"]])) { 
   med4[i] <- median_from_kaikyuusetaisuu(df4, i)
}
cpi <- read.server("http://webtobins.appspot.com/t/", "CPI", num=100)
years <- c(df3[["Year"]], df4[["Year"]])
meds <- c(med3, med4)
cpis <- 0
for(i in 1:length(years)){ cpis[i] <- cpi[["All items"]][cpi[["Year"]] == years[i]]}
real <- 100*meds/cpis

realGDP <- read.server("http://webtobins.appspot.com/t/", "実質GDP", num=100)
relativeIncome <- real/real[years==1980]
relativeIncome <- relativeIncome[(years)>=1980 & years < 2010]
gdpYear <- realGDP[["Fiscal Year"]][-1]
realGDPVal <- realGDP[["GDP(expenditure approach)"]][-1]
relativeGDP <- realGDPVal/realGDPVal[gdpYear==1980]
plot(gdpYear, relativeIncome, type="l")
lines(gdpYear, relativeGDP)

所得のmedianは1998年の値が異常値っぽくなっているのですが、理由は良く分かりませんでした。
元データをちらっと眺めた印象も、この年だけ前後と大きく違う。

relativeIncome[gdpYear == 1998]
などとテキストエリアの最後に打って行末でenterするとこの値が見れます。
1997が1.18, 1998が1.4, 1999が1.22です。


米国版ではちょっとしか成長していない、という話でしたが、日本版はむしろ低下している…
少し信じがたいので、計算違いとかかな?と思い、どこかに同じような計算をしたグラフは無いかしら?と思ったのですが、見つかりませんでした。
いろいろ検算等してみたのですが、どこが間違いか良く分からなかったのでとりあえず公開しておきます。

以上の結果が正しいとすると、なかなか驚きの結果ですね。
もともとのコンテキストでは日本や欧州に比べて英語圏は一部の富裕層にだけ富が分配されている、という話だったのですが、日本はもっと酷いという事になる。
1990年くらいまでは日本の方がむしろ良さそうですが、2000年後半からの中央値の低下は自分の印象よりも随分とひどい。

以下、この所得の中央値作成までの話をしていこうと思います。

ユーザー定義関数を実装(returnはまだ)
karino2
所得のmedianをplotしようと思ってデータを集めていた所、度数分布とか納税者数とか世帯数とか時代によってまちまちのデータから作らないといけなさそうです。

この辺の処理はやはりwebtobinqで出来る方がいいだろう、という事でmedianの計算を見ていた所、基本的なコードは

http://aoki2.si.gunma-u.ac.jp/R/median2.html

にありました。
さしあたってこれを実装する為には、sum, cumsum,簡単なユーザー定義関数があれば良さそうです。
returnはしばらく使わなさそうなので、実装は簡単ですが未実装のままにする事にしました。

四則演算が相変わらず変な優先順位なので括弧を補ったり、関数の最初の開き括弧の後の改行に未対応だったりとちまちまとした変更は必要でしたが、

median2 <- function(f, b, w) { cf <- cumsum(f);
n <- sum(f);
pos <- length(cf[cf < n/2]);
b+(w*pos)+(w*((n/2)-cf[pos]))/f[pos+1]
}
f <- c(2, 6, 39, 385, 888, 1729, 2240, 2007, 1233, 641, 201, 74, 14, 5, 1)
median2(f, 55.5, 8)

というコードが動く所まで来ました。
実際のmedianのプロットは次の機会に。

データ集めはじめ
karino2
WebTobinQ自体は大分出来てきたので、実際のデータを集め始めました。
とりあえず興味のあった失業率や人口統計を入れてみました。
テーブルが最初から扱いやすい形ならちょっとした修正でupload出来ますが、複雑なテーブルの場合は手作業ではつらいですね。
WebTobin自体よりも、この元データの加工の方で少し工夫が居る気がしています。

グラフとそのスクリプトをWikiに蓄積しはじめました。
http://www47.atwiki.jp/karino2/pages/36.html

当初はグラフを右クリックして画像で保存出来ると思っていたのですが、実際に保存してみると中がからっぽの画像…
割とhtmlで描いているようで。
http://code.google.com/p/clientsidegchart/issues/detail?id=31#c0

Chartライブラリを選ぶ時にサンプルが保存出来るか、は試したのですが、保存したファイルを開いていませんでした。

保存する方法は欲しいですね。
当面はprint screenで保存しておきます。
Google Chartにはimageに対してrenderingするような仕組みがあるっぽいので、これに乗り換えればいいのかもしれません。
ある程度使ってみてから考えます。スクリプト的には変わらないはずなので(Rサブセットですからね)

WebTobinServerのスキーマを変更
karino2
最初、一つのセルを一つのBigTableのrowとする形で、

class NumericCell(db.Model):
   table = db.StringProperty()
   row = db.IntegerProperty()
   col = db.IntegerProperty()
   val = db.FloatProperty()
という形にしていましたが、bigtableは小さいcellをたくさん返す、という類に弱いようで、パフォーマンスも遅く、fetchの最大数も1000という数字になっていました。
そこで、SQLぽくは無いのですが、csvの一行を一レコードにするように変更しました。
class NumericRow(db.Model):
   table = db.StringProperty()
   row = db.IntegerProperty()
   index = db.FloatProperty()
   vals = db.ListProperty(float)
経済統計の場合、最初の行はどうせYearが入るので、これによる絞込みをサポートすべくindexingするようにしました。
結果としてread.serverの引数のrangeは、fieldNameが無くなります。

少し試した印象だと大分速くなっていて、CPUクオータにもひっかかりにくくなっているとは思います。
実環境でいろいろ実験しているとすぐクオータにひっかかってしまうのでなかなか確認するのも難しいですが。

クオータにひっかかってしまった…
karino2
失業率のデータをアップロードしたり、途中までアップロードされたデータを削除したりすると、すぐにCPUクオータにひっかかってしまいます。
一日すれば復活するのですが、あまりにもすぐにひっかかってしまう…

データはwebからでは無くbulkloaderを使ってアップする方が現実的なんですかね?
不便なので出来たらwebからのオペレーションだけで一通りの作業が出来るようにしたいのですが。

WebTobinQ、サーバーからデータを取得出来るように
karino2
4つのマイルストーンの中で一番の山であった、サーバーからデータを取得する部分が終了しました。
これで当初意図していた形がだいたい見える状態にはなったと思います。
データのサーバーにはとりあえず内閣府の実質GDPのデータが入っています。

一体このWebTobinQが何なのか、というのはこの日記だけだと良く分からないかもしれないので、もし良かったら実際に試してみてください。


GDPと消費のグラフを書く手順

まず以下のurlにいき、
http://webtobinq.appspot.com/

以下のコードを一行づつコピペして行末でEnterを押していってください。
df <- read.server("http://webtobins.appspot.com/t/", "実質GDP", c("Fiscal Year", "GDP(expenditure approach)", "PrivateConsumption"), num=30)
plot(df[["Fiscal Year"]],df[["GDP(expenditure approach)"]],ylim=c(100000,600000))
lines(df[["Fiscal Year"]],df[["PrivateConsumption"]])

その他実際に動いたサンプルは以下にあります。
http://www47.atwiki.jp/karino2/pages/34.html
詳しくはWikiで。
http://www47.atwiki.jp/karino2/pages/31.html

read.serverが動いた!
karino2
Wed, 07 Sep 2011 17:19:51 +0900

サーバー側はまだ引数をほとんど全部無視ですが、大きな流れは完成しました!
サーバーの基本機能を足せばマイルストーンの3は終了です。
そこまでいけば、後は実際に使いながら必要な物を足していくフェーズになります。
いやぁ、この2週間くらい、なかなかこればっか作っていました(^^;