Hatena::ブログ(Diary)

驚異のアニヲタ社会復帰への道

Prima Project

2017-04-20

MikuHatsune2017-04-20

フーリエ変換でみる矩形波

矩形波フーリエ変換の要素数を増やすとどれくらい近似度があがるかという動画を見かけたのでやってみる。

矩形波

f(x)=¥frac{4}{¥pi}¥sum_{n=1}^¥infty¥sin((2n-1)x)

となる。

n が増えると近似度があがり、そのための小さな円が数珠つなぎに増える。

矩形波の平坦な部分はこの数珠がぐるぐる巻にされ、yが入れ替わるところで数珠がダイナミックに伸びるのがわかる。

f:id:MikuHatsune:20170420220633p:image

GIFにするとこんな感じ。

f:id:MikuHatsune:20170420220927g:image

f <- function(n, x) 4/pi*sin((2*n-1)*x) / (2*n-1)
x <- seq(0, 2*pi, length=1000)

p <- 4/pi
y <- f(1, x)
xl <- c(-6.5, 2*pi)
yl <- c(-1, 1)*p*2
xc <- 3.5
Ns <- c(2, 3, 10, 20)
for(i in seq(x)){
  x0 <- x[i]
  png(sprintf("%04d.png", i), 250, 600)
  par(mar=c(1, 5, 1, 1), mfcol=c(length(Ns), 1), cex.lab=1.5)
  for(N in Ns){
    y <- rowSums(mapply(function(z) f(z, x), 1:N))
    plot(x, y, xlim=xl, ylim=yl, type="l", xlab="", ylab=sprintf("n = %d のとき", N), col=gray(0.8))
    lines(x[seq(i)], y[seq(i)], col="red")
    abline(v=x0, lty=3)
    plot.circle(-xc, 0, p)
    x1 <- cos(x0)*p - xc
    y1 <- sin(x0)*p
    points(x1, y1, pch=16)
    for(n in 2:N){
      plot.circle(x1, y1, p/(2*n-1))
      x1 <- x1 + cos((2*n-1)*x0)/(2*n-1)
      y1 <- y1 + f(n, x0)
      points(x1, y1, pch=16)
    }
    segments(x1, y1, x0, lty=3)
    points(x0, y1, pch=16, col="red")
  }
  dev.off()
}


# 実際には lines 関数を呼ぶので,... には lines 関数が許容する引数を書くことができる。
plot.circle <- function(ox, oy, r, start=0, end=360, ...)
{
        plot.ellipse(ox, oy, r, r, 0, start, end, ...)
}

# 中心を (ox, oy) とする,半径 r の円を描き,内部を塗りつぶす。
# 実際には polygon 関数を呼ぶので,... には polygon 関数が許容する引数を書くことができる。
plot.circlef <- function(ox, oy, r, ...)
{
        plot.ellipse(ox, oy, r, r, 0, 0, 360, func=polygon, ...)
}

# 中心を (ox, oy) とする,長径 ra,短径 rb の楕円を描く。
# phi は楕円の傾き(長径が水平線となす角度),start, end には,描き始めと描き終わりの位置を指定できる。
# 長径を基準として,角度(度単位)で指定できる(0, 360 とすると,完全な楕円を描くことになる。
# 90, 270 とすると,左半分の楕円を描くことを指示することになる)。
# 実際には lines 関数を呼ぶので,... には lines 関数が許容する引数を書くことができる。
plot.ellipse <- function(ox, oy, ra, rb, phi=0, start=0, end=360, length=100, func=lines, ...)
{
        theta <- c(seq(radian(start), radian(end), length=length), radian(end))
        if (phi == 0) {
                func(ra*cos(theta)+ox, rb*sin(theta)+oy, ...)
        }
        else {
               x<- ra*cos(theta)
               <- rb*sin(theta)
                phi <- radian(phi)
                cosine <- cos(phi)
                sine <- sin(phi)
                func(cosine*x-sine*y+ox, sine*x+cosine*y+oy, ...)
        }
}

radian <- function(degree) {
        degree/180*pi
}

2017-04-14

MikuHatsune2017-04-14

beamer poster でプレゼンテーション

beamer でスライドプレゼンテーションをしたが、今回はポスター発表をbeamer で作った。

ネットでbeamer poster のsty ファイルを見てみたのだが、華美すぎて鬱陶しい。

個人的には究極グレースケールで終わらせたい。

こういうテンプレ、なんでみんな大好きなんだろうか…

 

tex とsty はgithub でもここでもどちらにでもあげてもいいのだが、うむ。

f:id:MikuHatsune:20170414233124p:image

