『ラブライブ!サンシャイン!!』 の手作り料理
世間が最も関心を寄せたのは・・・
シャイ煮?ヨキソバ?堕天使の涙?

ZIPfaceSheet.png
author: Unadon (見習い飯炊き兵)
動作環境:Mac OS Sierra 10.12.1; R version3.3.1; rstan 2.10.1
 
アニメに関するキャプチャ画像は、本稿内容の紹介に必要最低限な範囲で使用するものとして利用させていただきます。

はじめに

Stanアドベントカレンダー2016、9日目。

http://qiita.com/advent-calendar/2016/stan

入門者に向けた丁寧な記事。ハイレベルな技術面の記事、テーマ性のある面白い記事。盛り上がっておりますねアドベントカレンダー。

5日目の「Stanでパイの実」は特に笑いました。お気に入りです。

<本稿のねらい>

・GoogleTrendsからRでデータを取ってこれるようになる。

・そのデータをStanにわたして正規分布で平均値の推定が出来るようになる。

・ゼロ過剰ポアソン分布で「事象が生起した時」の「平均値」が推定できるようになる。

・Stanでアニメの舞台(静岡県沼津市)を盛り上げたい。

研究の背景

Stanユーザーのみなさまは、『ラブライブ!サンシャイン!!』というアニメーション作品をご存知でしょうか。

絶大な人気を誇るシリーズ「ラブライブ!」の続編作品であり、このアニメの放映により、舞台である静岡県沼津市に革命的ともいえる現象が起こっています。

ーーー以前はもみじマーク車ばかり

・・・ 今何が走ってんねん・・・きになるなぁ

そんな、沼津の近況(アニメに影響を受けたローカルな社会・経済状況)を細やかなところまで知りたい方はこのあたりをどうぞ。

リンク先設定のエラーを修正しました2016.12.09.12:30

http://www.slideshare.net/izunyan/2-67543738

さて、そんな大人気アニメの本編 第10話では、 ヒロインたちが海の家で手作り料理を売り出します。 まずは、本編に登場した3つの料理をご紹介しましょう。

(A)シャイ煮: 高級食材の洋風煮込み

Shini 2.png

(B)ヨキソバ: 定番オムそば

YokiSoba.png

(C)堕天使の涙: 呪術的なたこ焼き

Datenshi_no_Namida.png

研究の目的

これらの料理、アニメファンにとっては「ぜひとも一度食べてみたい!夢の料理」であることとおもいます。なにせ、ヒロインたちの手作りですからね。 ところで、上記3つの料理って、どれが最も世間的な注目を集めたのでしょうか。

本研究の目的は「ヨキソバ、堕天使の涙、シャイ煮の人気度を推定する」こととします。

私は西洋料理が好きなので、『シャイ煮』を応援していきたいところです。

それでは、Stanユーザー、アニメファン、舞台沼津市が更に盛り上がらんことを祈って、

データ分析 スタートです!  

 

方法1 : Scraping

まずは、Google Trendsで検索したときの、「キーワードの人気度」をスクレイピングする準備です。 以下のサイトを参照させていただきました。使用パッケージは”gtrendsR”。

https://www.karada-good.net/analyticsr/r-392/N

あわせて、stanやデータ成形用パッケージも読み込んでおきます。

library(rstan) 
library(gtrendsR) # Scraping Google Trends
library(ggplot2) 
library(reshape2) 
library(stringi)
library(stringr)

続いて、gtrendsRを用いてRからGoogleTrendsにアクセスして検索をかけ、データを引っ張ってきます。

皆様のGoogleアカウントアドレスとパスワード、さらには検索キーワードを下記のコード内に入れて実行してください。

(解析時に複数アカウントでログインしている場合などはうまくいかないことがあります) (パスワードを含んだコードをGitクローンしないよう気をつけてくださいね。)
#グーグルアカウントの "メール" と"パスワード"を入力
usr = " " #mail adress
psw = " " # PassWord

# 検索キーワードを設定
keyWords<-c("シャイ煮", "ヨキソバ", "堕天使の涙")

# データ取得期間の設定。機関3ヶ月未満だと日毎のデータ、それ以上は週ごとになってしまう
startDate = c("2016-08-01") #アニメ放送3週間くらい前、たぶん 
endDate = c("2016-11-01") #今回は3ヶ月間の情報を取得


