Hatena::ブログ(Diary)

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

Prima Project

2014-12-26

MikuHatsune2014-12-26

第26羽 ココア「チノちゃん、スリーサイズ教えて」チノ「嫌ですよ」

この記事はごちうさ住民 Advent Calendar 2014の第26日目の記事です(勝手に参加。

 

ひと目で尋常でない解析ネタだと見抜いたよ

f:id:MikuHatsune:20141225225419j:image:small:left「チノちゃん、スリーサイズ教えて」

 

 

f:id:MikuHatsune:20141225225504j:image:small:left「嫌ですよ」

 

 

f:id:MikuHatsune:20141225225419j:image:small:left(´・ω:;.:...

 

 

 

データを愛した少女と解析に愛された少女

f:id:MikuHatsune:20141225225419j:image:small:left「というわけでお馴染みのスリーサイズ解析をするよっ。

まずは公式設定から年齢と身長だけ入手できたから

アイマスデータセットと機械学習を用いて年齢と身長から、体重、スリーサイズを推定するよっ。

 

主成分分析をする

ブラのサイズ推定をする

SPADE法で分化解析をする

Trajectory detection で成長度推定をする

viSNEによる次元削減とクラスタリングをする

 

f:id:MikuHatsune:20141225225419j:image:small:leftという感じでやろう」

 

 

 

初めてデータを扱った日の事憶えてる? 他人のデータでスリーサイズ推定しようとしたわよね

f:id:MikuHatsune:20141225225505j:image:small:leftアイマスデータセットは年齢、身長、体重、スリーサイズの6次元のデータで

アイマスのみならずスクフェス、GF(仮)、WUG、ハナヤマタの登場人物を含む369人のデータセットだ」

 

f:id:MikuHatsune:20141225225506j:image:small:left「推定は面倒なのでrandomForest パッケージよりrandomforestを使用するわ〜。

交差検証や推定精度は面倒くさいのでやらないけど、適当な推定精度はあるわ〜」

 

f:id:MikuHatsune:20141225225505j:image:small:left「私の胸が79.1…だと…!?」

 

 

 

ラッキーアイテムはおっぱいと罪と罰

f:id:MikuHatsune:20141225225419j:image:small:leftスリーサイズが推定できたら、次はブラのサイズ推定をしようねっ。

身長、おっぱい、ウエストからなるおっぱい関数を用いてブラのサイズを推定するよっ」

 

name	age	height	weight	B	W	H	cup
ココア	15	154	42.8	78.8	56	80	B
チノ	13	144	37	72.3	53.1	75.1	AA
リゼ	16	160	45.2	79.1	57	79.2	B
チヤ	15	157	43.9	81.9	56.1	82	D
シャロ	15	151	40.7	78.6	55.9	80	B
マヤ	13	140	35.4	72.3	52.9	75	AA
メグ	13	145	37.9	73	53.1	75.8	A
青山ブルーマウンテン	25	163	46.2	84.2	57.8	85.1	D

 

f:id:MikuHatsune:20141225225505j:image:small:left「私のブラのサイズがB…だと…!?」

 

 

f:id:MikuHatsune:20141225225507j:image:small:left「リゼ先輩と同じサイズで嬉しいッ…」

 

 

f:id:MikuHatsune:20141225225419j:image:small:left「チノちゃんはまだまだ大きくなるよっ」

 

 

f:id:MikuHatsune:20141225225504j:image:small:left「余計なお世話です」

 

 

f:id:MikuHatsune:20141225225507j:image:small:left「このおっぱいがBカップなわけない…

ちなみに長内転筋腱がきちんと作画に描きこまれているのは高評価…」

 

f:id:MikuHatsune:20141225232451p:image

第1羽

 

f:id:MikuHatsune:20141225225419j:image:small:left「全然関係ないけど

ここ、わんこが濡れなくてよかったっ!!だからねっ!!」

 

f:id:MikuHatsune:20141225232452p:image

第4羽

 

ココアと悪意なき殺意

f:id:MikuHatsune:20141225225419j:image:small:left「チノちゃんチノちゃん、今は小さくてもまだまだ大きくなるはずだよっ」

 

 

f:id:MikuHatsune:20141225225504j:image:small:left「うるさいですね」

 

 

f:id:MikuHatsune:20141225225419j:image:small:left「チノちゃんがどれくらい大きくなるか解析してみようよ」

 

 

f:id:MikuHatsune:20141225225506j:image:small:left「とりあえずPCAね〜」

 

 

f:id:MikuHatsune:20141225232453p:image

f:id:MikuHatsune:20141225225505j:image:small:left機械学習で推定したデータを用いているから、結局密集してしまう感じだな」

 

 

 

解析をする解析

f:id:MikuHatsune:20141225225419j:image:small:left「SPADEというBioinformaticsの最新手法があるんだよ。とりあえずやっておいたよ。

スリーサイズデータを多次元データとみなして、各マーカーの発現量に応じたFCM解析と同様なものと見做すと、

幹細胞の分化がそのままチノちゃんの成長分化とそっくりそのまま考えられるんだよっ!!」

f:id:MikuHatsune:20141225232454p:image

f:id:MikuHatsune:20141225225505j:image:small:left「チマメ隊は見事にロリ集団だな」

 

 

f:id:MikuHatsune:20141225225504j:image:small:left「納得いきません」

 

 

 

Call Me Sister.

f:id:MikuHatsune:20141225225419j:image:small:left「チノちゃんチノちゃん、WanderlustというBioinformaticsの最新手法もやってみたよっ。

 

 

f:id:MikuHatsune:20141225225505j:image:small:left「多次元データ内のあるサンプルを分化開始の細胞と見做すと、そこからの分化度を定量化する手法だな。

SPADEと同じように考えると、一番ロリな女の子からどれくらい成長しているかが分かるな」

 

f:id:MikuHatsune:20141225232455p:image

f:id:MikuHatsune:20141225225419j:image:small:left「チノちゃんチノちゃん、私、チノちゃんのお姉ちゃんになったよっ」

 

 

f:id:MikuHatsune:20141225225504j:image:small:left「お姉さんはいりません」

 

 

f:id:MikuHatsune:20141225225507j:image:small:left「(リゼ先輩も私よりお姉ちゃんに…ゴクリンコ)」

 

 

 

プールに濡れて 雨に濡れて涙に濡れて

f:id:MikuHatsune:20141225225506j:image:small:left「この回はヘタなエロゲ原作アニメよりエロいわ〜」

 

 

f:id:MikuHatsune:20141225225505j:image:small:left「viSNE は生のパラメータと次元削減後のパラメータが似る確率をKL情報量を用いて推定する手法だな。

ぶっちゃけよくわかってないぞ」

 

f:id:MikuHatsune:20141226102111p:image

 

f:id:MikuHatsune:20141225225507j:image:small:left「(リゼ先輩と離れちゃった…)」

 

 

f:id:MikuHatsune:20141225225505j:image:small:left「チマメ隊はまた一緒だな」

 

 

f:id:MikuHatsune:20141225225504j:image:small:left「もうチマメ隊でいいです…」

 

 

 

青山スランプマウンテン

f:id:MikuHatsune:20141225225507j:image:small:left「(髪下ろした先輩のワンピース姿が可愛すぎる…(青ブルマン関係ない)」

 

 

 

対お姉ちゃん用決戦部隊、通称チマメ隊

f:id:MikuHatsune:20141225231128j:image:small:left「なんでメグだけAカップなんだよ!!」

 

 

f:id:MikuHatsune:20141225231129j:image:small:left「ワ,ワカンナイヨー」

 

 

f:id:MikuHatsune:20141225225504j:image:small:left「メグだけわずかに身長が私より高いので、体重もおっぱいも大きくなったようですね。私もがんばれば…」

 

 

f:id:MikuHatsune:20141225225419j:image:small:left「チノちゃんは今のままでいいよ〜」

 

 

f:id:MikuHatsune:20141225225505j:image:small:left「(なんで私がBカップなんだ…)」

 

 

 

社畜は白い外套を纏いウサギを駆りて聖夜の業務を征く

この解析は社畜の合間をぬって24,25日に行われました。

25日は所属部局でクリスマスコンサートと称して合唱をさせられました。

 

君のためなら精進する

f:id:MikuHatsune:20141225225419j:image:small:left「うーん…、やっぱりうまくいかないなあ」

 

 

f:id:MikuHatsune:20141225225504j:image:small:left「じゃあ…、実際に私で確かめてもいいですよ」

 

 

二人は夜の街へと消えていった…

 

結論:SSが雑!!




推定結果

f:id:MikuHatsune:20141225225419j:image:small:leftココア

15歳/154cm/42.8kg/78.8cm/56.0cm/80.0cm/B

 

f:id:MikuHatsune:20141225225504j:image:small:leftチノ

13歳/144cm/37.0kg/72.3cm/53.1cm/75.1cm/AA

 

f:id:MikuHatsune:20141225225505j:image:small:leftリゼ

16歳/160cm/45.2kg/79.1cm/57.0cm/79.2cm/B

 

f:id:MikuHatsune:20141225225506j:image:small:leftチヤ

15歳/157cm/43.9kg/81.9cm/56.1cm/82.0cm/D

 

f:id:MikuHatsune:20141225225507j:image:small:leftシャロ

15歳/151cm/40.7kg/78.6cm/55.9cm/80.0cm/B

 

画像は公式HPから。

 

スクリプト

# スリーサイズ推定
g <- data.frame(age=c(15,13,16,15,15,13,13,25), height=c(154,144,160,157,151,140,145,163))
rownames(g) <- c("ココア", "チノ", "リゼ", "チヤ", "シャロ", "マヤ", "メグ", "青山ブルーマウンテン")
dat <- read.delim("girl.txt")
rownames(dat) <- dat$name
dat <- dat[, -c(1, 8)]

library(randomForest)
k <- c("weight", "B", "W", "H")
# 年齢と身長だけから体重とBWHを推定する
rf <- mapply(function(x) randomForest(eval(parse(text=k[x])) ~ age + height, data=dat), seq(k), SIMPLIFY=FALSE)
p1 <- sapply(rf, predict, g)
colnames(p1) <- k
g1 <- cbind(g, p1)

# 推定したデータからもう一回推定する
rf <- mapply(function(x) randomForest(eval(parse(text=k[x])) ~ . - eval(parse(text=k[x])), data=dat), seq(k), SIMPLIFY=FALSE)
p2 <- sapply(rf, predict, g1)
colnames(p2) <- k
g2 <- cbind(g, p2) # 推定したデータ

# PCA
# 公式HPから Twitter画像を取っておく
pngs <- list.files(pattern="jpg")
library(png)
library(jpeg)
pics <- mapply(function(x) readJPEG(x, native=TRUE), pngs)
ra <- 1 #原点に近いところが潰れるので拡大したかったけど、等倍でやった。
xy0 <- sapply(pics, dim)[1:2, ] #pixel
xy0[] <- min(xy0)
xy0[2,] <- 200 # なんか横が潰れるのでテコ入れ
rownames(xy0) <- c("height", "width")
s0 <- 0.008 #拡大縮小率
data1 <- rbind(dat, g2)
pca_score <- scale(data1) %*% eigen(cor(data1))$vectors *sqrt(nrow(data1)/(nrow(data1) - 1))
par(mar=c(4.5, 5, 3, 2), cex.lab=2)
plot(pca_score[, 1:2], type="n", xlab="", ylab="← グラマー	年相応スタイル スレンダー →")
mtext("← 大人\tスタイル\tロリ →", 1, line=3, cex=2)
title("アニメキャラ分析", cex.main=2)
abline(h=0, v=0, lty=3, col=grey(0.5), lwd=2)
text(pca_score[, 1:2], rownames(pca_score))
lay0 <- pca_score[rownames(g2), 1:2]
for(i in seq(pics)){
	xleft=lay0[i, 1]*ra - xy0[2, i]/2*s0
	ybottom=lay0[i, 2]*ra - xy0[1, i]/2*s0
	xright=lay0[i, 1]*ra + xy0[2, i]/2*s0
	ytop=lay0[i, 2]*ra + xy0[1, i]/2*s0
	rasterImage(image=pics[[i]], xleft=xleft, ybottom=ybottom, xright=xright, ytop=ytop, xpd=TRUE)
}

# SPADE
# 前回と比べていろいろ仕様が変わってるんだけど…
library(spade)
library(flowViz)
dat <- read.csv("animegirls.csv") # 推定したごちうさメンバーを含むデータセット
dat1 <- subset(dat, select=-c(name, type))
fcs <- new("flowFrame", as.matrix(data1))
fcs <- new("flowFrame", as.matrix(round(data1*100))) # 強制的に整数化
write.FCS(fcs, "animedata.fcs", what="integer")

downsample_file_path <- paste(output_dir, "animedata.fcs.downsample.fcs", sep="/")
cells_file_path <- paste(output_dir, "/", "clusters.fcs", sep="")
clust_file_path <- paste(output_dir, "/", "clusters.table", sep="")
graph_file_path <- paste(output_dir, "/", "mst.gml", sep="")

set.seed(2)
data_file_path = "animedata.fcs"
output_dir <- "fcs"
SPADE.driver(data_file_path, out_dir=output_dir, k=30, clustering_samples = 5000)
density_file_path <- paste(output_dir, "animedata.fcs.density.fcs", sep="/")
SPADE.addDensityToFCS(data_file_path, density_file_path)
SPADE.FCSToTree(downsample_file_path, cells_file_path, graph_file_path, clust_file_path, k=30)

upsample_file_path <- "upsamle.fcs"
SPADE.addClusterToFCS(density_file_path, upsample_file_path, cells_file_path)
up <- read.FCS("upsamle.fcs")
write.table(up@exprs[, "cluster"], "clusterID.txt", row.names=FALSE, col.names=FALSE)
cl <- unlist(read.csv("clusterID.txt", header=FALSE)) # クラスター番号
table(cl) # クラスターに所属する女の子の数
tab <- paste(getwd(), "/", output_dir, "/tables/bySample/animedata.fcs.density.fcs.cluster.fcs.anno.Rsave_table.csv", sep="")
cvs <- read.csv(tab) # SPADE を実行した後のいろいろな統計量

mst_graph <- igraph:::read.graph(paste(output_dir,"mst.gml",sep=.Platform$file.sep),format="gml")
clust <- read.table(paste(output_dir, "/clusters.table", sep=""), sep=" ", header=TRUE)
lay0 <- read.table(paste(output_dir,"layout.table",sep=.Platform$file.sep))

pngs <- list.files(pattern="jpg")
library(png)
library(jpeg)
pics <- mapply(function(x) readJPEG(x, native=TRUE), pngs)
ra <- 1 #原点に近いところが潰れるので拡大したかったけど、等倍でやった。
xy0 <- sapply(pics, dim)[1:2, ] #pixel
xy0[] <- min(xy0)
xy0[2,] <- 200 # なんか横が潰れるのでテコ入れ

rownames(xy0) <- c("height", "width")
s0 <- 0.0025 #拡大縮小率
b <- list(c(2,6,7),1,c(3,5),c(4),c(),c(8)) # 画像を貼る順番

par(mfrow=c(2, 3), mar=c(0, 0, 3, 0))
for(m in seq(ncol(dat))){
	l <- colnames(dat)[m]
	cut0 <- seq(min(dat[, l]), max(dat[, l]), length=99) # ノードの色付け
	g0 <- mst_graph
	V(g0)$label <- NA
	V(g0)$size <- log(cvs$count, 1.15)
	V(g0)$frame.color <- "black"
	f <- tapply(dat[, l], cl, median)
	V(g0)$color <- bluered(length(f))[rank(f, ties.method="random")]
	plot(g0, layout=as.matrix(lay0))
	title(colnames(clust)[m], cex.main=3)
	for(i in b[[m]]){
		xleft=xy[i, 1]*ra - xy0[2, i]/2*s0
		ybottom=xy[i, 2]*ra - xy0[1, i]/2*s0
		xright=xy[i, 1]*ra + xy0[2, i]/2*s0
		ytop=xy[i, 2]*ra + xy0[1, i]/2*s0
		rasterImage(image=pics[[i]], xleft=xleft, ybottom=ybottom, xright=xright, ytop=ytop, xpd=TRUE)
	}
}

# Wanderlust, Trajectory detection
# s の設定
loli <- which(rownames(data1) == "横山千佳")
res <- trajectory_detection(data1, s=loli, k=35, l=15, nl=20, ng=30)

# おっぱい関数をやったあとで
a <- as.data.frame(t(mapply(function(i) unlist(oppai(dat$height[i], dat$B[i], dat$W[i] ,correct=TRUE)[2:4]), seq(nrow(dat)))))
rownames(a) <- rownames(dat)

ll <- tail(dat, 8)
ub <- as.numeric(as.vector(a$u_bust))
lv_idx <- replace(seq(12), 1:4, 4:1)
cup_size <- c("AA",LETTERS[1:4])
cols <- rainbow(length(cup_size))[match(ll$cup, cup_size)]
idx <- match(rownames(ll), rownames(dat)[order(dat$score)])

pngs <- list.files(pattern="jpg")
library(png)
library(jpeg)
pics <- mapply(function(x) readJPEG(x, native=TRUE), pngs)
ra <- 1 #原点に近いところが潰れるので拡大したかったけど、等倍でやった。
xy0 <- sapply(pics, dim)[1:2, ] #pixel
xy0[] <- min(xy0)
xy0[2,] <- 40000 # なんか横が潰れるのでテコ入れ
rownames(xy0) <- c("height", "width")
s0 <- 0.001 #拡大縮小率

par(mar=c(4.5, 4.5, 2, 2), cex.lab=1.6, cex.axis=1.6)
plot(sort(dat$score), pch=16, xlab="アニメキャラ", ylab="Trajectory score")
points(idx, ll$score, col=cols, pch=16, cex=2)
legend("topleft", legend=cup_size, col=rainbow(length(cup_size)), pch=16, bty="n", cex=1.5, title="ブラサイズ")
#lay0 <- matrix(unlist(locator(8)), nc=2)
for(i in seq(pics)){
	arrows(lay0[i,1], lay0[i,2], idx[i], ll$score[i], lwd=3, length=0.15, col=cols[i])
	xleft=lay0[i, 1]*ra - xy0[2, i]/2*s0
	ybottom=lay0[i, 2]*ra - xy0[1, i]/2*s0
	xright=lay0[i, 1]*ra + xy0[2, i]/2*s0
	ytop=lay0[i, 2]*ra + xy0[1, i]/2*s0
	rasterImage(image=pics[[i]], xleft=xleft, ybottom=ybottom, xright=xright, ytop=ytop, xpd=TRUE)
}

# viSNE
library(tsne)
# 論文で言っているパラメータに合わせる。
# 他のパラメータは、t-SNEの元の論文のデフォルトパラメータにしてある。
tsne_sub2 <- tsne(sub2, max_iter = 500, perplexity=30, whiten=FALSE)

library(gplots)
cup_size <- c("AAA未満","AAA","AA",LETTERS[1:7],"H以上")
cols <- colorpanel(length(cup_size), "blue", grey(0.9), "red")

pngs <- list.files(pattern="jpg")
library(png)
library(jpeg)
pics <- mapply(function(x) readJPEG(x, native=TRUE), pngs)
ra <- 1 #原点に近いところが潰れるので拡大したかったけど、等倍でやった。
xy0 <- sapply(pics, dim)[1:2, ] #pixel
xy0[] <- min(xy0)
xy0[2,] <- 160 # なんか横が潰れるのでテコ入れ
rownames(xy0) <- c("height", "width")
s0 <- 0.05 #拡大縮小率
#idx <- match(g$name, rownames(dat))
lay0 <- tsne_sub2[idx,]
par(cex.axis=1.6)
plot(tsne_sub2, pch=16, xlab="", ylab="", type="n")
abline(h=0, v=0, lty=3, lwd=2, col=grey(0.2))
points(tsne_sub2, pch=16, col=cols)
legend("topright", legend=cup_size, col=cols, pch=16, bty="n")
for(i in seq(pics)){
	xleft=lay0[i, 1]*ra - xy0[2, i]/2*s0
	ybottom=lay0[i, 2]*ra - xy0[1, i]/2*s0
	xright=lay0[i, 1]*ra + xy0[2, i]/2*s0
	ytop=lay0[i, 2]*ra + xy0[1, i]/2*s0
	rasterImage(image=pics[[i]], xleft=xleft, ybottom=ybottom, xright=xright, ytop=ytop, xpd=TRUE)
}

2014-03-22

MikuHatsune2014-03-22

Rで始めた医学・統計学・Bioinformatics

という本を書いた。

Rで始めた医学・統計学・Bioinformatics

Rで始めた医学・統計学・Bioinformatics

アマゾンKDPから買うと文字化けするため、こちらでは買わず代わりにとらのあなMelonbooks(審査中)で買ってください。

platexで書き始めたはいいが、電子書籍化するためにepubもしくはhtml化しようと思ったけど、epubpngしか受け付けずしかも数式、表はすべてレイアウトグダグダ、htmlはutf8にしてもKDPで文字化けがどうしても直らないので諦めた。

 

内容としてはこんな感じ。参考文献について自分周辺のものについてだけリンク貼っておいた。

 

グリコシミュレーション数理モデル

Rを始めたときにやったのがグリコ。そこから推移行列などで数理モデル化までやってみた初心者だったころの私。

プログラミングセミナー グリコ

京大入試数学2014

 

p値

統計のラボだったので、長らく統計について勉強させていただいた。プログラミングによるシミュレーションを併用すると理解が深まると思った。

ryamadaの遺伝学・遺伝統計学メモ マルチプルテスティングとFDR オミックス統計学入門2014

p値

サンプルサイズ

カプランマイヤー曲線のサンプル数

12/12 MIKUセミナー

数学いらずの医科統計学 第2版

数学いらずの医科統計学 第2版

 

嫁たちのスリーサイズ解析と予測

機械学習を勉強し始めてからの応用例。そのあとはバイオインフォマティクス手法なども用いてみたり。

ラブライブ!というアニメを観ていて

ガールフレンド(仮)のキャラ分析

ボディイメージ

(^q^)「くおえうえーーーるえうおおおwwwwwwwwwwwwwwwwwwwww」

機械学習を用いて双葉杏のスリーサイズ推定問題を解く

ミス・モノクロームのプロフィール推定

ガールフレンド(仮)のキャラ分析

Wake Up, Girls!のキャラ分析

Rでフローサイトメトリー(FCM/FACS)

ryamadaの遺伝学・遺伝統計学メモ cytoSPADE

SPADEを使いやすくする

二次元キャラの分化(成長)を多次元データ解析法で真面目にやってみる

viSNEによる次元削減とプロット

アイマス, ラブライブ, WUG, GFのキャラ分析

FRaC Feature Regression and Classification

外れ値となるアイマスメンバーを探す

デンドログラムをヒートマップのまわりに描き足す

 

ぼくのかんがえたさいきょうのせいゆうキャスティング

日本声優統計学会に声をかけていただいてから、解析がさらにおもしろくなった。ベタな線形回帰を使って、声優出演情報から円盤売上予測をしようという話。本書ではデータ取りの苦労話は割愛されているのでブログを読んで欲しい。次回もし参加することになっても、忙しすぎて原稿を書く時間がなければこの章がコピペで出されると思う。

ぼくのかんがえた さいきょうの せいゆう キャスティング

データ

重回帰法の落とし穴

分散拡大係数 VIF

機械学習の精度評価

交絡

 

これならわかる Excelで楽に学ぶ多変量解析

これならわかる Excelで楽に学ぶ多変量解析

 

百合ネットワーク解析

一時期ネットワーク解析にハマって、声優統計にも書いたネタ。

pixivのタグ頻度から考えるラブライブのカップリング

GRAph ALigner Algorithm (GRAAL)

GRAph ALigner Algorithm (GRAAL) を真面目にやる

ラブライブ!各話の百合ネットワークの動的変化

ネットワークの動的進化

R+igraphではじめる生物ネットワーク解析

 

アニメキャラで考える遺伝学

統計遺伝学的な話。ベイズ的思考も取り入れて事後確率の変動をどう考えるかという話もいれた。

海の人間と陸の人間の遺伝学

伊藤誠で作成しようと思ったら沢越止のほうが鬼畜だと気づいた

ベイズ的思考な遺伝相談

遺伝統計学の基礎―Rによる遺伝因子解析・遺伝子機能解析―

遺伝統計学の基礎―Rによる遺伝因子解析・遺伝子機能解析―

 

 

医師国家試験で考える勉強態度

マークシート試験で適当に解答しても合格基準を超えるのではないかという疑問があって、どれくらい勉強したらいいかを考えた。得点予想など時系列解析も含む。

試験の合格基準は6割です

試験に合格するためにぎりぎりを目指すのもいいけどどれくらいがんばればいいかわからない人に

試験問題が多いと大変ですね()

挑戦回数による合格率の変化

項目応答理論

 

線形計画法ストーキング

シンプレックス法と並列計算によるシミュレーションで、班分けを考えた話。並列計算はこのときに勉強したが、まさかこんなことで使うと思わなかったし、このときの経験が今後ものすごい活きてくるとはまったく思わなかった。

シンプレックス法を用いたコース割当シミュレーション

線形計画法を用いて当直の最適な割り当てを考える

ryamadaの遺伝学・遺伝統計学メモ 2次計画問題

安定マッチング問題

就職活動希望先の地域差

沖縄だけ絶妙に例外的に描く

 

ド素人が始めるテキストマイニング

テキストマイニングもやってみると面白かった。周りの人にはこれをやっている人は少ないけれども、この応用範囲は広くて、これ自体に興味を持っている人が多かったので勉強になったし、話のネタにもなった。

MeCabの辞書をはてなキーワードで充実させるのにものすごい苦労した話

声優統計第二号 トピックモデルを用いたニコニコ動画コメントデータの声優トピック流行推移解析

初音ミクの流行解析をDTMで

トピックモデルを使ってラブライブの歌を解析する

ラブライブ スクフェスの楽曲属性をCTMで予測する

47の心得シリーズをトピックモデルで分類する

 

声優の声に魅せられて

声優統計の音声解析に興味を持って、Rでもできるらしいと知ってやってみた。やってみたはいいがかなり難しく、これ以上の勉強が進んていないのは残念。

声優統計第三号 複数の声優によるセリフの音響的類似性の考察:不愉快です

 

RでGIS

空間的なものの処理とか統計とかそういう知識はほとんどないのだけれども、

祇園祭の山鉾の最適巡回経路

RGoogleMapでGoogleと連携してiPhone5sの繋がりやすさをプロットする

 

Rでお絵描き

Rの真骨頂は柔軟なプロット機能だと思っているので、Rでアイコンを作成することなぞ造作も無い。

みくみくにしてやんよ

みくみくにしてあげる

初音ミク関数

おっぱい関数

KABIRAの日記 3dの拡散

KABIRAの日記 3dでの拡散

細胞カウント

 

あとがき

バイオインフォマティクスといいながらあんまりバイオインフォマティクス成分がなかった件。UCSC Genome Browser をRのplotで再現するということも可能。

これからはこちらを参考にBUGSをやってみようかなと思う。とりあえず

RStanを入れないといけないっぽい。

2013-11-05

MikuHatsune2013-11-05

初音ミクの流行解析をDTMで

昔、声優統計でDTMをしたのだが、その下準備に初音ミクでDTMをしようとしてしたはいいけど結果を書いてなかったので書く。

結果としてはよくわからんがこんな感じのトピックを抽出した。

 

Topic 6: 元気な歌?

,意味,とか,かご,分かる,,気がつく,se,なんか,control,大切,,考える,それでも,たり,,死ぬ,鼓動,すぎる,ちゃう,気分,くらい,カナ,),じゃう,聞く,,呼吸,気づく,こす,かも,,能力,違う,捨てる,,閉じる,文字,ミックス,醒める,),だけど,仲間,ぼる,痛み,ちゃ,キライ,冗談,強い,欲しい

Topic 8: 切ない想いを伝える別れの歌?

伝える,!!!!,一緒,,とき,会う,幼い,燃える,会える,通り過ぎる,絶対,素敵,,さようなら,朝焼け,願う,悲しい,出会う,言える,...?,思い,,呟く,,そろそろ,過ぎる,,,見つめる,続く,季節,切ない,歩む,癒す,関係,これから,届け,どんな,濡らす,願い,幾つ,楽しい,おしゃれ,,旅立つ,溢れる,仕草,約束,人形,

Topic 13: また元気な歌?

うま,,震わせる,求める,重なる,,入り込む,燃やす,,映る,溶かす,!!!,夕焼け,無限,揺れる,見つめる,草原,ヒトリ,語りかける,凍える,カラ,,帰る,近づく,,叫ぶ,動ける,ドリル,高鳴る,溢れる,願う,重ねる,おいしい,定め,ほど,沈む,お願い,濡らす,姿,目と,,決める,オレ,震える,,かご,悲しみ,一つ,白い,意味

Topic 9: みんなでドンチャンな歌?

なぁ,みんな,,ちゃう,ちょっと,,クマ,ぼる,,動画,とか,たのしい,下手,,寝る,光景,全然,飲む,こける,大切,ニコニコ,涼しい,ちゃ,展開,,いったい,たどり着く,だけど,書く,戻る,,キライ,とても,刻み付ける,なんか,たり,大根,ちゃんと,兄ちゃん,マンガ,過ぎ行く,居場所,,,でる,儚い,,すぎる,,じゃう

Topic 18: ボーカロイドの運命的な歌?

,~,,ボーカロイド,オチ,ハーモニー,!!,,ふる,リン,誰か,またがる,うまい,!」,,j,ちゃう,うた,るり,すぎ,たれる,儚い,,!(,,ブレーキ,,仲良く,強引,+,はるか,あんた,,アイドル,,よる,:,いう,w,大根,~♪,クール,,えっ,f,カラダ,押入る,兄ちゃん,みんな,混ぜる

Topic 11: 中二病的な歌?

,妊娠,革命,ちょい,教え,一度,もしかして,シマ,塗る,,何もかも,人気,,記載,,,冗談,香港,愛す,見飽きる,),使徒,遅れる,履く,,,臨界,ステレオ,,スキ,,,実体,タイムリミット,テーブル,いびつ,jealousy,超絶,点く,まもる,,転載,吹き荒れる,,,賭ける,ぎゅっと,ふらり,プルプル,からくり

Topic 2: 大自然でオシャレな歌?

輝く,...?,高い,,描く,濡らす,,日差し,進む,ふわふわ,続く,草原,おしゃれ,向こう,浮く,バラ,作れる,見上げる,っていう,希望,目指す,限り,,届ける,新しい,ファン,麦わら帽子,もうすぐ,歩む,叫べる,燃やす,ささやか,きらきら,満ちる,見つける,,,決める,さあ,溶け合う,もぐ,願い,果て,乗せる,やっと,,羽ばたく,すぐ,勇気,打ち寄せる

f:id:MikuHatsune:20131105224423p:image

 

初音ミク解析を実行してmusic_info.txtがあるとして

data0 <- read.csv(paste(wd, "music_info.txt", sep=""), header=FALSE, encoding="utf-8") # data 読み込み

submit_data <- read.delim(paste(wd, "count_info.txt", sep=""), header=FALSE)
colnames(submit_data) <- c("nicoID", "view", "comment", "mylist", "post")
submit_date <- as.Date(gsub("/", "-", submit_data[, 5]))

music_info <- read.csv(paste(wd, "music_info.csv", sep=""), header=FALSE)
colnames(music_info) <- c("ID", "title", "nicoID", "song", "music", "arrenge", "vocaloid", "lyric")
matchID <- match(submit_data$nicoID, music_info[,3])

# 統合したデータ
data1 <- cbind(submit_data, music_info[matchID,])
data1 <-  data1[order(data1$post),] # 日付順に並び替え

lyrics <- as.character(data1$lyric)
corpus <- lexicalize(lyrics)

vocab0 <- as.character(tfidf0$word[tfidf0$nscore > 0.4])
corpus0 <- lexicalize(doclines=lyrics, vocab=vocab0) #TFIDFスコアで選出した語彙
corpus <- list(documents=corpus0, vocab=vocab0)

multi_convert <- function(lex){
	paste(apply(lex, 2, paste, collapse=":"), collapse=" ")
}
input <- lapply(corpus$document, multi_convert)
word_count <- sapply(corpus$document, ncol)
multi_dat <- mapply(paste, word_count, input, sep=" ")

write.csv(matrix(unlist(multi_dat), nc=1), file=paste("miku_multi_dat.txt", sep=""), row.names=FALSE, quote=FALSE) # DTM用データ

# seq.dat の作成
days <- sort(as.Date(submit_data[,5]))
birth <- as.Date("07-08-31")
qbreak <- birth + seq(0, 6, by=1/4)*365 # 4半期
hbreak <- birth + seq(0, 6, by=1/2)*365 # 半年
ybreak <- birth + seq(0, 6, by=1)*365 # 一年
# 1行目に時系列の数、2行目以降 i 番目の時系列の数
ta <- table(cut(days, hbreak))
file0 <- "seq_half.dat"
write.table(matrix(c(length(ta), ta), nc=1), file=file0, row.names=FALSE, col.names=FALSE)

# TF-IDF っぽいスコアを計算する関数
TFIDF <- function(corpus, progress=FALSE){ # lexicalize した corpus を使用
	res <- matrix(0, nr=length(corpus$vocab), nc=4)
	dimnames(res) <- list(corpus$vocab, c("documents", "count", "freq", "score"))
	res[, "documents"] <- length(corpus$documents)
	wordset <- mapply(function(x) x[1,], corpus$documents) # documents中の単語リスト
	wordfreq <- table(unlist(wordset)) # すべての単語の、全documents中の出現頻度
	for(v in seq(corpus$vocab)){ # vocab と i は 1 ずれているので注意
		count_docs <- sum(sapply(lapply(wordset, "==", v-1), any)) # その単語が出現する文章の数
		res[v, "freq"] <- count_docs
		if(progress){ # Linux用。プログレスバーを付ける
			pb <- txtProgressBar(min=1, max=length(corpus$vocab), style=3)
			setTxtProgressBar(pb, v)
		}
	}
	res[, "count"] <- wordfreq
	res[, "score"] <- log(res[, "count"]) * log(res[, "documents"]/res[, "freq"])
	return(res)
}
tfidf0 <- TFIDF(corpus, progress=TRUE)
write.(tfidf0, file="TFIDF_score.txt")
# DTMの解析
lyrics <- as.character(data1$lyric)
time0 <- read.csv("miku-seq.dat", header=FALSE)$V1
n_times <- time0[1] # 時系列の数
timeseg <- time0[-1]
timerange <- range(as.Date(data1$post))
birth <- as.Date("07-08-31")
timeidx <- seq(birth, timerange[2], length=12)

tfidf0 <- read.csv(paste(wd, "TFIDF_score.csv", sep=""))
word0 <- as.character(read.csv(paste(wd, "score04word.txt", sep=""))$x) # TFIDFでスコア0.4以上のものを選んでいるはず


# e-log のあるディレクトリ
file0 <- list.files()[grep("e-log", list.files())]
topicdtm <- array(0, c(length(word0), n_times, length(file0))) # term, time, topic
for(i in seq(file0)){
	topicdtm[,, i] <- exp(matrix(scan(file0[i]), ncol=n_times, byrow=TRUE))
}

corpus <- lexicalize(lyrics, vocab=word0)
lda0 <- lda.collapsed.gibbs.sampler(corpus, 30, word0, 25, 0.1, 0.1)
toptopis <- top.topic.words(lda0$topics, 50, by.score=TRUE)
tmptopic <- match(toptopis[, 1], cor0$vocab)
# トピックごとの上位単語
mapply(function(x) apply(apply(topicdtm[, , x], 2, order, decreasing=TRUE), 2, head, 50), seq(n_times))

L <- list(word=array(0, c(50, n_times, 30)), prob=array(0, c(50, n_times, 30)))
for(i in seq(n_times)){
	for(j in seq(30)){
		idx <- head(order(topicdtm[, i, j], decreasing=TRUE), 50)
		L$prob[, i, j] <- topicdtm[idx, i, j]
		L$word[, i, j] <- word0[idx]
}}

# 各トピック、時系列で再生数の最も多い代表曲を確かめる
n_music <- 2 # 再生数の多い上位2曲
result <- array(0, c(n_music, 3, n_times, 30))
saisei <- 393939
for(t0 in seq(30)){
for(d0 in seq(timeseg)){
	label0 <- which(docgroup == d0)
	topmusic <- apply(e.theta[label0,], 2, order, decreasing=TRUE)
	y0 <- head(which(data1[label0,][topmusic[,t0], "view"] > saisei), n_music) # 再生数の多い上位2曲
	#print(as.character(data1[label0,][topmusic[,1], "title"][y0]))
	result[,,d0,t0] <- as.matrix(data1[label0,][topmusic[,t0], c("view", "post", "title")][y0,])
}}

# プロットする
plottopic <- 7
cols <- rainbow(plottopic)
topicdiffidx <- order(apply(apply(probsum, 2, range), 2, diff), decreasing=TRUE)
yl <- range(c(0, probsum[, head(topicdiffidx, plottopic)]))
matplot(probsum[, head(topicdiffidx, plottopic)], xlim=c(0.2, n_times+1.5), ylim=yl, col=cols, lty=1, lwd=5, type="l", xaxt="n", ylab="topic proportion", bty="l")
axis(1, at=seq(n_times), label=timeidx, las=2)
y0 <- probsum[n_times, head(topicdiffidx, plottopic)]
y0[c(4,6)] <- y0[c(4,6)] + c(0.005, -0.003)
text(rep(par()$usr[2], plottopic), y0, paste("Topic", head(topicdiffidx, plottopic)), xpd=TRUE, adj=0.7, col=cols, cex=1.2)
title("Change of topic proportion", cex.main=2)
for(t0 in head(topicdiffidx, 7)){
for(d0 in seq(n_times)){
	text(d0, probsum[d0, t0], paste(result[,3,d0,t0], collapse="\n"), cex=0.8, xpd=TRUE) # 上位の歌
}}

2013-10-01

MikuHatsune2013-10-01

時系列解析

セミナーで時系列を学んだのだが、初音ミクの投稿数に適応してみる。

Kalman filterで過去から現在までたどって、平滑化で現在から過去へ遡ることでパラメータを最適化していく。

これを使うと欠損値補完もできる。

そして前回はやらなかった、今後の投稿数予測もやる。

時系列データにはTrend(上昇気味か減少気味か), Seasonal(所定の周期における規則的変動), Cycle(任意の周期における規則的変動), Residuals(誤差)から構成されるもんだと考える。これらをそれぞれ推定して、組み合わせたら元の時系列データだと考える。

予測をしたら95%CIだとものすごいばらつく。こんなもんなのか…Trendとしては増加気味なのに予測投稿数がオワコン化しているのでこれは次数選択がミスっているのではないかという指摘を受けた。

2000くらいの観測点で6時間くらい計算に時間がかかってるのでどこか最適化したい。

f:id:MikuHatsune:20131001144438p:image

f:id:MikuHatsune:20131001174625p:image

 

y <- time.series.data #解析したい時系列データ。name()に日付が入っているとよい。
up <- var(y) #分散のオーダー
N <- length(y)
range(y)

#オプションの設定
limy <- 1000      #計測値の上限(range(y)を見て設定)
MJ <- 6          #Cycle成分:次数の最大値
p <- 12           #Seasonal成分:季節変動の周期(月次データなので12)
k <- 2            #Trend成分:次数(1次関数のトレンドを想定)
MM <- 5         #Cycle成分:係数推定のための反復計算の回数

####Preparation for maximizing likelihood and optimization(Defining functions)#####

#Cycle成分:PARCORからAR係数を計算する関数(係数決定のためのアルゴリズム)
#詳しくは「ベイズ統計データ解析」p.148-149を参照のこと
ARCOEF <- function(par,q) {
  aa <- numeric(50)
  al <- numeric(q)
  if (q == 1) {
    al <- par
  } else {
    for (II in 1:q) {
      al[II] <- par[II]
      aa[II] <- par[II]
      if (II > 1) {
        for (J in 1:(II-1)) {
          al[J] <- aa[J] - par[II]*aa[II-J]
        }
        if (II < q) {
          for (J in 1:(II-1)) {
            aa[J] <- al[J]
          }          
        }
      }
    }
  }
  return(al)
}

#全体:状態空間モデルの行列設定のための関数(資料に掲載した数式に対応)
FGHset2 <- function(al,k,p,q) {
  m <- k + p + q - 1
  if (q > 0) {G <- matrix(0,m,3)} else {G <- matrix(0,m,2)}
  F <- matrix(0,m,m)
  H <- matrix(0,1,m)
  G[1,1] <- 1
  H[1,1] <- 1
  if (k == 1) {F[1,1] <- 1}
  if (k == 2) {F[1,1] <- 2;F[1,2] <- -1;F[2,1] <- 1}
  if (k == 3) {F[1,1] <- 3;F[1,2] <- -3;F[1,3] <- 1;F[2,1] <- 1;F[3,2] <- 1}
  LS <- k
  NS <- 2
  G[LS+1,NS] <- 1
  H[1,LS+1] <- 1
  for (i in 1:(p-1)) {F[LS+1,LS+i] <- -1}
  for (i in 1:(p-2)) {F[LS+i+1,LS+i] <- 1}
  LS <- LS + p - 1
  if (q > 0) {
    NS <- NS + 1
    G[LS+1,NS] <- 1
    H[1,LS+1] <-1
    for (i in 1:q) {F[LS+1,LS+i] <- al[i]}
    if (q > 1) {
      for (i in 1:(q-1)){F[LS+i+1,LS+i] <- 1}
    }
  }
  return(list(m=m,MatF=F,MatG=G,MatH=H))
}

#状態空間モデルにおける行列Qの設定(資料に掲載した数式に対応)
Qset <- function(TAU12,TAU22,TAU32,k,p,q) {
  NS <- 0
  if (k > 0) {NS <- NS + 1}
  if (p > 0) {NS <- NS + 1}
  if (q > 0) {NS <- NS + 1}
  Q <- diag(NS)
  NS <- 0
  if (k > 0) {NS <- NS + 1; Q[NS,NS] <- TAU12}
  if (p > 0) {NS <- NS + 1; Q[NS,NS] <- TAU22}
  if (q > 0) {NS <- NS + 1; Q[NS,NS] <- TAU32}
  return(Q)
}

#カルマンフィルタの関数(資料に掲載した数式に対応)
#"%*%"は行列の積を表す記号
KF <- function(y,XFO,VFO,F,H,G,Q,R,limy,ISW,OSW,m,N) {
  if(OSW == 1) {
    XPS <- matrix(0,m,N)       #平均xn|n-1を格納する行列
    XFS <- matrix(0,m,N)       #平均xn|nを格納する行列
    VPS <- array(dim=c(m,m,N)) #分散共分散行列vn|n-1を格納する配列
    VFS <- array(dim=c(m,m,N)) #分散共分散行列vn|nを格納する配列
  }
  XF <- XFO #x0|0
  VF <- VFO #v0|0
  NSUM <- 0 #フィルタの実行回数カウンター
  SIG2 <- 0 #分散係数の最尤推定値(の初期値)
  LDET <- 0 #最大対数尤度の第2項(の初期値)
  for(n in 1:N) {
    #1期先予測
    XP <- F %*% XF                             #xn-1|n-1からxn|n-1を求める
    VP <- F %*% VF %*% t(F) + G %*% Q %*% t(G) #Vn-1|n-1からVn|n-1を求める
    #フィルタ
    if(y[n] < limy) { #y[i]が上限値(limy)を超えていない場合
      NSUM <- NSUM + 1            #フィルタの実行回数のカウント
      B <- H %*% VP %*% t(H) + R       #ynの分散共分散行列
      B1 <- solve(B)             #Bの逆行列を求める
      K <- VP %*% t(H) %*% B1           #Kはカルマンゲイン
      e <- y[n] - H %*% XP          #ynの予測誤差
      XF <- XP + K %*% e           #xn|n-1とynの予測誤差からxn|nを求める
      VF <- VP - K %*% H %*% VP        #Vn|n-1からVn|nを求める
      SIG2 <- SIG2 + t(e) %*% B1 %*% e #SIG2の最尤推定値
      LDET <- LDET + log(det(B))        #最大対数尤度の第2項の計算
    } else { #y[i]が上限値(limy)を超えている場合
      XF <- XP                         #xn|n=xn|n-1とする
      VF <- VP                         #vn|n=vn|n-1とする
    }
    if (OSW == 1) { #OSW==1:状態フィルタ分布を計算する場合
      XPS[,n] <- XP  #xn|n-1を格納
      XFS[,n] <- XF  #xn|nを格納
      VPS[,,n] <- VP #vn|n-1を格納
      VFS[,,n] <- VF #vn|nを格納
    }
  }
  SIG2 <- SIG2/NSUM
  if (ISW == 0) {  
    FF <- -0.5 * (NSUM * (log(2*pi*SIG2)+1)+LDET) #SIG2を未知パラメータとして最尤推定する場合の最大対数尤度
  } else {
    FF <- -0.5 * (NSUM * (log(2*pi)+SIG2)+LDET)   #SIG2を所与の値とした場合の最大対数尤度
  }
  if (OSW == 0) { #OSW==0:最大対数尤度を計算する場合
    return(list(LLF=FF,Ovar=SIG2))
  }
  if (OSW == 1) { #OSW==1:状態フィルタ分布を計算する場合
    return(list(XPS=XPS,XFS=XFS,VPS=VPS,VFS=VFS,Ovar=SIG2))
  }
}

#平滑化の関数(資料に掲載した数式に対応)
SMO <- function(XPS,XFS,VPS,VFS,F,GSIG2,k,p,q,m,N) {
  XSS <- matrix(0,m,N)     #平滑化の結果のうち,Xn|Nを格納する入れ物
  VSS <- array(dim=c(m,m,N))  #平滑化の結果のうち,Vn|Nを格納する入れ物
  XS1 <- XFS[,N]  #xN|NをXS1に代入
  VS1 <- VFS[,,N] #VN|NをVS1に代入
  XSS[,N] <- XS1
  VSS[,,N] <- VS1
  for (n1 in 1:(N-1)) {
    n <- N - n1      #N-n
    XP <- XPS[,n+1]  #XN-n+1|N
    XF <- XFS[,n]    #XN-n|N
    VP <- VPS[,,n+1] #VN-n+1|N
    VF <- VFS[,,n]   #VN-n|N
    VPI <- solve(VP) #VN-n+1|Nの逆行列
    A <- VF %*% t(F) %*% VPI              #AN = Vn|n * F^t_n+1 * F^-1n+1|nに対応
    XS2 <- XF + A %*% (XS1 - XP)          #Xn|N = xn|n + An (xn+1|N - xn+1|n)に対応
    VS2 <- VF + A %*% (VS1 - VP) %*% t(A) #Vn|N = Vn|n + An (Vn+1|N - Vn+1|n)*A^t_nに対応
    XS1 <- XS2 #次の回に備えてXS1を作っておく
    VS1 <- VS2 #次の回に備えてVS1を作っておく
    XSS[,n] <- XS1  #Xn|Nを格納
    VSS[,,n] <- VS1  #Vn|Nを格納
  }
  t <- numeric(N)     #tnの係数の平均値の収束値を格納するスペース
  s <- numeric(N)   #snの係数の平均値の収束値を格納するスペース
  r <- numeric(N)     #rnの係数の平均値の習得値を格納するスペース
  tv <- numeric(N)    #tnの係数の分散の収束値を格納するスペース
  sv <- numeric(N)  #snの係数の分散の収束値を格納するスペース
  rv <- numeric(N)    #rnの係数の分散の収束値を格納するスペース
  if (q == 0) {
    for (n in 1:N) {
      t[n] <- XSS[1,n]
      s[n] <- XSS[k+1,n]
      tv[n] <- GSIG2 * VSS[1,1,n]
      sv[n] <- GSIG2 * VSS[k+1,k+1,n]
    }
  } else {
    for (n in 1:N) {
      t[n] <- XSS[1,n]
      s[n] <- XSS[k+1,n]
      r[n] <- XSS[k+p,n]
      tv[n] <- GSIG2 * VSS[1,1,n]
      sv[n] <- GSIG2 * VSS[k+1,k+1,n]
      rv[n] <- GSIG2 * VSS[k+p,k+p,n]
    }
  }
  return(list(trd=t,sea=s,arc=r,trv=tv,sev=sv,arv=rv))
}

#カルマンフィルタおよび平滑化の関連計算
TSRest <- function(y,TAU12,TAU22,TAU32,al,OSW,limy,k,p,q,N) {
  MAT <- FGHset2(al,k,p,q)
  m <- MAT$m
  F <- MAT$MatF
  G <- MAT$MatG
  H <- MAT$MatH
  ISW <- 0
  Q <- Qset(TAU12,TAU22,TAU32,k,p,q)
  R <- diag(1)
  XFO <- numeric(m)
  VFO <- 100*diag(m)
  for (i in 1:k) {XFO[i] <- y[1]}
  if (q > 0) {
    for (i in 1:q) {
      VFO[k+p+i-1,k+p+i-1] <- 10
    }
  }
  LLF <- KF(y,XFO,VFO,F,H,G,Q,R,limy,ISW,OSW,m,N)
  if (OSW == 1) {
    XPS <- LLF$XPS
    XFS <- LLF$XFS
    VPS <- LLF$VPS
    VFS <- LLF$VFS
    SIG2 <- LLF$Ovar
    LL - LLF$LLF
    XVS <- SMO(XPS,XFS,VPS,VFS,F,SIG2,k,p,q,m,N)
    xt <- XVS$trd
    xs <- XVS$sea
    xr <- XVS$arc
    return(list(LLF=LL,xt=xt,xs=xs,xr=xr,SIG2=SIG2,F=F,G=G,H=H,Q=Q,R=R,
                XPS=XPS,XFS=XFS,VPS=VPS,VFS=VFS))
  } else {
    LL <- LLF$LLF
    SIG2 <- LLF$Ovar
  }
  return(list(LLF=LL,SIG2=SIG2))
}

#Trend成分:TAU1^2(分散)の尤度関数を最大化する関数
#計算を容易にするため対数をとっています
LogLI1 <- function(theta,y,limy,k,p,N) {
  TAU12 <- theta
  #TAU2^2に関する最適化
  #オーダーを踏まえて、upperを指定
  LLF <- optimize(LogLI2,lower=1e-2,upper=up,maximum=TRUE,y=y,
                  TAU12=TAU12,limy=limy,k=k,p=p,N=N)
  LL <- LLF$objective
  return(LL)
}

#Seasonal成分:TAU2^2(分散)の尤度関数を最大化する関数
#上記でTAU1^2が決定した後にこの関数を用います
LogLI2 <- function(theta,y,TAU12,limy,k,p,N) {
  TAU22 <- theta
  LLF <- TSRest(y,TAU12,TAU22,0,0,0,limy,k,p,0,N)
  LL <- LLF$LLF
  return(LL)
}

#Cycle成分:TAU3^2(分散)の尤度関数を最大化する関数
#詳しくは「ベイズ統計データ解析」p.148-149を参照のこと
#上記でTAU1^2,TAU^2が決定した後にこの関数を用います
LogL1 <- function(theta,y,TAU12,TAU22,Spar,limy,k,p,q,N) {
  TAU32 <- theta
  #PARCOR(偏自己相関係数)に関する最適化
  if (q == 1) {
    LLF <- optimize(LogL2,lower=-0.95,upper=0.95,maximum=TRUE,y=y,
                    TAU12=TAU12,TAU22=TAU22,TAU32=TAU32,limy=limy,k=k,p=p,N=N)
  } else {
    LLF <- optimize(LogL3,lower=-0.95,upper=0.95,maximum=TRUE,y=y,
                    TAU12=TAU12,TAU22=TAU22,TAU32=TAU32,Spar=Spar,limy=limy,
                    k=k,p=p,q=q,N=N)
  }
  LL <- LLF$objective
  return(LL)
}

#Cycle成分:ARモデルの係数αを決定するためのアルゴリズム
#上記でTAU1^2,TAU2^2,TAU3^2が決定した後にこの関数を用います
#アルゴリズムの詳細は「ベイズ統計解析」p.148-149を参照のこと
LogL2 <- function(theta,y,TAU12,TAU22,TAU32,limy,k,p,N) {
  al <- numeric(1)
  al[1] <- theta
  LLF <- TSRest(y,TAU12,TAU22,TAU32,al,0,limy,k,p,1,N)
  LL <- LLF$LLF
  return(LL)
}

#Cycle成分:ARモデルの係数αを決定するためのアルゴリズム
#アルゴリズムの詳細は「ベイズ統計解析」p.148-149を参照のこと
#(q > 1)次以上のPARCORの対数尤度
#(q-1)次までのPARCOR,TAU1^2,TAU2^2,TAU3^2は所与
LogL3 <- function(theta,y,TAU12,TAU22,TAU32,Spar,limy,k,p,q,N) {
  par <- numeric(q)
  par[1:(q-1)] <- Spar[1:(q-1)]
  par[q] <- theta
  al <- ARCOEF(par,q)
  LLF <- TSRest(y,TAU12,TAU22,TAU32,al,0,limy,k,p,q,N)
  LL <- LLF$LLF
  return(LL)
}

#Cycle成分:ARモデルの係数αを決定するためのアルゴリズム
#アルゴリズムの詳細は「ベイズ統計解析」p.148-149を参照のこと
#TAU1^2の対数尤度(al,TAU3^2は所与)
LogL4 <- function(theta,y,TAU32,al,limy,k,p,q,N) {
  TAU12 <- theta
  #オーダーを踏まえて、upperを指定
  LLF <- optimize(LogL5,lower=1e-2,upper=up,maximum=TRUE,y=y,
                  TAU12=TAU12,TAU32=TAU32,al=al,limy=limy,k=k,p=p,q=q,N=N)
  LL <- LLF$objective
  return(LL)
}

#Cycle成分:ARモデルの係数αを決定するためのアルゴリズム
#アルゴリズムの詳細は「ベイズ統計解析」p.148-149を参照のこと
#TAU2^2の対数尤度(TAU1^2,TAU3^2は所与)
LogL5 <- function(theta,y,TAU12,TAU32,al,limy,k,p,q,N) {
  TAU22 <- theta
  LLF <- TSRest(y,TAU12,TAU22,TAU32,al,0,limy,k,p,q,N)
  LL <- LLF$LLF
  return(LL)
}

####Executing algorithm maximizing likelihood and optimization#####

#Trend成分:TAU1^2の最尤推定
#オーダーを踏まえて、upperを指定
LLF1 <- optimize(LogLI1,lower=1e-2,upper=up,maximum=TRUE,y=y,limy=limy,k=k,p=p,N=N) 
ITAU12 <- LLF1$maximum

#Sesonal成分:TAU2^2の最尤推定
#オーダーを踏まえて、upperを指定
LLF2 <- optimize(LogLI2,lower=1e-2,upper=up,maximum=TRUE,y=y,TAU12=ITAU12,limy=limy,k=k,p=p,N=N)

ITAU22 <- LLF2$maximum

STAU12 <- numeric(MJ)     #トレンド成分モデルの誤差分散の保存スペース
STAU22 <- numeric(MJ)     #季節成分モデルの誤差分散の保存スペース
STAU32 <- numeric(MJ)     #AR成分モデルの誤差分散の保存スペース
SSIG2 <- numeric(MJ)      #観測ノイズの分散の保存スペース
SLL <- numeric(MJ)        #最大対数尤度の保存スペース
SAIC <- numeric(MJ)       #AICの保存スペース
Sal <- matrix(0,MJ,MJ)    #AR成分モデルの係数の保存スペース(各行毎)
Spar <- numeric(MJ)       #PARCORの保存スペース
MAIC <- 1e10              

#Cycle成分:ARモデルの次数選択,係数αの最尤推定
for (J in 1:MJ) {
  q <- J
  PAR <- numeric(q)
  for (II in 1:MM) {
    if ((q == 1)&&(II == 1)) {TAU12 <- ITAU12; TAU22 <- ITAU22}
    #TAU3^2の最尤推定
    #オーダーを踏まえて、upperを指定
    LLF3 <- optimize(LogL1,lower=1e-4,upper=up,maximum=TRUE,y=y,
                     TAU12=TAU12,TAU22=TAU22,Spar=Spar,limy=limy,
                     k=k,p=p,q=q,N=N)
    TAU32 <- LLF3$maximum
    #PARCORに関する最適化(係数αの最尤推定)
    if (q == 1) {
      LLF4 <- optimize(LogL2,lower=-0.95,upper=0.95,maximum=TRUE,y=y,
                       TAU12=TAU12,TAU22=TAU22,TAU32=TAU32,limy=limy,
                       k=k,p=p,N=N)
      par <- LLF4$maximum
      PAR[1] <- par
      al <- numeric(1)
      al[1] <- par
    } else {
      LLF4 <- optimize(LogL3,lower=-0.95,upper=0.95,maximum=TRUE,y=y,
                       TAU12=TAU12,TAU22=TAU22,TAU32=TAU32,Spar=Spar,
                       limy=limy,k=k,p=p,q=q,N=N)
      par <- LLF4$maximum
      PAR[1:(q-1)] <- Spar[1:(q-1)]
      PAR[q] <- par
      al <- ARCOEF(PAR,q)
    }
    #TAU1^2に関する最適化
    #オーダーを踏まえて、upperを指定
    LLF5 <- optimize(LogL4,lower=1e-2,upper=up,maximum=TRUE,y=y,
                     TAU32=TAU32,al=al,limy=limy,k=k,p=p,q=q,N=N)
    TAU12 <- LLF5$maximum
    #TAU2^2に関する最適化
    #オーダーを踏まえて、upperを指定
    LLF6 <- optimize(LogL5,lower=1e-2,upper=up,maximum=TRUE,y=y,
                     TAU12=TAU12,TAU32=TAU32,al=al,limy=limy,k=k,p=p,q=q,N=N)
    TAU22 <- LLF6$maximum
    print(paste("J =",J,",II =",II,"done"))
  }
  LLF7 <- TSRest(y,TAU12,TAU22,TAU32,al,0,limy,k,p,q,N)
  LL <- LLF7$LLF
  SIG2 <- LLF7$SIG2
  AIC <- -2*LL + 2*(q+4)
  if (AIC < MAIC) {
    MAIC <- AIC     #AICの最小値を更新
    Mq <- q     #最小のAICを与える次数qを更新
  }
  Sal[q,1:q] <- al
  Spar[q] <- PAR[q]
  STAU12[q] <- TAU12
  STAU22[q] <- TAU22
  STAU32[q] <- TAU32
  SSIG2[q] <- SIG2
  SLL[q] <- LL
  SAIC[q] <- AIC
  print(paste("J =",J,"done"))
} #JはCycle成分の次数,IIはTAU最尤推定に関するアルゴリズムの実行回数

TAU12 <- STAU12[Mq]   #Trend成分:最小のAICを与えるTAU1^2
TAU22 <- STAU22[Mq]   #Seasonal成分:最小のAICを与えるTAU2^2
TAU32 <- STAU32[Mq]   #Cycle成分:最小のAICを与えるTAU3^2
al <- Sal[Mq,1:Mq]    #Cycle成分:最小のAICを与えるARモデルの係数
LLF8 <- TSRest(y,TAU12,TAU22,TAU32,al,1,limy,k,p,Mq,N)
xt <- LLF8$xt         #Trend成分
xs <- LLF8$xs         #Seasonal成分
xr <- LLF8$xr         #Cycle成分
w <- y - xt - xs - xr #観測ノイズ(residuals)
yad <- y - xs         #季節調整済み系列

#モデル推定の結果表示
print(SLL)   #モデル全体:最大対数尤度(左から順にCycle成分の次数が1,2,…,6の場合)
print(SAIC)    #モデル全体:AIC(左から順にCycle成分の次数が1,2,…,6の場合)
print(STAU12)  #Trend成分:ノイズ項の分散(左から順にCycle成分の次数が1,2,…,6の場合)
print(STAU22)  #Sesonal成分:ノイズ項の分散(左から順にCycle成分の次数が1,2,…,6の場合)
print(STAU32)  #Cycle成分:ノイズ項の分散(左から順にCycle成分の次数が1,2,…,6の場合)
print(SSIG2)   #Residual成分:分散(左から順にCycle成分の次数が1,2,…,6の場合)
print(Sal)     #Cycle成分:係数(1行目から順にCycle成分の次数が1,2,…,6の場合)
print(Spar)    #Cycle成分:係数推定のために用いた偏自己相関係数(左から順にCycle成分の次数が1,2,…,6の場合)
print(Mq)      #Cycle成分の次数の最終結果(AICが最も小さい次数が選択される)

#トレンドと時系列データの折れ線グラフの作成
t <- c(1:N)
par(mfrow=c(4,1), mar=c(2, 4, 0.5, 2))
plot(y, type="l", xlab="", ylab="Posted musics with trend", xaxt="n")
lines(t, xt, col=3, lwd=3)
axis(1, at=t, labels=names(y), tick=FALSE)

plot(t, xs, type="l", xlab="", ylab="Seasonal", xaxt="n", lwd=1)
axis(1, at=t, labels=names(y), tick=FALSE)

plot(t, xr, type="l", xlab="", ylab="Cycle", xaxt="n", lwd=1)
axis(1, at=t, labels=names(y), tick=FALSE)

plot(t,  w, type="l", xlab="", ylab="Residuals", xaxt="n", lwd=1)
axis(1, at=t, labels=names(y), tick=FALSE)


#####Prediction#####
#予測を行う関数
KFP <- function(XFO,VFO,F,H,G,Q,R,m,N) {
  XPSP <- matrix(0,m,N)       #平均xn|n-1を格納する行列
  XFSP <- matrix(0,m,N)       #平均xn|nを格納する行列
  VPSP <- array(dim=c(m,m,N)) #分散共分散行列vn|n-1を格納する配列
  VFSP <- array(dim=c(m,m,N)) #分散共分散行列vn|nを格納する配列
  XF <- XFO                   #予測時の初期値(xN|N)
  VF <- VFO                   #予測時の初期値(vN|N)
  ypred <- rep(0,N)           #yの予測値を格納するスペース
  ysig2 <- rep(0,N)           #yの予測値を格納するスペース
  for(n in 1:N) {
    XP <- F %*% XF                             #xn-1|n-1からxn|n-1を求める
    VP <- F %*% VF %*% t(F) + G %*% Q %*% t(G) #Vn-1|n-1からVn|n-1を求める
    ypred[n] <- H %*% XP                       #xn|n-1からynを求める
    ysig2[n] <- H %*% VP %*% t(H) + R        #Vn|n-1からUn|n-1を求める
    XF <- XP                                   #xn|n=xn|n-1とする
    VF <- VP                                   #vn|n=vn|n-1とする
    XPSP[,n] <- XP
    XFSP[,n] <- XF
    VPSP[,,n] <- VP
    VFSP[,,n] <- VF
  }
  return(list(XPSP=XPSP,XFSP=XFSP,VPSP=VPSP,VFSP=VFSP,ypred=ypred,ysig2=ysig2))
}

#100日分を予測する
nmonth <- 100
Pred <- KFP(XFO=LLF8$XFS[,length(y)],VFO=LLF8$VFS[,,length(y)],
            F=LLF8$F,H=LLF8$H,G=LLF8$G,Q=LLF8$Q,R=LLF8$R,m=k+p+Mq-1,N=nmonth)
yhat <- Pred$ypred
sig <- sqrt(Pred$ysig2 + SSIG2[Mq]) #Residuals成分の分散も合わせている(使用上は吟味すること)
yhatl <- Pred$ypred + 2*sig
yhatu <- Pred$ypred - 2*sig
yhatl2 <- Pred$ypred + sig
yhatu2 <- Pred$ypred - sig
t <- c(1:length(y))
t_p <- c((length(y)+1):(length(y)+nmonth))
plot(t,y,xlim=c(1,t_p[nmonth]),ylim=c(min(y),max(yhatl,y)),xaxt="n",
     type="l",xlab="",ylab="Posted musics",lwd=1)
lines(t, xt, col=2, lwd=3)
axis(1, at=t, labels=names(y), tick=FALSE)

lines(t_p,yhat,col=2,lwd=3)
lines(t_p,yhatl2,lty=1,col=3,lwd=3)
lines(t_p,yhatu2,lty=1,col=3,lwd=3)
lines(t_p,yhatl,lty=1,col=4,lwd=3)
lines(t_p,yhatu,lty=1,col=4,lwd=3)
legend("topleft",legend=c("Prediction","70%CI","95%CI"),col=2:4,lty=1,lwd=3)