2017-04-10

MikuHatsune2017-04-10

推計人口 2053年に1億人下回る

という記事があった。

国立社会保障・人口問題研究所日本の将来推計人口(平成29年推計)というものがあったので、表1-1 総人口,年齢3区分(0〜14歳,15〜64歳,65歳以上)別人口及び年齢構造係数:[出生中位(死亡中位)推計]というネ申エクセルデータを図示してみる。

f:id:MikuHatsune:20170410213523p:image

 

これだけだとアレなので出生中位、高位(合計特殊出生率1.65)、低位(1.25)の3パターンをやってみた。中位は明言されていないようだがこの調子だと1.45 の仮定だろう。

実数でみると、いずれの出生率仮定でも大差ないように見える。

f:id:MikuHatsune:20170410220404p:image

 

しかし、総人口比率にすると、「前回推計と比較して人口減少の速度や高齢化の進行度合いは緩和」とはいうが、それでも老年人口は増え続ける。出生高位では2050年後半から頭打ちになって減少傾向になるようだが、合計特殊出生率1.65になるとは到底思えない() 出生低位仮定はかなり悲惨で、老年人口はぐんぐん上がって40% 超えるし、生産年齢人口は減るし、何より若年人口の減少がやばい。

f:id:MikuHatsune:20170410220431p:image

 

上の実数棒グラフを見ると、若年人口、生産年齢人口は目に見えて減っているのに、老年人口の帯幅は50年一定である。つまり、老年になる人口と、死んでここから消えていく人口がトントンになっている、ということである。「平均寿命は、平成27(2015)年男性80.75年、女性86.98年から、平成77年(2065)年に男性84.95年、女性91.35年に伸長(中位仮定)」ということなので、やったね長寿ニッポン万歳()

 

その他、このネ申エクセルデータには、2060年以降の長期予測結果や、出生率をさらに広い幅で振った予測、死亡率を変えた予測などあるので、ネ申エクセルの扱いに慣れたい人や、記述統計を頑張りたい人はがんばろう。

library(gdata)
library(stringr)
# ファイルは表1-1, 1-2, 1-3 を対象にする
files <- list.files(pattern="pp29gt\\d{4}")
dat <- mapply(read.xls, files, skip=2, SIMPLIFY=FALSE)
dat <- lapply(dat, head, -1)

# 年度を処理
year <- as.numeric(str_extract(dat[[1]][,2], "\\d+"))

# 生数値を処理
dat <- mapply(function(z) apply(apply(z[,4:6], 2, gsub, pattern=",", replace=""), 2, as.numeric)/100, dat, SIMPLIFY=FALSE)
for(i in seq(dat)) colnames(dat[[i]]) <- c("0-14", "15-64", "65-")
names(dat) <- sprintf("出生%s位", c("中", "高", "低"))

# 実数値棒グラフ
cols <- c("orange", "yellow", "grey")
par(mfrow=c(1, 3))
for(i in seq(dat)){
  barplot(t(dat[[i]]), horiz=TRUE, col=cols, main=names(dat)[i])
}


# 総人口に対する割合
yl <- c(0, 0.6)
par(mfrow=c(1, 3))
for(i in seq(dat)){
  y <- dat[[i]]/rowSums(dat[[i]])
  matplot(year, y, type="l", col=cols, lwd=5, lty=1, ylim=yl, ylab="総人口割合", main=names(dat)[i])
  for(j in seq(dat)) lines(year, y[,j])
  legend("right", legend=colnames(dat[[1]]), col=cols, lwd=3, lty=1)
}

2017-04-05

MikuHatsune2017-04-05

MeshLab を使って曲率を取得するのをバッチ処理して自動化する

3次元物体の曲率を以前やっていたのだが、自前のスクリプトに重大なミスが発覚したので、プログラミング力の高い人のプログラムから借用することを考える。

離散的な3D オブジェクトから、平均曲率とガウス曲率を計算したいが、意外とない。

と思ったらMeshLab にあるようである。

.obj を読みこんで、Filters > Color creation and processing > Discrete curvature と選択すると、mean かgaussian か選べと出るので、どちらか選んでapply すると、こんな感じになる。

色は Filter > Quality measure and computations > Quality mapper applier から、カラースケールを適当にいじる。

Render > Show vert quality histogram を選ぶと、今回は頂点について曲率を計算しているので、頂点の値に応じたヒストグラムが出てくる。これはGUI では背景が濃い青だったので白抜きの文字が本当は見えている。

f:id:MikuHatsune:20170405201643p:image

 