#login
gconnect(usr, psw)


#Scraping with gtrends
GetTrend = gtrends(query = stri_encode(keyWords, "", "utf-8"), 
                    geo = "JP", # get data of "Japan Regions: JP"
                    start_date = startDate, 
                    end_date = endDate
)

# ggplot2-------------------------------------------

#以下使用する全てのggplotにtheme_bwと日本語表記可能な設定を適用する
ggplot()+theme_set(theme_bw(base_family="HiraKakuProN-W3"))

#geom_line
Dat = as.data.frame(GetTrend$trend)
ggplot(Dat, aes(x = start, y = hits, colour = keyword))+
  geom_line(size = 1.0) +labs(colour = "検索ワード", x = "日付(月)", y = "検索人気度", title = GetTrend$meta) 

データを取り出してくることができました。Plotしていきましょう。

本稿では、”Plotly”で可視化しています。

ライン上にカーソルをあわせて、データがどんなもんかご確認ください。

グラフ右側の各検索ワードをクリック(タップ)すると、当該データの表示/非表示が切り替えられます。

方法2 : 正規分布による平均の推定

さて、上記データから、”それぞれの人気度平均”を求めていきたいと思います。

まずは、3つの料理それぞれの人気度が正規分布に従うと仮定して、人気度平均を求めてみます。 平均パラメータμには、それぞれの料理名をとりました。

# Data ----------------------------------------------------
#上の検索結果Datから、料理別でヒット数を抜いてくる
Shini<-Dat$hits[1:c(length(Dat$hits)/3)]
Datenshi<-Dat$hits[c(length(Dat$hits)/3+1):(2*c(length(Dat$hits)/3))]
Yoki<-Dat$hits[(2*c(length(Dat$hits)/3)+1):c(length(Dat$hits))]
N<-c(length(Dat$hits)/3)

#For Stan----------------------------------------------------------
#shini, datenshi, yokiをそれぞれGaussianの平均パラメータとする
parameters<- c("shini","datenshi","yoki","sigma")

# 各料理の検索人気は[0,100]の整数データ。Nは料理別での取得データ数。
#これをStanに渡すデータとする。
datastan=list(Shini=Shini,Datenshi=Datenshi,Yoki=Yoki,N=N)

#Gaussian Means
model = '
data {
int N;
int Shini[N];
int Datenshi[N];
int Yoki[N];
}

parameters{
real <lower=0, upper=100>shini;
real <lower=0, upper=100>datenshi;
real <lower=0, upper=100>yoki;
vector <lower=0>[3]sigma;

}

model {
//Gaussian means and SD
shini ~ normal(0, 100);
datenshi ~ normal(0, 100);
yoki ~ normal(0, 100);
sigma ~ cauchy(0,25);

//Mean Estimation
for(i in 1:N){
Shini[i] ~ normal(shini, sigma[1]); 
Datenshi[i] ~ normal(datenshi, sigma[2]); 
Yoki[i] ~ normal(yoki, sigma[3]); 
}
}

'

#Stanfit
fit = stan(model_code = model,data=datastan,pars = parameters)

(´Д`)ハァハァ 
(・・・Stanすると(´Д`)ハァハァっていうの、どなたが最初に言い始めたものなのでしょうか・・・きになります)

でました。推定結果です。

正規分布を仮定した際の、検索人気度の平均および標準偏差は以下の通り。

print(summary(fit)$summary[,c(1,4,8,9,10)])
                 mean         2.5%       97.5%    n_eff      Rhat
shini       5.5502215    3.0395662    8.068889 3583.988 0.9999276
datenshi    1.3717772    0.8449075    1.894418 2680.855 1.0000612
yoki        0.8342583    0.3712833    1.278270 2129.883 0.9998592
sigma[1]   12.2028121   10.5247844   14.216162 3400.325 1.0006626
sigma[2]    2.5360877    2.1903405    2.933439 3057.129 0.9996223
sigma[3]    2.2331777    1.9276959    2.597941 3663.233 0.9997976
lp__     -519.3415320 -524.1032132 -516.729989 1685.838 1.0000500

Plotして見てみましょう。

#extract Samples
d.ext = as.data.frame(extract(fit,permuted=T))  #rstan::extract()のほうがいいかも
shini = d.ext$shini[1:4000]
datenshi<-d.ext$datenshi[1:4000]
yoki = d.ext$yoki[1:4000]