ここで、この値情報を保持したままファイルに書き出してその値が欲しいのだが、値の情報は .ply 形式しか保持してくれないようである。なので、ファイルの保存は .ply にする。

ここで、.ply にはbinary 形式とascii 形式があるので、export のときにadditional parameters のbinary encoding のチェックをはずさないと、その後のデータ処理で困る。

 

というのがGUI での一連の操作である。ここで、この処理をしないといけないファイルが複数あるとき、for loop に入れたい。

この時に使えるのがmeshlabserver コマンドである。

まず、このGUI での操作を、Filters > Show current filter script に行って、マクロ化したスクリプト .mlx ファイルを手にいれる。ここでは、Filter タブでの操作(のみ、たぶん)が記録される。

こうして、GUI 操作がマクロ化されたので、

meshlabserver -i input.obj -o output.ply -s procedure.mlx -om vq > /dev/null 2>&1

とする。入力はなんでもよく、出力に .ply を指定すれば勝手に .ply 形式で出力してくれる。-s でやりたい操作.mlx を指定して、今回は頂点の曲率の値が欲しいので、-vq で vertices quality 値を .ply ファイルに書き出しておく。

ここで、meshlabserver の標準出力がどうしても邪魔なので、/dev/null 2>&1 であの世へ送る。

 

この場合に問題になるのが、-o output.ply がどうしてもbinary 形式の .ply ファイルしか書き出さないようで、けっこうなmeshlab ユーザーが困っているようである。というわけで、binary 形式で出力した .ply をascii 形式の .ply に変換したいのだが、どうもいいのがないようなのでR で頑張る。

R のRvcg パッケージに .ply ファイルをbinary 形式でも読める関数があるので、これで .ply ファイルを読み込み、目的のquality 値だけ取り出して出力できるようにバッチ処理にしてみる。

# quality.R
library(Rvcg)
args <- commandArgs(TRUE)
ply <- vcgPlyRead(args)
cat(sprintf("%f", ply$quality))
cat("\n")

これをバッチ処理すると、コンソールにquality 値がずらずらと出力されるので、適当にファイルに保存する。

Rscript --slave --vanilla quality.R $ply.ply > `printf %04d $ply`.tmp

R でバッチ処理したものを出力すると、

[1] 1

というのが残るので、cat して改行コードを付け足しておく。

 

とすればなんとか取得できる。

f:id:MikuHatsune:20170405204033p:image

2017-04-01

ゼロから作るDeep Learning Pythonで学ぶディープラーニングの理論と実装

読んだ。

COI:自費で買った。著者は知り合いでない。

 

本当にゼロから始めるDeep learning. Python のみで完結する。

パーセプトロンニューラルネットでどのように線形分類、非線形分類ができるかを解説し、最初はパラメータを自分で与えないといけなかったが、伝播法(順・逆)により自動で最適なパラメータ設定を見つける過程を実装している。

この本のいいところは、最初に「対象とする読者」のほかに「対象としない読者」を定義しており、CNN を使ったMNIST の分類のタスクをするが、それ以外はお断り(というか詳細には解説しないので他をあたれ)と書いている。とはいっても、分類問題だけでなく、回帰問題を解くときにどうしたらいいか、は書いてある。本書は画像解析だけなので、畳み込みの原理などを解説している。

 

個人的には、数式読めないマンなので、勾配法により最適パラメータを見つけ出すときのステップで、

を読んでいたけど、逆伝播はあ〜微分がi 層からi-1 層について書けるんだな〜くらいにしか理解していなかったが、Python でゼロから作る、というだけあって、逆伝播の仕組みがよくわかった。

逆に、よろしくなかったのが、逆伝播の説明の時に「掛け算ならひっくり返して…」みたいな書き方をしていて、そこでむしろわかりにくかった、気がする。

しかし、単著なので、全章に渡って一貫した感じの書き方なので、わかりやすかった。

 

ゼロから作るとはいえ、本書だけではスクリプトはカバーできず、後半戦はほとんどサイトの○○のスクリプトを参照してください、という感じである。対象者にPython が初めての人、とはあるが、決してプログラミングが初めての人向けではないことは書いておきたい。というのも、最低限の線形代数知識は必要だし、Python の扱いも、スクリプトを書いて関数を定義しておいて、それを呼び出して組み合わせながら、という感じなので、まあ自分がクラスとかそういう概念がいまだにあやふやプログラミング初心者というのもあるのだが、ゼロから作る(ゼロとは言っていない 状態に陥る可能性はあるのはある。せめて困った時に超初歩を聞ける玄人が近くにいれば、本当にプログラミングもDeep learning も初心者です、という人でも出来ないことはないかも知れない(適当

 

安いので買い。