#result Plot 
extData<-data.frame(shini=shini,datenshi=datenshi,yoki=yoki)
extDf<-melt(extData,id=c(),variable.name = "Dishes")
extDf$value<-as.numeric(extDf$value)
ggplot(data=extDf, aes(x=value, fill=Dishes)) + geom_density(alpha=0.4,position = "identity")

この結果をみると、ヨキソバと堕天使の涙より、シャイ煮のほうが注目を集めていたということがわかりますね。

方法3 : ゼロ過剰ポアソン分布のポアソンの部分で人気度平均推定

いやぁ、シャイ煮は強いですね。

しかし、ちょっとまった。

今回採用したデータをよくよく見てみますと、検索人気度ゼロの日が結構多いことがわかります。

悲しいかな、ヒロインたちの手作り料理をキーワードとした検索は、コンスタントに生起しない事象のようなのです。

ゼロが多いと、正規分布の平均は、当然の事ながらゼロの値に大きな影響を受けます。

さらに、本編放送前のデータが入ってしまっていることもあいまって、このゼロたちの影響はちょっと取り除いて考えたい。

「検索があった時」にどれくらい注目を集めていたのか、を考えることにしましょう。

ここからは、まず「検索があるかないか」を切り分けた上で、

「検索があった時に どの程度注目が集まったのか」考えることができる、

ゼロ過剰ポアソン分布(Zero-Inflated Poisson)で推定していきたいと思います(詳細は検索していただけると幸いです)。

「検索があった日」を切りだしてポアソン分布にわたし、その平均パラメータで「検索があった日に平均的にどの程度の人気度があったのか」を求めていきます。

平均と分散が等しいという強い仮定をもつポアソン分布なので、 上記のヒストグラムを見る限り、特に『シャイ煮』の事後予測がうまくなせない可能性がありますが。

まぁ勉強です。とりあえずやってみましょう。

各料理のポアソン平均lambdaは、わかりやすいように料理名で示します。

#最初の方の検索結果Datから、料理別でヒット数を一応もう一回抜いてくる
Shini<-Dat$hits[1:c(length(Dat$hits)/3)]
Datenshi<-Dat$hits[c(length(Dat$hits)/3+1):(2*c(length(Dat$hits)/3))]
Yoki<-Dat$hits[(2*c(length(Dat$hits)/3)+1):c(length(Dat$hits))]
N<-c(length(Dat$hits)/3)

#For Stan----------------------------------------------------------
#shini, datenshi, yokiをそれぞれポワソン平均パラメータとする
#vector [3]thetaはそれぞれの料理の検索人気度がゼロじゃない確率。
parameters<- c("shini","datenshi","yoki","theta")

# 各料理の検索人気は[0,100]の整数データ。Nは料理別での取得データ数。
#これをStanに渡すデータとする。
datastan=list(Shini=Shini,Datenshi=Datenshi,Yoki=Yoki,N=N)


#zero-inflated poisson estimation model
model<-'
data {
int N;
int Shini[N];
int Datenshi[N];
int Yoki[N];
}

parameters{
real <lower=0, upper=100>shini;
real <lower=0, upper=100>datenshi;
real <lower=0, upper=100>yoki;
vector <lower=0, upper=1>[3]theta;

}

model {
for(i in 1:N){
theta[1] ~ beta(1, 1);
shini ~ uniform(0, 100);
theta[2] ~ beta(1, 1);
datenshi ~ uniform(0, 100);
theta[3] ~ beta(1, 1);
yoki ~ uniform(0, 100);

if (Shini[i] == 0) {
target+=log((1 - theta[1]) + theta[1] * exp(-shini));
} else {
target+=bernoulli_lpmf(1| theta[1])+ poisson_lpmf(Shini[i]|shini);
}

if (Datenshi[i] == 0) {
target+=log((1 - theta[2]) + theta[2] * exp(-datenshi));
} else {
target+=bernoulli_lpmf(1| theta[2])+ poisson_lpmf(Datenshi[i]| datenshi);
}
if (Yoki[i] == 0) {
target+=log((1 - theta[3]) + theta[3] * exp(-yoki));
} else {
target+=bernoulli_lpmf(1| theta[3])+ poisson_lpmf(Yoki[i]|yoki);
}
}
}

'

#Stanfit
fit<- stan(model_code = model,data=datastan,pars=parameters)

(´Д`)ハァハァ

推定結果は以下の通りです。

summary(fit)$summary[,c(1,3,4,8,10)]
>                  mean         sd          2.5%         97.5%      Rhat
>shini        9.9035930 0.43777692     9.0741888    10.7674637 0.9996099
>datenshi     4.8692234 0.44152881     4.0312414     5.7577792 1.0005611
>yoki         5.1567139 0.60543839     4.0145961     6.3688178 0.9993616
>theta[1]     0.5208945 0.03000860     0.4626146     0.5777840 1.0001533
>theta[2]     0.4289083 0.03034486     0.3698290     0.4879257 0.9995255
>theta[3]     0.3894506 0.02979860     0.3309735     0.4486778 0.9994020
>lp__     -1017.7564495 1.76691621 -1022.1444147 -1015.3761634 1.0010550

シャイ煮と、他の料理との人気度の差がよりはっきりと推定されたかんじでしょうか。

図でみてみたいので、先ほどと同じggplotのコードを使ってプロットしてみましょう(本稿ではplotlyにしてます)

#extract Samples
d.ext = as.data.frame(extract(fit,permuted=T))  #rstan::extract()のほうがいいかも
shiniZIP = d.ext$shini[1:4000]
datenshiZIP<-d.ext$datenshi[1:4000]
yokiZIP = d.ext$yoki[1:4000]

#result Plot 
extData<-data.frame(shini=shiniZIP,datenshi=datenshiZIP,yoki=yokiZIP)
extDf<-melt(extData,id=c(),variable.name = "Dishes")
extDf$value<-as.numeric(extDf$value)
ggplot(data=extDf, aes(x=value, fill=Dishes)) + geom_density(alpha=0.4,position = "identity")

ゼロ過剰ポアソン推定結果

結果からの考察

正規分布でも、ゼロ過剰ポワソンで考えた場合でも、シャイ煮の人気度は他の料理よりも高い!という結果になりました。

95%確信(信用)区間をみても、やはりシャイ煮は他の料理より世間の注目を集めていた、といえそうです。

本格的にシャイ煮を再現して売り出しているケースは今のところないようですし・・・(いくつか素敵な試みはある。)

「シャイ煮」を売り出すことで、聖地沼津はさらに盛り上がるのではないでしょうか!? 

議論

さて、Stanユーザーのみなさまともあれば、上記の提案にたいしていくつかの疑問を持たれることと思います。

Q. まずは、「シャイ煮っていう他の料理とか同名の何かがあるんじゃないの? 検索人気度は、それが混じってたんじゃないの?」というご指摘でしょう。

その点については、、Google検索、フランス料理大辞典ラルース、私が信頼する論文データベースPubMedとGoogleScholarを検索しましたがヒットしませんでしたので、たぶん大丈夫とおもいます。堕天使の涙は、宝塚歌劇団の公演名の一部に含まれていました。

Q. 「よく検索されたからといって、『食べたい』と思ってるわけではないんじゃないか?」というご指摘も考えられます。

もしそうなら、確かに売れません。

その点にかんしては、Rパッケージ”TwitteR”により「シャイ煮」をキーワードとした過去の772ツイートをテキストマイニング(共起分析)した結果でお答えいたします。

ちゃんとデータクリーニングしなかったので、余計なものが混じってますが。

TextMinningResult.png

ラブライブ!ファンのみなさまは、シャイ煮もヨキソバも堕天使の涙も、「食べたいやつはリツイート!!!」と、欲望を拡散している模様でございます。

これが日本の現状です。

Q. 「シャイ煮を売り出したところで、本当に関心を持っているのは一部のファンなんじゃないの?それで売上見込めるの?」という実利的な観点からの疑問を持たれるからもいらっしゃるでしょう。

Google Trendの地域別検索人気度ランキングを以下に示します。

まずは、ヨキソバと堕天使の涙。

YokiDatenshi.png

上記の結果は、関東に在住する一部のコアなファンたちからの関心が強いことをうかがわせます。

一方、シャイ煮は・・・といいますと・・・

MAPshi_ni.png

もはや関心は全国区なのであります。

結論

シャイ煮は人気。みんな食べたそうだから、売りましょう。きっと全国からのお客さんが増えることと思います。

以上、Stanをつかったデータ分析の一例でした。

Enjoy!!

 

Written on December 9, 2016