kamonohashi – ページ 8 – かものはしの分析ブログ

某洋楽ヒットチャートの週次ランキングデータをBeautiful Soupで集めてみた

はじめに

知人より、洋楽の流行りに疎いのでキャッチアップしたいという要望があり、某洋楽ヒットチャートの週次ランキングとTop100のデータを大量に集めてみようと思うに至りました。今回は深い考察を行うには至っていませんが、簡単にRにて集計・可視化を行います。

データ収集

Webスクレイピング対象の某洋楽ヒットチャートの週次ランキングは今週の順位・先週の順位・アーティスト名・曲名・詳細ページへのリンクなどが載せられおり、毎週土曜日更新されています。サイト内から導線はありませんが、URLのパラメータに法則があるため、うまく収集できます。今回は2010年8月〜2017年6月の約7年分のデータを集めます。

URLのリストをCSVで読み込み、BeautifulSoupでタグを指定して抽出します。

#ライブラリの読み込み
import urllib
from bs4 import BeautifulSoup
from urllib.request import urlopen
import requests
from urllib.error import HTTPError
import csv, time

#URLのリストを読み込む
f = open('url_list.csv', 'r')
dataReader = csv.reader(f)

#結果の出力用のリストを作る
data01 =[] #URL
data02 =[] #順位
data03 =[] #title&artist
data04 =[] #link
data05 =[] #artist

for row in dataReader:
       for url in row:
            time.sleep(10.0) #sleep(秒指定)
            try:
                    r = requests.get(url)
                    soup =  BeautifulSoup(r.content, 'html.parser')
                    
                    for body in soup.findAll("tbody"):
                        for detail in body.findAll("tr"):
                            for ranking in detail.findAll("td",{'class':'rank_td'}):
                                for content in detail.findAll("div",{'class':'name_detail'}):
                                    for artist_name in content.findAll("span"):
                                        for link in artist_name.findAll("a"):
                                            data01.append(url)
                                            data02.append(''.join(ranking.findAll(text=True)))
                                            data03.append(''.join(content.findAll(text=True)))
                                            data04.append(link.get("href"))
                                            data05.append(''.join(artist_name.findAll(text=True)))
                                            data = zip(data01,data02,data03,data04,data05)
                                            #CSV出力
                                            with open('result.csv','wt',errors='backslashreplace') as fout:
                                                writecsv = csv.writer(fout,lineterminator='\n')
                                                writecsv.writerows(data)                                   
                                
                                    
            except HTTPError as e:
                print(e.code)

#ライブラリの読み込み

import urllib

from bs4 import BeautifulSoup

from urllib.request import urlopen

import requests

from urllib.error import HTTPError

import csv, time

#URLのリストを読み込む

f = open('url_list.csv', 'r')

dataReader = csv.reader(f)

#結果の出力用のリストを作る

data01 =[] #URL

data02 =[] #順位

data03 =[] #title&artist

data04 =[] #link

data05 =[] #artist

for row in dataReader:

for url in row:

time.sleep(10.0) #sleep(秒指定)

try:

r = requests.get(url)

soup = BeautifulSoup(r.content, 'html.parser')

for body in soup.findAll("tbody"):

for detail in body.findAll("tr"):

for ranking in detail.findAll("td",{'class':'rank_td'}):

for content in detail.findAll("div",{'class':'name_detail'}):

for artist_name in content.findAll("span"):

for link in artist_name.findAll("a"):

data01.append(url)

data02.append(''.join(ranking.findAll(text=True)))

data03.append(''.join(content.findAll(text=True)))

data04.append(link.get("href"))

data05.append(''.join(artist_name.findAll(text=True)))

data = zip(data01,data02,data03,data04,data05)

#CSV出力

with open('result.csv','wt',errors='backslashreplace') as fout:

writecsv = csv.writer(fout,lineterminator='\n')

writecsv.writerows(data)

except HTTPError as e:

print(e.code)

データ取得後は簡単にpandasのstr.replaceで整形すると、以下のような結果になります。今週の順位と先週の順位が引っ付いてしまっています。

#Webスクレイピングしたデータをpandasで整形
import pandas as pd

#データ読み込み
data_list = pd.read_csv("result.csv", header = None, delimiter=",", encoding='utf-8')
data_list.columns = ['url', 'ranking','title&artist','link','artist']

#文字列の置換
data_list['url'] = data_list['url'].str.replace('置換前','置換後')
data_list['title&artist'] = data_list['title&artist'].str.replace('置換前','置換後')
data_list['ranking'] = data_list['ranking'].str.replace('置換前','置換後')

#CSVに保存
data_list.to_csv("dataset.csv",index=False)

#Webスクレイピングしたデータをpandasで整形

import pandas as pd

#データ読み込み

data_list = pd.read_csv("result.csv", header = None, delimiter=",", encoding='utf-8')

data_list.columns = ['url', 'ranking','title&artist','link','artist']

#文字列の置換

data_list['url'] = data_list['url'].str.replace('置換前','置換後')

data_list['title&artist'] = data_list['title&artist'].str.replace('置換前','置換後')

data_list['ranking'] = data_list['ranking'].str.replace('置換前','置換後')

#CSVに保存

data_list.to_csv("dataset.csv",index=False)

ここから横着してRで整形し、各週の順位データを作成しました。

#引っ付いたランキングに関するデータを分離して、データテーブルに加える。
dataset <- strsplit(ranking_dataset$ranking, "from")
rankings <- data.frame(current_ranking=as.integer(),
                       previous_ranking=as.integer())
for( i in 1:nrow(ranking_dataset)){
  rankings[i,1] <- as.integer(dataset[[i]][1])
  rankings[i,2] <- as.integer(dataset[[i]][2])
}

ranking_dataset <- cbind(ranking_dataset,rankings)

ranking_dataset <- ranking_dataset %>% select(url,current_ranking,previous_ranking,title.artist,artist,link)
colnames(ranking_dataset) <- c("date","current_ranking","previous_ranking","title",
                               "artist","link")

#引っ付いたランキングに関するデータを分離して、データテーブルに加える。

dataset <- strsplit(ranking_dataset$ranking, "from")

rankings <- data.frame(current_ranking=as.integer(),

previous_ranking=as.integer())

for( i in 1:nrow(ranking_dataset)){

rankings[i,1] <- as.integer(dataset[[i]][1])

rankings[i,2] <- as.integer(dataset[[i]][2])

}

ranking_dataset <- cbind(ranking_dataset,rankings)

ranking_dataset <- ranking_dataset %>% select(url,current_ranking,previous_ranking,title.artist,artist,link)

colnames(ranking_dataset) <- c("date","current_ranking","previous_ranking","title",

"artist","link")

データ確認

データ構造はこのような形です。

#データ構造を確認する
> str(ranking_dataset)
'data.frame':	22254 obs. of  6 variables:
 $ date            : chr  "2010/08/21" "2010/08/21" "2010/08/21" "2010/08/21" ...
 $ current_ranking : int  1 2 3 4 5 6 7 8 10 11 ...
 $ previous_ranking: int  1 3 NA 2 5 4 6 9 11 8 ...
 $ title           : chr  "Love The Way You Lie" "Dynamite" "Mine" "California Gurls" ...
 $ artist          : chr  "Eminem Featuring Rihanna" "Taio Cruz" "Taylor Swift" "Katy Perry Featuring Snoop Dogg" ...
 $ link            : chr  "/artists/detail/306903" "/artists/detail/464754" "/artists/detail/319433" "/artists/detail/450934" ...

#データ構造を確認する

> str(ranking_dataset)

'data.frame': 22254 obs. of 6 variables:

$ date : chr "2010/08/21" "2010/08/21" "2010/08/21" "2010/08/21" ...

$ current_ranking : int 1 2 3 4 5 6 7 8 10 11 ...

$ previous_ranking: int 1 3 NA 2 5 4 6 9 11 8 ...

$ title : chr "Love The Way You Lie" "Dynamite" "Mine" "California Gurls" ...

$ artist : chr "Eminem Featuring Rihanna" "Taio Cruz" "Taylor Swift" "Katy Perry Featuring Snoop Dogg" ...

$ link : chr "/artists/detail/306903" "/artists/detail/464754" "/artists/detail/319433" "/artists/detail/450934" ...

まず、どんな楽曲やアーティストがランキングに入っているのかを簡単に確認してみます。ほとんど聞いたことない人の名前、曲名ですが。。

#楽曲ごとのランクイン回数の上位20人
> ranking_dataset %>% count(title) %>% arrange(desc(n)) %>% print(n=20)
# A tibble: 1,785 × 2
                          title     n
                          <chr> <int>
1                   Radioactive    91
2                          Sail    78
3             Party Rock Anthem    68
4                Counting Stars    67
5                       Animals    66
6           Rolling In The Deep    65
7                        Ho Hey    62
8                         Sorry    62
9  Somebody That I Used To Know    61
10                       Demons    60
11                    All Of Me    58
12                       Lights    57
13 Dark Horse Featuring Juicy J    56
14                  Some Nights    56
15                 Stay With Me    55
16                        Hello    54
17      Can't Stop The Feeling!    52
18     Cheap Thrills Featuring     52
19            Don't Let Me Down    52
20                      Pompeii    52
# ... with 1,765 more rows

#楽曲ごとのランクイン回数の上位20人

> ranking_dataset %>% count(title) %>% arrange(desc(n)) %>% print(n=20)

# A tibble: 1,785 × 2

title n

1 Radioactive 91

2 Sail 78

3 Party Rock Anthem 68

4 Counting Stars 67

5 Animals 66

6 Rolling In The Deep 65

7 Ho Hey 62

8 Sorry 62

9 Somebody That I Used To Know 61

10 Demons 60

11 All Of Me 58

12 Lights 57

13 Dark Horse Featuring Juicy J 56

14 Some Nights 56

15 Stay With Me 55

16 Hello 54

17 Can't Stop The Feeling! 52

18 Cheap Thrills Featuring 52

19 Don't Let Me Down 52

20 Pompeii 52

# ... with 1,765 more rows

#アーティストごとのランクイン回数の上位20位
> ranking_dataset %>% count(artist) %>% arrange(desc(n)) %>% print(n=20)
# A tibble: 876 × 2
             artist     n
              <chr> <int>
1             Drake   350
2        Bruno Mars   325
3           Rihanna   301
4      Taylor Swift   299
5             Adele   294
6   Imagine Dragons   258
7     One Direction   254
8        Katy Perry   237
9     Justin Bieber   204
10      Chris Brown   195
11 Carrie Underwood   192
12  Lady Antebellum   184
13      Nicki Minaj   183
14      OneRepublic   183
15   Ellie Goulding   182
16         Maroon 5   182
17       The Weeknd   182
18       Ed Sheeran   169
19    Blake Shelton   167
20          Beyonce   152
# ... with 856 more rows

#アーティストごとのランクイン回数の上位20位

> ranking_dataset %>% count(artist) %>% arrange(desc(n)) %>% print(n=20)

# A tibble: 876 × 2

artist n

1 Drake 350

2 Bruno Mars 325

3 Rihanna 301

4 Taylor Swift 299

5 Adele 294

6 Imagine Dragons 258

7 One Direction 254

8 Katy Perry 237

9 Justin Bieber 204

10 Chris Brown 195

11 Carrie Underwood 192

12 Lady Antebellum 184

13 Nicki Minaj 183

14 OneRepublic 183

15 Ellie Goulding 182

16 Maroon 5 182

17 The Weeknd 182

18 Ed Sheeran 169

19 Blake Shelton 167

20 Beyonce 152

# ... with 856 more rows

続いて、2010年8月〜2017年6月の間に100位以内に入った数を楽曲ごとにヒストグラムにしてみます。べき乗分布な形かと思いきや、20回前後で盛り上がっているのが気になりますね。

#タイトルごとのランクイン回数の集計値をヒストグラムに示す
title_ranking <- ranking_dataset %>% count(title) %>% arrange(desc(n))
ggplot(data = title_ranking,aes(x = n)) + geom_histogram() + xlab("the number of Appearance of Title")

#タイトルごとのランクイン回数の集計値をヒストグラムに示す

title_ranking <- ranking_dataset %>% count(title) %>% arrange(desc(n))

ggplot(data = title_ranking,aes(x = n)) + geom_histogram() + xlab("the number of Appearance of Title")

#楽曲のランクイン回数のサマリー
> summary(title_ranking$n)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   1.00    1.00    9.00   12.47   20.00   91.00

#楽曲のランクイン回数のサマリー

> summary(title_ranking$n)

Min. 1st Qu. Median Mean 3rd Qu. Max.

1.00 1.00 9.00 12.47 20.00 91.00

中央値が9週間なので、意外と長い期間Top100には入っています。

続いて、100位以内に入った数をアーティストごとにヒストグラムにしてみます。こちらはべき乗分布のような形になっています。

#アーティストごとのランクイン回数の集計値をヒストグラムに示す
artist_ranking <- ranking_dataset %>% count(artist) %>% arrange(desc(n))
ggplot(data = artist_ranking,aes(x = n)) + geom_histogram() + xlab("the number of Appearance of Artist")

#アーティストごとのランクイン回数の集計値をヒストグラムに示す

artist_ranking <- ranking_dataset %>% count(artist) %>% arrange(desc(n))

ggplot(data = artist_ranking,aes(x = n)) + geom_histogram() + xlab("the number of Appearance of Artist")

Top10入りの楽曲の実態

Top10に入っている楽曲のみに絞って、ヒストグラムを描いてみます。

#Top10入りのデータに絞ってヒストグラムを描く
top_10 <- ranking_dataset %>% filter(current_ranking<=10) %>% count(title) %>% arrange(desc(n))
ggplot(data = top_10,aes(x = n)) + geom_histogram() + xlab("the number of Appearance of Title in Top10")

#Top10入りのデータに絞ってヒストグラムを描く

top_10 <- ranking_dataset %>% filter(current_ranking<=10) %>% count(title) %>% arrange(desc(n))

ggplot(data = top_10,aes(x = n)) + geom_histogram() + xlab("the number of Appearance of Title in Top10")

#Top10に入っている楽曲の10位以内ランクイン回数サマリー
> summary(top_10$n)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1.000   3.000  10.000   9.779  15.000  32.000

#Top10に入っている楽曲の10位以内ランクイン回数サマリー

> summary(top_10$n)

Min. 1st Qu. Median Mean 3rd Qu. Max.

1.000 3.000 10.000 9.779 15.000 32.000

Top10に入ったら、10週近くは10位以内に含まれるようです。上位はすぐに取って代わられるのかと思いきや、人気が人気を呼ぶとかなのでしょうか。確か、某ECサイトの方が、生キャラメルは売れるから売れたんだとか言っていた気がします。

順位の推移

100位以内にランクインした回数が最も多かった楽曲のTop10に関して、時系列プロットをしてみます。

#100位内ランクイン回数が上位の10曲に関して時系列プロットを行う。
title_ranking_top10 <- title_ranking[1:10,]
title_ranking_top10 <- ranking_dataset %>% filter(title %in% title_ranking_top10$title)
title_ranking_top10$date <- as.Date(title_ranking_top10$date)
ggplot(data = title_ranking_top10,aes(x = date,
                                      y = current_ranking,
                                      colour=title)) + geom_line() + 
  　　　　　　　　　　　　　　　　　　ylab("current ranking") + 
  　　　　　　　　　　　　　　　　　　ggtitle("Time Series Plot of Top10")

#100位内ランクイン回数が上位の10曲に関して時系列プロットを行う。

title_ranking_top10 <- title_ranking[1:10,]

title_ranking_top10 <- ranking_dataset %>% filter(title %in% title_ranking_top10$title)

title_ranking_top10$date <- as.Date(title_ranking_top10$date)

ggplot(data = title_ranking_top10,aes(x = date,

y = current_ranking,

colour=title)) + geom_line() +

　　　　　　　　　　　　　　　　　　ylab("current ranking") +

　　　　　　　　　　　　　　　　　　ggtitle("Time Series Plot of Top10")

初回に上位にランクインして、後は下がるだけの楽曲や、じわじわとランキングを上げていく楽曲などが観察されています。楽曲の消費のサイクルみたいなものがあるのでしょう。

今後について

せっかく面白そうなデータが手に入ったので、リンク先も辿って、どのような楽曲やアーティストの特徴が人気に繋がりうるのか見てみるのも良いですね。あと、もう少し洋楽聴いてみようと思います。私はクラシック音楽とジャズしかウォークマンに入っていないので。

参考文献

Pythonクローリング&スクレイピング -データ収集・解析のための実践開発ガイド-

ラーメン二郎の某飲食店レビューサイトデータに対して共分散構造分析をしてみる

データ

ラーメン二郎に関して、某飲食店レビューサイトのデータをWebスクレイピングしたもので、「料理・味」・「サービス」・「雰囲気」・「CP」・「酒・ドリンク」の評価項目に関して、1~5の実数値が割り振られています。ラーメン二郎の店舗数40店のうち、欠損のなかった37店舗の評価データとなります。

記述統計量・ヒストグラム

まずは記述統計量とヒストグラムを見てみます。

> describe(jirou_dataset[,-1])
                 vars  n mean   sd median trimmed  mad  min  max range  skew kurtosis   se
taste               1 37 3.58 0.09   3.59    3.59 0.01 3.10 3.73  0.63 -3.77    17.97 0.02
service             2 37 3.51 0.11   3.55    3.53 0.04 3.10 3.59  0.49 -2.64     6.41 0.02
atmosphere          3 37 3.47 0.16   3.53    3.50 0.03 3.00 3.59  0.59 -2.11     3.06 0.03
cost_performance    4 37 3.61 0.12   3.59    3.61 0.01 3.07 3.95  0.88 -1.46     9.16 0.02
drink               5 37 3.08 0.09   3.07    3.07 0.10 2.99 3.33  0.34  1.00     0.39 0.01

> describe(jirou_dataset[,-1])

vars n mean sd median trimmed mad min max range skew kurtosis se

taste 1 37 3.58 0.09 3.59 3.59 0.01 3.10 3.73 0.63 -3.77 17.97 0.02

service 2 37 3.51 0.11 3.55 3.53 0.04 3.10 3.59 0.49 -2.64 6.41 0.02

atmosphere 3 37 3.47 0.16 3.53 3.50 0.03 3.00 3.59 0.59 -2.11 3.06 0.03

cost_performance 4 37 3.61 0.12 3.59 3.61 0.01 3.07 3.95 0.88 -1.46 9.16 0.02

drink 5 37 3.08 0.09 3.07 3.07 0.10 2.99 3.33 0.34 1.00 0.39 0.01

#各評価項目に関してヒストグラムを描く
hist_jirou_dataset <- melt(jirou_dataset)
g <- ggplot(data = hist_jirou_dataset,
            aes(x = value,
                y = ..density..)) +
  geom_histogram(alpha = 0.5,position = "identity") +
  geom_density(alpha = 0)
g + facet_wrap(~variable,nrow=5)

#各評価項目に関してヒストグラムを描く

hist_jirou_dataset <- melt(jirou_dataset)

g <- ggplot(data = hist_jirou_dataset,

aes(x = value,

y = ..density..)) +

geom_histogram(alpha = 0.5,position = "identity") +

geom_density(alpha = 0)

g + facet_wrap(~variable,nrow=5)

まぁ、お酒を飲むところではないので、drinkは低いですよね。しかしながら、総じて3点台なのは驚きです。昔行ったことがあって、雰囲気は決して良くはないはずなので。

共分散構造分析

今回、検証したい仮説は「ラーメン店としての質が二郎愛につながるかどうか」です。

ラーメン店としての質に繋がりそうな評価項目
「料理・味」
「CP」
「酒・ドリンク」

二郎愛に繋がりそうな評価項目
「料理・味」
「サービス」
「雰囲気」

加えて、料理・味とCPは関係していそうな項目なので、その点もパスにおいて考慮しておきます。

以下の図は仮説のイメージです。

パスを描写する際の細かい指定は参考文献を参照されると良いと思います。
さっそく、モデルを作って実行します。

#lavaanパッケージのためにモデルを描写する。
library(lavaan)

model_jirou <- 'f1 =~ cost_performance + drink + taste
                f2 =~ taste + atmosphere + service
                f2 ~ f1
                cost_performance ~~ taste'

fit <- sem(model = model_jirou,
           data = jirou_dataset[,-1],
           estimator="MLM")

summary(fit)

#lavaanパッケージのためにモデルを描写する。

library(lavaan)

model_jirou <- 'f1 =~ cost_performance + drink + taste

f2 =~ taste + atmosphere + service

f2 ~ f1

cost_performance ~~ taste'

fit <- sem(model = model_jirou,

data = jirou_dataset[,-1],

estimator="MLM")

summary(fit)

こちらが推定結果です。

> summary(fit)
lavaan (0.5-23.1097) converged normally after 100 iterations

  Number of observations                            37

  Estimator                                         ML      Robust
  Minimum Function Test Statistic                2.057       1.965
  Degrees of freedom                                 2           2
  P-value (Chi-square)                           0.358       0.374
  Scaling correction factor                                  1.047
    for the Satorra-Bentler correction

Parameter Estimates:

  Information                                 Expected
  Standard Errors                           Robust.sem

Latent Variables:
                   Estimate  Std.Err  z-value  P(>|z|)
  f1 =~                                               
    cost_performnc    1.000                           
    drink            -0.153    0.134   -1.139    0.255
    taste             0.381    0.194    1.958    0.050
  f2 =~                                               
    taste             1.000                           
    atmosphere        4.008    1.972    2.032    0.042
    service           3.147    1.423    2.211    0.027

Regressions:
                   Estimate  Std.Err  z-value  P(>|z|)
  f2 ~                                                
    f1                0.095    0.130    0.731    0.465

Covariances:
                      Estimate  Std.Err  z-value  P(>|z|)
 .cost_performance ~~                                    
   .taste               -0.001    0.007   -0.115    0.908

Intercepts:
                   Estimate  Std.Err  z-value  P(>|z|)
   .cost_performnc    3.606    0.020  179.760    0.000
   .drink             3.081    0.014  214.293    0.000
   .taste             3.582    0.015  237.571    0.000
   .atmosphere        3.471    0.027  129.316    0.000
   .service           3.509    0.019  187.760    0.000
    f1                0.000                           
   .f2                0.000                           

Variances:
                   Estimate  Std.Err  z-value  P(>|z|)
   .cost_performnc   -0.007    0.017   -0.390    0.696
   .drink             0.007    0.002    3.170    0.002
   .taste             0.002    0.003    0.712    0.477
   .atmosphere        0.005    0.004    1.273    0.203
   .service          -0.000    0.001   -0.329    0.742
    f1                0.021    0.017    1.242    0.214
   .f2                0.001    0.001    1.282    0.200

> summary(fit)

lavaan (0.5-23.1097) converged normally after 100 iterations

Number of observations 37

Estimator ML Robust

Minimum Function Test Statistic 2.057 1.965

Degrees of freedom 2 2

P-value (Chi-square) 0.358 0.374

Scaling correction factor 1.047

for the Satorra-Bentler correction

Parameter Estimates:

Information Expected

Standard Errors Robust.sem

Latent Variables:

Estimate Std.Err z-value P(>|z|)

f1 =~

cost_performnc 1.000

drink -0.153 0.134 -1.139 0.255

taste 0.381 0.194 1.958 0.050

f2 =~

taste 1.000

atmosphere 4.008 1.972 2.032 0.042

service 3.147 1.423 2.211 0.027

Regressions:

Estimate Std.Err z-value P(>|z|)

f2 ~

f1 0.095 0.130 0.731 0.465

Covariances:

Estimate Std.Err z-value P(>|z|)

.cost_performance ~~

.taste -0.001 0.007 -0.115 0.908

Intercepts:

Estimate Std.Err z-value P(>|z|)

.cost_performnc 3.606 0.020 179.760 0.000

.drink 3.081 0.014 214.293 0.000

.taste 3.582 0.015 237.571 0.000

.atmosphere 3.471 0.027 129.316 0.000

.service 3.509 0.019 187.760 0.000

f1 0.000

.f2 0.000

Variances:

Estimate Std.Err z-value P(>|z|)

.cost_performnc -0.007 0.017 -0.390 0.696

.drink 0.007 0.002 3.170 0.002

.taste 0.002 0.003 0.712 0.477

.atmosphere 0.005 0.004 1.273 0.203

.service -0.000 0.001 -0.329 0.742

f1 0.021 0.017 1.242 0.214

.f2 0.001 0.001 1.282 0.200

lavaan (0.5-23.1097) converged normally after 100 iterationsとあるので、適切に推定されたようです。自由度が2とデータ数が少ないためかなりギリギリな推定となっています。ただ、推定すべき母数の数よりもデータ数が一応多い状態ではあります。二郎インスパイア系の店のデータも集めた方がいいかもしれないですね。

モデルの評価

モデルの評価として適合度と母数の推定に関して見ていきます。

適合度
・適合度指標であるCFI(Comparative Fit Index)が1なので、適合度に関しては良さそうです。
・同じく適合度指標であるTLI(Tucker-Lewis Index)が0.998なので、適合度に関しては良さそうです。
・0.05以下であれば当てはまりが良いとされるRMSEAは0.028なので、当てはまりは良さそうです。
・0に近いほど良いとされるSRMRは0.024となっています。

母数の推定

f2 ~ f1はラーメン店としての質と二郎愛の関係を想定したものですが、p値が0.47と全然だめでした。ラーメン店としての質が二郎愛に繋がるという仮説は正しいとは言えないです。

参考文献にもあるよう、係数の解釈を行いやすくするために標準化推定値を求めます。

> standardizedsolution(fit)
                lhs op              rhs est.std     se      z pvalue
1                f1 =~ cost_performance   1.208  0.524  2.303  0.021
2                f1 =~            drink  -0.257  0.178 -1.448  0.148
3                f1 =~            taste   0.612  0.325  1.883  0.060
4                f2 =~            taste   0.400  0.120  3.327  0.001
5                f2 =~       atmosphere   0.901  0.066 13.723  0.000
6                f2 =~          service   1.016  0.049 20.694  0.000
7                f2  ~               f1   0.382  0.338  1.132  0.258
8  cost_performance ~~            taste  -0.201  1.651 -0.122  0.903
9  cost_performance ~~ cost_performance  -0.459  1.267 -0.362  0.717
10            drink ~~            drink   0.934  0.091 10.226  0.000
11            taste ~~            taste   0.278  0.399  0.698  0.485
12       atmosphere ~~       atmosphere   0.187  0.118  1.582  0.114
13          service ~~          service  -0.033  0.100 -0.333  0.740
14               f1 ~~               f1   1.000  0.000     NA     NA
15               f2 ~~               f2   0.854  0.258  3.309  0.001
16 cost_performance ~1                   29.960  8.666  3.457  0.001
17            drink ~1                   35.715  4.676  7.638  0.000
18            taste ~1                   39.595 15.318  2.585  0.010
19       atmosphere ~1                   21.553  4.333  4.974  0.000
20          service ~1                   31.293  7.949  3.937  0.000
21               f1 ~1                    0.000  0.000     NA     NA
22               f2 ~1                    0.000  0.000     NA     NA

> standardizedsolution(fit)

lhs op rhs est.std se z pvalue

1 f1 =~ cost_performance 1.208 0.524 2.303 0.021

2 f1 =~ drink -0.257 0.178 -1.448 0.148

3 f1 =~ taste 0.612 0.325 1.883 0.060

4 f2 =~ taste 0.400 0.120 3.327 0.001

5 f2 =~ atmosphere 0.901 0.066 13.723 0.000

6 f2 =~ service 1.016 0.049 20.694 0.000

7 f2 ~ f1 0.382 0.338 1.132 0.258

8 cost_performance ~~ taste -0.201 1.651 -0.122 0.903

9 cost_performance ~~ cost_performance -0.459 1.267 -0.362 0.717

10 drink ~~ drink 0.934 0.091 10.226 0.000

11 taste ~~ taste 0.278 0.399 0.698 0.485

12 atmosphere ~~ atmosphere 0.187 0.118 1.582 0.114

13 service ~~ service -0.033 0.100 -0.333 0.740

14 f1 ~~ f1 1.000 0.000 NA NA

15 f2 ~~ f2 0.854 0.258 3.309 0.001

16 cost_performance ~1 29.960 8.666 3.457 0.001

17 drink ~1 35.715 4.676 7.638 0.000

18 taste ~1 39.595 15.318 2.585 0.010

19 atmosphere ~1 21.553 4.333 4.974 0.000

20 service ~1 31.293 7.949 3.937 0.000

21 f1 ~1 0.000 0.000 NA NA

22 f2 ~1 0.000 0.000 NA NA

これによると、f1(ラーメン店としての質)に対してはコストパフォーマンスが最も影響を与えるようです。味よりもコストパフォーマンスが勝っているという考察になりますが、それはそれで面白いですね。他方、f2(二郎愛)に対しては味よりも雰囲気・サービスが影響を与えるようです。

semPlotパッケージを用いてパス図を出力してみます。

#パス図の出力
library(semPlot)

semPaths(fit,whatLabels = "stand",optimizeLatRes = T)

#パス図の出力

library(semPlot)

semPaths(fit,whatLabels = "stand",optimizeLatRes = T)

変数名が3文字に省略されているのですが、修正方法がパッと見つからなかったので、そのまま載せています。

初歩の初歩ですが、一通りの進め方がわかったので、今後も共分散構造分析にチャレンジしてみたいと思います。

参考文献

共分散構造分析 R編―構造方程式モデリング

RStanで学部時代の研究を振り返ってみる

研究概要

大学時代に実験経済学で行った実験結果のデータがUSBに入っていたので、振り返って分析などをしてみたいと思います。

研究目的
　ピア効果に関して、競争相手が自分よりも秀でた人がいいのか劣った人がいいのかを確かめる。

実験方法
・1分間で100マス計算を2セット解いてもらう。（めちゃ速い人には3枚目も渡した）
・実験開始後、実験対象のクラスによって、途中で「平均的なクラスは○○マスまで進んでいます！」とアナウンスします。アナウンスすることで、競争相手のレベルを知り、焦るなり余裕を感じるなりしてもらおうという計画です。
なお、対照群はアナウンスをしていません。アナウンス内容は「平均告知（１８秒）、上告知（１５秒）、超上告知（１２秒）、下告知（２０秒）」と4パターンとなります。
・計算が間違っているものは加点しません。

実験対象
　某国立大学の経済学部生の1~2年の必修科目履修者217名（先生に交渉して授業の開始5分を頂いて実験を行いました。）
　内部進学やスポ専などがない分、計算能力的にある程度近い集団ではないかと思われます。

検証方法
　アナウンスごとに100マス計算の点数の水準が変わりうるのかを回帰分析などで判断。

データ可視化

以下、実験カテゴリごとの略記です。
下告知（２０秒）・・・slow_20
上告知（１５秒）・・・fast_15
超上告知（１２秒）・・・fastest_12
平均告知（１８秒）・・・average_18
対照群・・・baseline

データ構造の確認です。

> str(dataset)
'data.frame':	217 obs. of  4 variables:
 $ categories    : chr  "baseline" "baseline" "baseline" "baseline" ...
 $ points        : int  72 79 81 81 98 99 100 101 102 104 ...
 $ errors        : int  0 0 0 0 2 1 0 0 1 0 ...
 $ genuine_points: int  72 79 81 81 96 98 100 101 101 104 ...

> str(dataset)

'data.frame': 217 obs. of 4 variables:

$ categories : chr "baseline" "baseline" "baseline" "baseline" ...

$ points : int 72 79 81 81 98 99 100 101 102 104 ...

$ errors : int 0 0 0 0 2 1 0 0 1 0 ...

$ genuine_points: int 72 79 81 81 96 98 100 101 101 104 ...

平均値、中央値、標準偏差、サンプルサイズを出してみます。

> dataset %>% group_by(categories) %>% summarise(average=mean(genuine_points),
+                                                median=median(genuine_points),
+                                                stdev=sd(genuine_points),
+                                                sample=n())
# A tibble: 5 × 5
  categories  average median    stdev sample
       <chr>    <dbl>  <dbl>    <dbl>  <int>
1 average_18 120.8605    117 28.08455     43
2   baseline 121.0000    121 22.19109     46
3    fast_15 123.0357    126 24.32355     56
4 fastest_12 123.6774    125 24.09756     31
5    slow_20 126.3902    126 23.20547     41

> dataset %>% group_by(categories) %>% summarise(average=mean(genuine_points),

+ median=median(genuine_points),

+ stdev=sd(genuine_points),

+ sample=n())

# A tibble: 5 × 5

categories average median stdev sample

1 average_18 120.8605 117 28.08455 43

2 baseline 121.0000 121 22.19109 46

3 fast_15 123.0357 126 24.32355 56

4 fastest_12 123.6774 125 24.09756 31

5 slow_20 126.3902 126 23.20547 41

中央値で見てみると、baselineに対してわずかですが点数に違いがありそうに見えます。

実験種別で点数に関するヒストグラムと確率密度関数を確認してみます。

library(ggplot2)

g <- ggplot(data = dataset,
            aes(x = genuine_points,
                y = ..density..)) +
            geom_histogram(alpha = 0.5,position = "identity") +
            geom_density(alpha = 0)
g + facet_wrap(~categories,nrow=5)

library(ggplot2)

g <- ggplot(data = dataset,

aes(x = genuine_points,

y = ..density..)) +

geom_histogram(alpha = 0.5,position = "identity") +

geom_density(alpha = 0)

g + facet_wrap(~categories,nrow=5)

baselineが多峰性がありそうなのが気になります。average_18は低そうに見えますね。

RStanで重回帰

『StanとRでベイズ統計モデリング』にあるコードを参考にしています。正規分布を事前分布にした線形回帰モデルです。
被説明変数が点数、説明変数が実験種別のダミー変数だけからなります。

library(rstan)
library(dummies)

dummies <- dummy.data.frame(dataset, sep = "_", names = c("categories"))

analytical_dataset <- dummies %>% select(categories_average_18,
                                         categories_fast_15,
                                         categories_fastest_12,
                                         categories_slow_20,
                                         genuine_points)


data <- list(N=nrow(analytical_dataset),
             average_18=analytical_dataset$categories_average_18,
             fast_15=analytical_dataset$categories_fast_15,
             fastest_12=analytical_dataset$categories_fastest_12,
             slow_20=analytical_dataset$categories_slow_20,
             genuine_points=analytical_dataset$genuine_points)


stan_code <- "
data{
int N; //the number of student
int<lower=0> genuine_points[N];
real<lower=0, upper=1> average_18[N];
real<lower=0, upper=1> fast_15[N];
real<lower=0, upper=1> fastest_12[N];
real<lower=0, upper=1> slow_20[N];
}

parameters{
real b1;
real b2;
real b3;
real b4;
real b5;
real<lower=0> sigma;
}

transformed parameters{
real mu[N];
for(n in 1:N)
mu[n] = b1 + b2*average_18[n] + b3*fast_15[n] + b4*fastest_12[n] + b5*slow_20[n];
}

model{
for(n in 1:N)
genuine_points[n] ~ normal(mu[n], sigma);
}
"

fit <- stan(model_code =stan_code, data=data, seed=1234)
fit.summary <-data.frame(summary(fit)$summary)
head(fit.summary,6)

library(rstan)

library(dummies)

dummies <- dummy.data.frame(dataset, sep = "_", names = c("categories"))

analytical_dataset <- dummies %>% select(categories_average_18,

categories_fast_15,

categories_fastest_12,

categories_slow_20,

genuine_points)

data <- list(N=nrow(analytical_dataset),

average_18=analytical_dataset$categories_average_18,

fast_15=analytical_dataset$categories_fast_15,

fastest_12=analytical_dataset$categories_fastest_12,

slow_20=analytical_dataset$categories_slow_20,

genuine_points=analytical_dataset$genuine_points)

stan_code <- "

data{

int N; //the number of student

int<lower=0> genuine_points[N];

real<lower=0, upper=1> average_18[N];

real<lower=0, upper=1> fast_15[N];

real<lower=0, upper=1> fastest_12[N];

real<lower=0, upper=1> slow_20[N];

}

parameters{

real b1;

real b2;

real b3;

real b4;

real b5;

real<lower=0> sigma;

}

transformed parameters{

real mu[N];

for(n in 1:N)

mu[n] = b1 + b2*average_18[n] + b3*fast_15[n] + b4*fastest_12[n] + b5*slow_20[n];

}

model{

for(n in 1:N)

genuine_points[n] ~ normal(mu[n], sigma);

}

fit <- stan(model_code =stan_code, data=data, seed=1234)

fit.summary <-data.frame(summary(fit)$summary)

head(fit.summary,6)

結果

traceplot(fit)でMCMCのサンプリング結果を確認する。

収束しているように見えます。

以下推定結果ですが、残念ながらベイズ予測区間において符号の逆転が起きていないものはなかったので、アナウンスによる効果があるとは言えないようです。ただ、slow_20の係数がおしいですね。少なくとも他の実験種別よりも、アナウンス効果があるかもしれないという考察に止まりそうです。

> head(fit.summary,6)
             mean   se_mean       sd      X2.5.       X25.        X50.       X75.     X97.5.    n_eff      Rhat
b1    121.0253094 0.1034836 3.460792 114.214104 118.793988 121.0140383 123.213054 127.902759 1118.429 0.9998152
b2     -0.1195318 0.1326887 4.957153 -10.009827  -3.372795  -0.1304577   3.117195   9.694184 1395.715 1.0004034
b3      1.9440496 0.1248433 4.661571  -7.494857  -1.043060   1.9756929   5.007503  10.827640 1394.228 1.0000830
b4      2.5460694 0.1436941 5.447930  -8.026616  -1.128587   2.4997953   6.220891  13.044467 1437.426 1.0008688
b5      5.3029965 0.1351662 4.988323  -4.479865   2.004280   5.3108476   8.543292  15.458379 1361.988 1.0004238
sigma  24.6158717 0.0211252 1.150481  22.501612  23.811858  24.5900106  25.368057  26.970197 2965.905 1.0017414

> head(fit.summary,6)

mean se_mean sd X2.5. X25. X50. X75. X97.5. n_eff Rhat

b1 121.0253094 0.1034836 3.460792 114.214104 118.793988 121.0140383 123.213054 127.902759 1118.429 0.9998152

b2 -0.1195318 0.1326887 4.957153 -10.009827 -3.372795 -0.1304577 3.117195 9.694184 1395.715 1.0004034

b3 1.9440496 0.1248433 4.661571 -7.494857 -1.043060 1.9756929 5.007503 10.827640 1394.228 1.0000830

b4 2.5460694 0.1436941 5.447930 -8.026616 -1.128587 2.4997953 6.220891 13.044467 1437.426 1.0008688

b5 5.3029965 0.1351662 4.988323 -4.479865 2.004280 5.3108476 8.543292 15.458379 1361.988 1.0004238

sigma 24.6158717 0.0211252 1.150481 22.501612 23.811858 24.5900106 25.368057 26.970197 2965.905 1.0017414

分布でも見てみます。アナウンス効果が0を確実に超えているとは言えないですね。

library (reshape)
library(dplyr)
library(ggplot2)
post   <- extract (fit, permuted = F)
m.post <- melt (post)
m.post <- m.post %>% filter(parameters %in% c("b1","b2","b3","b4","b5"))
graph  <- ggplot (m.post, aes(x = value))
graph  <- graph + geom_density () + facet_grid(. ~ parameters, scales = "free") + theme_bw() 
plot (graph)

library (reshape)

library(dplyr)

library(ggplot2)

post <- extract (fit, permuted = F)

m.post <- melt (post)

m.post <- m.post %>% filter(parameters %in% c("b1","b2","b3","b4","b5"))

graph <- ggplot (m.post, aes(x = value))

graph <- graph + geom_density () + facet_grid(. ~ parameters, scales = "free") + theme_bw()

plot (graph)

結局、学部時代のレポートと結論は変わらないのですが、係数が0よりも大きい確率という観点で結果に向き合えたのは良かったと思います。

参考文献

StanとRでベイズ統計モデリング (Wonderful R)

データマイニングに関するSlideShareを大量に集めてみた#2

前回の「SlideShareの機械学習に関するスライドを大量に集めてみた」でSlideShareの検索機能に不満のある方が、やはりおられたようなので、他のバージョンも作っていきます。今回はGoogleのsiteコマンド検索（　「データマイニング site:slideshare.net　」）でヒットした520件を対象にWebスクレイピングし、タイトルで検索できるようになっています。前回からの改良点は、URLがリンクになっている点です。（スマホユーザーも嬉しい）
ビュー数順に並んでいますが、TablePressのビュー数でソートするとどうやら整数と認識されていないようで、うまく並び替えができないようです。

たまに企業での分析事例共有などのスライドもあるので、仕事のアイデアも手に入ったりします。書籍になっていなかったりするし、Google検索だけで見つけようとすると大変なので、この試みは当分続けていこうと思います。本文内のテキストを取得してタグ付けしておくのも行う予定です。

[table id=4 /]

SlideShareの機械学習に関するスライドを大量に集めてみた

SlideShareの検索機能に満足できなかったので、Googleのsiteコマンド検索（　「機械学習 site:slideshare.net　」）で集めてき
た、機械学習に関するスライド520件のURLに対して、タイトル・ビュー数・公開日をWebスクレイピングにより取得してみました。

以下の表はビュー数の順番で並んでいます。WordPressのTablepressというプラグインを使っているので、簡単に検索ないし100件表示なども行えます。数の多いものから見ていこうと思います。

(取得用の稚拙なPythonコードは最下部に載っています。)
(URLがリンクになるように改良しました。)

[table id=3 /]

一応、SlideShareのURLのCSVファイル（N行1列のデータ）があれば取得できるPythonコード（2系）を載せておきます。

import urllib, urllib2, csv
from bs4 import BeautifulSoup
from BeautifulSoup import BeautifulSoup
from urllib2 import Request, urlopen, URLError, HTTPError
import sys, string, codecs, time, re
from httplib import BadStatusLine, IncompleteRead

#CSVでリストを読み込む。
f = open('url_list_ss.csv', 'rb')
dataReader = csv.reader(f)

#結果の出力用のリストを作る。
data01 =[]
data02 =[]
data03 =[]
data04 =[]
data05 =[]

for row in dataReader:
       for url in row:
              try:
                     res = urllib2.urlopen(url)
                     soup = BeautifulSoup(res.read(),fromEncoding="utf-8")
                     time.sleep(5.0) #sleep(秒指定)
                     
                     for title in soup.findAll("div",{"class":"small-10 columns"}):
                            for views in soup.findAll("div",{"class":"small-2 columns text-right format-views"}):
                                for times in soup.findAll("time",{"itemprop":"datePublished"}):
                                            data01.append(url)
                                            data02.append(''.join(title.findAll(text=True)))
                                            data03.append(''.join(views.findAll(text=True)))
                                            data04.append(''.join(times.findAll(text=True)))
                                            data05.append("contentslist")
                                            
                                            #データ結合
                                            data = zip(data01,data02,data03,data04,data05)

                                            #CSV出力
                                            f= open('ss_extract_result.csv','w')
                                            writecsv = csv.writer(f,lineterminator='\n')
                                            for x in range(len(data)):
                                                writecsv.writerows([data[x]])
                                                                        

              except HTTPError, e:
                     print e.code
                     
              except BadStatusLine:
                     print "could not fetch"
                     
              except IncompleteRead:
                     print "IncompleteRead"
                     
              except IndexError:
                     print "IndexError"

import urllib, urllib2, csv

from bs4 import BeautifulSoup

from BeautifulSoup import BeautifulSoup

from urllib2 import Request, urlopen, URLError, HTTPError

import sys, string, codecs, time, re

from httplib import BadStatusLine, IncompleteRead

#CSVでリストを読み込む。

f = open('url_list_ss.csv', 'rb')

dataReader = csv.reader(f)

#結果の出力用のリストを作る。

data01 =[]

data02 =[]

data03 =[]

data04 =[]

data05 =[]

for row in dataReader:

for url in row:

try:

res = urllib2.urlopen(url)

soup = BeautifulSoup(res.read(),fromEncoding="utf-8")

time.sleep(5.0) #sleep(秒指定)

for title in soup.findAll("div",{"class":"small-10 columns"}):

for views in soup.findAll("div",{"class":"small-2 columns text-right format-views"}):

for times in soup.findAll("time",{"itemprop":"datePublished"}):

data01.append(url)

data02.append(''.join(title.findAll(text=True)))

data03.append(''.join(views.findAll(text=True)))

data04.append(''.join(times.findAll(text=True)))

data05.append("contentslist")

#データ結合

data = zip(data01,data02,data03,data04,data05)

#CSV出力

f= open('ss_extract_result.csv','w')

writecsv = csv.writer(f,lineterminator='\n')

for x in range(len(data)):

writecsv.writerows([data[x]])

except HTTPError, e:

print e.code

except BadStatusLine:

print "could not fetch"

except IncompleteRead:

print "IncompleteRead"

except IndexError:

print "IndexError"

言語処理学会第23回年次大会〜Webマーケティングで使えそうな発表資料まとめ〜

行きたいけど行けなかった言語処理学会第23回年次大会の発表内容がPDFで見れるということで、発表内容の中でWebマーケティングなどの仕事で役に立つかもしれない12件の研究を独断と偏見でまとめています。

プログラムはこちらのリンクから見れます。
言語処理学会第23回年次大会(NLP2017) プログラム

今回取り上げるのは以下の12件です。

・Wikipediaのカテゴリ構造を特徴ベクトルに用いたRandom Forestによるショートメッセージ分類
・NMFを用いた為替ドル円レートの変動要因分析とAdaboostを用いた予測システム
・Are Deep Learning Methods Better for Twitter Sentiment Analysis?
・化粧品レビューサイトにおけるクチコミの有用性判定
・機械学習を用いたQAサイト質問文のカテゴリの類推
・不満調査データセットの素性ベクトル化
・検索エンジンを用いた記載欠落箇所の補完
・生コーパスからの単語難易度関連指標の予測
・常識から外れた雑談応答の検出
・店舗レビューには何が書かれているか？–調査及びその自動分類–
・トップダウン型共起グラフを用いたブログからの観光地の行動分析
・社内データに基づくイノベータ人財のピックアップ

Wikipediaのカテゴリ構造を特徴ベクトルに用いたRandom Forestによるショートメッセージ分類

http://www.anlp.jp/proceedings/annual_meeting/2017/pdf_dir/P16-4.pdf

Wikipediaのカテゴリ構造を学習させたナイーブベイズにより特徴ベクトルを生成し、それをTweetに付与して、ランダムフォレストを用いてメッセージ分類をしているようです。直接Tweetのデータを使わないでWikipediaのデータを活用するところにトリックがあるようです。ビジネスにおいて関心のある対象ユーザーのTweetのモニタリングや対象ユーザーの特定に使えるかもしれません。

NMFを用いた為替ドル円レートの変動要因分析とAdaboostを用いた予測システム

http://www.anlp.jp/proceedings/annual_meeting/2017/pdf_dir/P15-5.pdf

『金融経済月報』や『日本経済展望』のテキストデータに対して、PCAとNMFの2つの手法を用いて特徴抽出を行い、Adaboostを用いて為替レートの予測を行っているようです。為替レートを予測せずとも、テキストデータを用いた市場動向の予測などの参考になるかもしれません。

Are Deep Learning Methods Better for Twitter Sentiment Analysis?

http://www.anlp.jp/proceedings/annual_meeting/2017/pdf_dir/C5-1.pdf

Tweetの感情分析で深層学習を用いて、SVMなどと比較している研究のようです。Conclusionを見る限りは、SVMよりも精度が高いとは言えず（ネットワーク構造やデータセットの大きさに依存する）、少ないデータセットなどではSVMの方が依然として精度が高いようです。

化粧品レビューサイトにおけるクチコミの有用性判定

http://www.anlp.jp/proceedings/annual_meeting/2017/pdf_dir/C5-4.pdf

@cosmeにおける、Likeの数をクチコミの有用度と定義して、クチコミデータに対してTF-IDFやLDAを用いて構造解析や意味解析などで特徴を抽出し、その特徴を説明変数とした回帰モデル（SVR）を構築しています。構築した回帰モデルを用い、クチコミのLikeの数を推定することを目的としています。

ここでの知見が直接に自社のコンテンツなどのナレッジになるわけではないですが、LDAやSVM（SVR）を駆使しているなど、分析手法の参考にはなると思います。

機械学習を用いたQAサイト質問文のカテゴリの類推

http://www.anlp.jp/proceedings/annual_meeting/2017/pdf_dir/P9-6.pdf

深層学習（Stacked Denoising Autoencoders(SdA) と Deep Belief　Network(DBN)）を用いて質問文に適したカテゴリを類推する研究です。従来手法のMLPやSVMと比べても精度が高かったようです。ユーザーの欲しい商品や、働きたい職場、住みたい場所などに関しても類推できるとマーケティングなどでの活用もできるかもしれません。

不満調査データセットの素性ベクトル化

http://www.anlp.jp/proceedings/annual_meeting/2017/pdf_dir/P9-2.pdf

不満に関するデータにタグ付けされた情報やJUMANの付属辞書から得られた意味情報も素性とすることで、不満の検索のしやすさを高める研究のようです。不満のテキストデータから特徴量に変えるまでのプロセスが記されているのが参考になりそうです。研究自体は、形態素解析の精度などにより、まだ課題があるようです。

検索エンジンを用いた記載欠落箇所の補完

http://www.anlp.jp/proceedings/annual_meeting/2017/pdf_dir/P16-5.pdf
「検索エンジンを用いて記載欠落箇所に適切な情報を埋める研究」とあり、コンテンツマーケティングなどの際に、作成した記事に説明不足な箇所を見つけることができるので、コンテンツの質を担保する上で役に立つかもしれません。

生コーパスからの単語難易度関連指標の予測

http://www.anlp.jp/proceedings/annual_meeting/2017/pdf_dir/E5-3.pdf

「均衡コーパスや深いアノテーションを施した言語資源に頼らないことで他言語への適用が容易な単語難易度関連指標の予測手法を提案」
「単語難易度関連指標の予測精度はLDAの各トピックからの単語出現確率を基にした素性を用いることで、従来のようにコーパスからの単語頻度を素性を用いる場合と比べ大幅に向上」
「word2vecのような単語のベクトル表現より、LDAの各トピックからの単語出現確率の方が予測精度の向上に有効」
とあり、文章のリーダビリティをよりカジュアルに計算することが可能になるようです。文章のリーダビリティとWebサイトの直帰率やスクロール率などとの関係を研究してみたいですね。

常識から外れた雑談応答の検出

http://www.anlp.jp/proceedings/annual_meeting/2017/pdf_dir/P15-3.pdf

「Word2Vecモデルの学習においては、匿名掲示板である「おーぷん2ちゃんねる」から取得した人気記事7959記事を用いている。」と2ちゃんねるの表現から非常識語を見つけるという試みです。自社でUGCを運営している場合は役に立ちそうな研究です。

店舗レビューには何が書かれているか？–調査及びその自動分類–

http://www.anlp.jp/proceedings/annual_meeting/2017/pdf_dir/W-3.pdf

「店舗レビュー1,510件（5,266文）に対して人手でアノテーションを行い、このデータに基づいてトピックおよびその評価極性を分類するモデルを構築」と、恐ろしい数の文書を人出で分類したそうです。これもユーザーアンケートなどの研究に使えそうです。

トップダウン型共起グラフを用いたブログからの観光地の行動分析

http://www.anlp.jp/proceedings/annual_meeting/2017/pdf_dir/P20-3.pdf

「寄り道をする観光地を指定すると、その特徴を可視化する手法の開発を目的」とされています。マーケターであれば、関心のある商材を指定すると、その商材に関する特徴を可視化できるのではないでしょうか。テキストマイニングで言うところの共起グラフによる分析となります。

提案されている手法のステップは以下の通りで、自分の関心のある市場のデータを使えば同様のことができると思われます。
(1) 商材を指定しブログ記事集を取得する。
(2) 単語対の共起スコアを同時に出現した記事数をベースとして算出する。具体的には最低出現数を4回として単語を限定し、Jaccard 係数で共起スコアを算出する。
(3) 単語を頂点、共起スコアの得られた単語対を辺として共起グラフを生成する。
(4) 分析目的に応じて注目する頂点を指定する。
(5) 指定した頂点から1または2の距離で到達可能、かつ共起スコアが一定範囲内という辺および頂点を一定数を限度に残す。

社内データに基づくイノベータ人財のピックアップ

http://www.anlp.jp/proceedings/annual_meeting/2017/pdf_dir/P12-1.pdf

番外ではありますが、面白そうな取り組みなので紹介します。
社員プロフィールや業務報告（日報）や適性検査や360度調査やアンケートなどのデータをもとに、9種類の指標を定め、一部の指標のスコアリングを単語ベースで行うに際して、スコアリング用の単語を定めるためにWord2Vecを用いています。これからピープルアナリティクスを目指す企業には参考になる情報だと思います。

PythonやRを用いたアルゴリズム取引・株価分析まとめ

はじめに

まわりでシステムトレードや株価の機械学習による予測などに関心が高まってきたので、私も少し調べてみようと思いPythonやRで行われた分析・実装の事例を集めてみました。自分の資産を突っ込む気にはなれないですが、事例を知っておくだけ知っておきたいですね。

調査法

Google検索において以下のクエリで上位に表示されたサイトを中心にまとめました。
「python 機械学習株価」
「python 機械学習為替」
「python アルゴリズム取引」
「python machine learning stock price」
（同様にRも調べました。）

アルゴリズムトレードの理論から学ぶ

システムトレードで億万長者になるぞ! coursera で Computational Investing Part I を受けた
こちらはCourseraで開かれていたシステムトレード向けの講座を受けた方のブログです。Active Portfolio Management: A Quantitative Approach for Producing Superior Returns and Selecting Superior Returns and Controlling Riskはもっとも重要なアルゴリズムトレードの本とされているようです。出版年が1999年と結構古いですが。
Active Portfolio Management: A Quantitative Approach for Producing Superior Returns and Selecting Superior Returns and Controlling Risk (McGraw-Hill Library of Investment and Finance)

Rを用いた株価データのモニタリング

Stock Analysis using R
こちらはアルゴリズム取引とかではないですが、Rで株価を簡単に取得できるquantmodライブラリについて紹介されています。

library(quantmod)
getSymbols("MSFT")
chartSeries(MSFT, subset="last 3 months")
addBBands()

library(quantmod)

getSymbols("MSFT")

chartSeries(MSFT, subset="last 3 months")

addBBands()

以上のコードを実行しただけで時系列データがボックスプロットで描写されました。非常に便利そうです。複数のグラフを並べてディスプレーで見てみたいですね。図ではテクニカル分析で使われるボリンジャーバンドを一発で出してくれています。傾向を掴むための分析に役立つのではないでしょうか。

RSIで株価の連動性を見る

過去のデータからビッグデータ分析で株価を予測する
ここでは、主にpandasを用いて株価の分析を行っています。RSI（Relative Strength Index）という「一定の期間変動幅の中でどれ位株価が上昇しているのか、下落しているのかをはかるもの」を計算して日経平均との相関を見ています。データさえあれば、個別株の市場連動性を見るぶんにはpandasで十分に分析できそうです。

アルゴリズムトレードのシステム構成

自動トレードボット
こちらでは、仮想通貨取引のための自動トレードボットの作成のための手順などが書かれています。
工程としては
・仮想通貨価格データ取得
・バックテストの実施
・明日の価格予想
・学習パラメータの最適化
・結果のメール送信
・APIを使った成行トレード
・ジョブのスケジューリング
などとなっているようです。
完全自動化しようと思うと、作るのは大変そうですね。

WEB屋の自分が機械学習株価予想プログラムを開発した結果
こちらはよりシンプルで
「チャートを見る限り、誰がどうみても今日値上がりする銘柄」を検索して、毎朝Slackに通知してくれるシステム
をPythonで作成されています。最終的に資産を溶かす形になっているようです。

アルゴリズムトレード向けのPythonライブラリ

pythonのアルゴリズムトレードライブラリ
こちらでは
zipline（ http://www.zipline.io/ ）
PyAlgoTrade( http://gbeced.github.io/pyalgotrade/ )
pybacktest( https://github.com/ematvey/pybacktest )
backtrader( https://github.com/mementum/backtrader )
というライブラリが紹介されています。

遺伝的アルゴリズムで為替の自動売買

pythonと遺伝的アルゴリズムで作るFX自動売買システムその1
遺伝的アルゴリズムでFX自動売買その2 進化する売買AIの実装
 遺伝的アルゴリズムでFX自動売買その3 OandaAPIで実際に取引
こちらは、外国為替の取引をPythonで自動化させた試みです。最終的に負けてしまってはいますがシステム構成についても詳しく書かれているので、勉強になりました。データはOANDAというサービスが提供しているAPIを用いて取得し、GA（遺伝的アルゴリズム）を使って為替の売買タイミングを決めているようです。実際の売買にもOandaAPIというものを利用して完全に自動化させています。

決定木による株価リターンの予測

機械学習で未来を予測する – scikit-learn の決定木で未来の株価を予測
こちらはPythonでscikit-learnを用いて、決定木による株価の予測をされています。目的変数としてはリターンインデックスをおいています。前処理においてpandasが使われています。

ランダムフォレストによる株価の予測

Pythonで機械学習を使った株価予測のコードを書こう
こちらでは、ランダムフォレストを用いた機械学習で、ETFなどのデータを特徴量にして個別の株価予測を行っています。

アカデミックサイドでの研究事例

Stock Price Forecasting Using Information from Yahoo Finance and Google Trend
UC　Berkeleyの経済学部生の研究です。こちらはR言語で、ヤフーファイナンスとGoogleトレンドの情報を用いて株価を予測している研究です。これまでの時系列手法よりもパフォーマンスが良いそうです。

最後に

事例を集めてみて、どこまでのレベルのものに手を出すべきか悩ましいと思いました。完全自動化は維持するコストがかかりそうですし、ロスが増大しないか心配です。自分としては情報を自動で取得し、リターンが発生する確率の高そうな銘柄をサジェストしてくれるレベルで十分な気がします。

追記

人工知能学会の全国大会2017でAIを用いた株式運用に関する研究がなされていました。深層学習を用いたものが多そうです。
株価変動パターンの類似性を用いた株価予測
 深層学習と高頻度データを用いた株式注文状況の推定
 LSTM を用いた株価変動予測
 深層学習を用いた株価予測の分析

Rで株価の時系列データを簡単に集計する

はじめに

かねてからYahoo!ファイナンスなどの情報に、株価について期間を指定したリターンや標準偏差や回帰分析のあてはまりなどが無いので、そのようなデータを簡単に見れるようにしたいという願望がありました。加えて、できるだけ多くの銘柄の情報をみたいという願望もありました。Rを使えばその全てが簡単に実現できるのでチャレンジしてみました。

・データ収集
・データテーブルの整形
・成長率の計算
・標準偏差の計算
・自己回帰の計算
と非常にシンプルな工程ですが、3600銘柄分のこれらの情報を一気に手に入れれるというのは魅力的だなぁと思います。

データの収集

RFinanceYJというパッケージをインストールして、以下のrfinance_module.Rのスクリプトを置いておきます。

library(RFinanceYJ)

#API
quoteStockTsData &lt;- function(x, since=NULL,start.num=0,date.end=NULL,time.interval='daily')
{
  time.interval &lt;- substr(time.interval,1,1)
  function.stock &lt;- function(quote.table.item){
    if( xmlSize(quote.table.item) &lt; 5) return(NULL) 
    d &lt;- convertToDate(xmlValue(quote.table.item[[1]]),time.interval)
    o &lt;- as.number(xmlValue(quote.table.item[[2]]))
    h &lt;- as.number(xmlValue(quote.table.item[[3]]))
    l &lt;- as.number(xmlValue(quote.table.item[[4]]))
    c &lt;- as.number(xmlValue(quote.table.item[[5]]))
    v &lt;- ifelse(xmlSize(quote.table.item) &gt;= 6,as.number(xmlValue(quote.table.item[[6]])),0)
    a &lt;- ifelse(xmlSize(quote.table.item) &gt;= 7,as.number(xmlValue(quote.table.item[[7]])),0)
    return(data.frame(date=d,open=o,high=h,low=l,close=c,volume=v, adj_close=a))
  }
  return(quoteTsData(x,function.stock,since,start.num,date.end,time.interval,type="stock"))
}
quoteFundTsData &lt;- function(x, since=NULL,start.num=0,date.end=NULL,time.interval='daily')
{
  time.interval &lt;- substr(time.interval,1,1)
  function.fund &lt;- function(quote.table.item){
    d &lt;- convertToDate(xmlValue(quote.table.item[[1]]),time.interval)
    if(time.interval=='monthly'){
      d &lt;- endOfMonth(d)
    }
    c &lt;- as.number(xmlValue(quote.table.item[[2]]))
    v &lt;- as.number(xmlValue(quote.table.item[[3]]))
    return(data.frame(date=d,constant.value=c,NAV=v))
  }
  return(quoteTsData(x,function.fund,since,start.num,date.end,time.interval,type="fund"))
}
quoteFXTsData &lt;- function(x, since=NULL,start.num=0,date.end=NULL,time.interval='daily')
{
  time.interval &lt;- substr(time.interval,1,1)
  function.fx &lt;- function(quote.table.item){
    d &lt;- convertToDate(xmlValue(quote.table.item[[1]]),time.interval)
    o &lt;- as.number(xmlValue(quote.table.item[[2]]))
    h &lt;- as.number(xmlValue(quote.table.item[[3]]))
    l &lt;- as.number(xmlValue(quote.table.item[[4]]))
    c &lt;- as.number(xmlValue(quote.table.item[[5]]))
    return(data.frame(date=d,open=o,high=h,low=l,close=c))
  }
  return(quoteTsData(x,function.fx,since,start.num,date.end,time.interval,type="fx"))
}
######  private functions  #####
#get time series data from Yahoo! Finance.
quoteTsData &lt;- function(x,function.financialproduct,since,start.num,date.end,time.interval,type="stock"){
  r &lt;- NULL
  result.num &lt;- 51
  financial.data &lt;- data.frame(NULL)
  start &lt;- (gsub("([0-9]{4,4})-([0-9]{2,2})-([0-9]{2,2})","&amp;sy=\\1&amp;sm=\\2&amp;sd=\\3",since))
  end   &lt;- (gsub("([0-9]{4,4})-([0-9]{2,2})-([0-9]{2,2})","&amp;ey=\\1&amp;em=\\2&amp;ed=\\3",date.end))
  
  if(!any(time.interval==c('d','w','m'))) stop("Invalid time.interval value")
  
  extractQuoteTable &lt;- function(r,type){
    if(type %in% c("fund","fx")){
      tbl &lt;- r[[2]][[2]][[7]][[3]][[3]][[9]][[2]]
    }
    else{
      tbl &lt;- r[[2]][[2]][[7]][[3]][[3]][[10]][[2]]
    }
    return(tbl)
  }
  
  #while( result.num &gt;= 51 ){
  while(1){
    start.num &lt;- start.num + 1
    quote.table &lt;- NULL
    quote.url &lt;- paste('http://info.finance.yahoo.co.jp/history/?code=',x,start,end,'&amp;p=',start.num,'&amp;tm=',substr(time.interval,1,1),sep="")
    try(r&lt;-htmlParse(quote.url))
    if( is.null(r) ) stop(paste("Can not access :", quote.url))
    
    try( quote.table &lt;- xpathApply(r,"//table")[[2]], TRUE )
    
    quote.size&lt;-xmlSize(quote.table)
    #cat(paste("size:",quote.size))
    if(xmlSize(quote.table)&lt;=1){
      return (financial.data)
    }
    if( is.null(quote.table) ){
      if( is.null(financial.data) ){
        stop(paste("Can not quote :", x))
      }else{
        financial.data &lt;- financial.data[order(financial.data$date),]
        return(financial.data)
      }
    }
    
    size &lt;- xmlSize(quote.table)
    for(i in 2:size){
      financial.data &lt;- rbind(financial.data,function.financialproduct(quote.table[[i]]))
    }
    
    #result.num &lt;- xmlSize(quote.table)
    Sys.sleep(1)
  }
  financial.data &lt;- financial.data[order(financial.data$date),]
  return(financial.data)  
}
#convert string formart date to POSIXct object
convertToDate &lt;- function(date.string,time.interval)
{
  #data format is different between monthly and dialy or weekly
  if(any(time.interval==c('d','w'))){
    result &lt;- gsub("^([0-9]{4})([^0-9]+)([0-9]{1,2})([^0-9]+)([0-9]{1,2})([^0-9]+)","\\1-\\3-\\5",date.string)
  }else if(time.interval=='m'){
    result &lt;- gsub("^([0-9]{4})([^0-9]+)([0-9]{1,2})([^0-9]+)","\\1-\\3-01",date.string)
  }
  return(as.POSIXct(result))
}
#convert string to number.
as.number &lt;- function(string)
{
  return(as.double(as.character(gsub("[^0-9.]", "",string))))
}
#return end of month date.
endOfMonth &lt;- function(date.obj)
{
  startOfMonth     &lt;- as.Date(format(date.obj,"%Y%m01"),"%Y%m%d")
  startOfNextMonth &lt;- as.Date(format(startOfMonth+31,"%Y%m01"),"%Y%m%d")
  return(startOfNextMonth-1)
}</pre>

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

library(RFinanceYJ)

#API

quoteStockTsData <- function(x, since=NULL,start.num=0,date.end=NULL,time.interval='daily')

{

time.interval <- substr(time.interval,1,1)

function.stock <- function(quote.table.item){

if( xmlSize(quote.table.item) < 5) return(NULL)

d <- convertToDate(xmlValue(quote.table.item[[1]]),time.interval)

o <- as.number(xmlValue(quote.table.item[[2]]))

h <- as.number(xmlValue(quote.table.item[[3]]))

l <- as.number(xmlValue(quote.table.item[[4]]))

c <- as.number(xmlValue(quote.table.item[[5]]))

v <- ifelse(xmlSize(quote.table.item) >= 6,as.number(xmlValue(quote.table.item[[6]])),0)

a <- ifelse(xmlSize(quote.table.item) >= 7,as.number(xmlValue(quote.table.item[[7]])),0)

return(data.frame(date=d,open=o,high=h,low=l,close=c,volume=v, adj_close=a))

}

return(quoteTsData(x,function.stock,since,start.num,date.end,time.interval,type="stock"))

}

quoteFundTsData <- function(x, since=NULL,start.num=0,date.end=NULL,time.interval='daily')

{

time.interval <- substr(time.interval,1,1)

function.fund <- function(quote.table.item){

d <- convertToDate(xmlValue(quote.table.item[[1]]),time.interval)

if(time.interval=='monthly'){

d <- endOfMonth(d)

}

c <- as.number(xmlValue(quote.table.item[[2]]))

v <- as.number(xmlValue(quote.table.item[[3]]))

return(data.frame(date=d,constant.value=c,NAV=v))

}

return(quoteTsData(x,function.fund,since,start.num,date.end,time.interval,type="fund"))

}

quoteFXTsData <- function(x, since=NULL,start.num=0,date.end=NULL,time.interval='daily')

{

time.interval <- substr(time.interval,1,1)

function.fx <- function(quote.table.item){

d <- convertToDate(xmlValue(quote.table.item[[1]]),time.interval)

o <- as.number(xmlValue(quote.table.item[[2]]))

h <- as.number(xmlValue(quote.table.item[[3]]))

l <- as.number(xmlValue(quote.table.item[[4]]))

c <- as.number(xmlValue(quote.table.item[[5]]))

return(data.frame(date=d,open=o,high=h,low=l,close=c))

}

return(quoteTsData(x,function.fx,since,start.num,date.end,time.interval,type="fx"))

}

###### private functions #####

#get time series data from Yahoo! Finance.

quoteTsData <- function(x,function.financialproduct,since,start.num,date.end,time.interval,type="stock"){

r <- NULL

result.num <- 51

financial.data <- data.frame(NULL)

start <- (gsub("([0-9]{4,4})-([0-9]{2,2})-([0-9]{2,2})","&sy=\\1&sm=\\2&sd=\\3",since))

end <- (gsub("([0-9]{4,4})-([0-9]{2,2})-([0-9]{2,2})","&ey=\\1&em=\\2&ed=\\3",date.end))

if(!any(time.interval==c('d','w','m'))) stop("Invalid time.interval value")

extractQuoteTable <- function(r,type){

if(type %in% c("fund","fx")){

tbl <- r[[2]][[2]][[7]][[3]][[3]][[9]][[2]]

}

else{

tbl <- r[[2]][[2]][[7]][[3]][[3]][[10]][[2]]

}

return(tbl)

}

#while( result.num >= 51 ){

while(1){

start.num <- start.num + 1

quote.table <- NULL

quote.url <- paste('http://info.finance.yahoo.co.jp/history/?code=',x,start,end,'&p=',start.num,'&tm=',substr(time.interval,1,1),sep="")

try(r<-htmlParse(quote.url))

if( is.null(r) ) stop(paste("Can not access :", quote.url))

try( quote.table <- xpathApply(r,"//table")[[2]], TRUE )

quote.size<-xmlSize(quote.table)

#cat(paste("size:",quote.size))

if(xmlSize(quote.table)<=1){

return (financial.data)

}

if( is.null(quote.table) ){

if( is.null(financial.data) ){

stop(paste("Can not quote :", x))

}else{

financial.data <- financial.data[order(financial.data$date),]

return(financial.data)

}

size <- xmlSize(quote.table)

for(i in 2:size){

financial.data <- rbind(financial.data,function.financialproduct(quote.table[[i]]))

}

#result.num <- xmlSize(quote.table)

Sys.sleep(1)

}

financial.data <- financial.data[order(financial.data$date),]

return(financial.data)

}

#convert string formart date to POSIXct object

convertToDate <- function(date.string,time.interval)

{

#data format is different between monthly and dialy or weekly

if(any(time.interval==c('d','w'))){

result <- gsub("^([0-9]{4})([^0-9]+)([0-9]{1,2})([^0-9]+)([0-9]{1,2})([^0-9]+)","\\1-\\3-\\5",date.string)

}else if(time.interval=='m'){

result <- gsub("^([0-9]{4})([^0-9]+)([0-9]{1,2})([^0-9]+)","\\1-\\3-01",date.string)

}

return(as.POSIXct(result))

}

#convert string to number.

as.number <- function(string)

{

return(as.double(as.character(gsub("[^0-9.]", "",string))))

}

#return end of month date.

endOfMonth <- function(date.obj)

{

startOfMonth <- as.Date(format(date.obj,"%Y%m01"),"%Y%m%d")

startOfNextMonth <- as.Date(format(startOfMonth+31,"%Y%m01"),"%Y%m%d")

return(startOfNextMonth-1)

}</pre>

以下のコードで2016-01-01から現在に至るまでの、指定した株価コードの株価情報を集めてきます。3600銘柄の株価コードを用意する必要があるので、集めてcodelist.csvで保存しておきます。銘柄コードの情報はここでは掲載しません。（適当なサイトにころがっています。）

source("rfinance_module.R")

codelist <- read.csv("codelist.csv")
since    <- "2016-01-01"

stock_price_database <- 0
for(i in 1:nrow(codelist)){
  if(i ==1){
    dataset <- quoteStockTsData(paste(codelist[i,1],".t",sep = ""),since=since)
    code    <- matrix(codelist[i,1], nrow(dataset), 1)
    dataset <- cbind(code, dataset) 
    stock_price_database <- dataset
  }
  else{
    dataset <- quoteStockTsData(paste(codelist[i,1],".t",sep = ""),since=since)
    code    <- matrix(codelist[i,1], nrow(dataset), 1)
    dataset <- cbind(code, dataset) 
    stock_price_database <- rbind(stock_price_database,dataset)    
  }
}

source("rfinance_module.R")

codelist <- read.csv("codelist.csv")

since <- "2016-01-01"

stock_price_database <- 0

for(i in 1:nrow(codelist)){

if(i ==1){

dataset <- quoteStockTsData(paste(codelist[i,1],".t",sep = ""),since=since)

code <- matrix(codelist[i,1], nrow(dataset), 1)

dataset <- cbind(code, dataset)

stock_price_database <- dataset

}

else{

dataset <- quoteStockTsData(paste(codelist[i,1],".t",sep = ""),since=since)

code <- matrix(codelist[i,1], nrow(dataset), 1)

dataset <- cbind(code, dataset)

stock_price_database <- rbind(stock_price_database,dataset)

}

以下のコードを実行して、元のデータを調整済み終値について、縦が日付で横が銘柄のテーブルに整形します。

#make stock price data table
library(dplyr)
code_list <- data.frame(code=distinct(stock_price_database,code)$code)
date_list <- data.frame(date=distinct(stock_price_database,date)$date)

stock_dataset <- date_list
stock_price_database_adj <- stock_price_database[,which(colnames(stock_price_database) %in% c("code","date","adj_close"))]

for( i in 1:nrow(code_list)){ #nrow(code_list)
  each_stock <- filter(stock_price_database_adj, code == code_list$code[i])
  each_stock <- each_stock[,-1]
  each_stock <- distinct(each_stock, date,.keep_all = TRUE)
  colnames(each_stock)[2] <- paste("sdj",code_list$code[i],sep = "") 
  stock_dataset <- left_join(stock_dataset, each_stock, by = "date")
}

#make stock price data table

library(dplyr)

code_list <- data.frame(code=distinct(stock_price_database,code)$code)

date_list <- data.frame(date=distinct(stock_price_database,date)$date)

stock_dataset <- date_list

stock_price_database_adj <- stock_price_database[,which(colnames(stock_price_database) %in% c("code","date","adj_close"))]

for( i in 1:nrow(code_list)){ #nrow(code_list)

each_stock <- filter(stock_price_database_adj, code == code_list$code[i])

each_stock <- each_stock[,-1]

each_stock <- distinct(each_stock, date,.keep_all = TRUE)

colnames(each_stock)[2] <- paste("sdj",code_list$code[i],sep = "")

stock_dataset <- left_join(stock_dataset, each_stock, by = "date")

}

成長率

指定した期間での成長率を計算します。

#calculate growth rate

#start_date
start_date <- "2016-01-04"
#end_date
end_date <- "2017-02-17"

#return rate datatable
return_data <- filter(stock_dataset, date==start_date |date==end_date)
return_data <- return_data[,-1]
return_data <- rbind(return_data, (as.numeric(return_data[1,]) / as.numeric(return_data[2,])-1)*100)
return_data <- data.frame(code=colnames(return_data),t(return_data))
colnames(return_data) <- c("code",end_date,start_date,"return_rate")
return_data <- return_data[order(return_data$return_rate,decreasing = TRUE),]

#calculate growth rate

#start_date

start_date <- "2016-01-04"

#end_date

end_date <- "2017-02-17"

#return rate datatable

return_data <- filter(stock_dataset, date==start_date |date==end_date)

return_data <- return_data[,-1]

return_data <- rbind(return_data, (as.numeric(return_data[1,]) / as.numeric(return_data[2,])-1)*100)

return_data <- data.frame(code=colnames(return_data),t(return_data))

colnames(return_data) <- c("code",end_date,start_date,"return_rate")

return_data <- return_data[order(return_data$return_rate,decreasing = TRUE),]

標準偏差

ここではシンプルに標準偏差の計算を行います。

#volatility
volatility_data <- stock_dataset[,-1]
volatility_data <- log(volatility_data[-nrow(volatility_data),]) - log(volatility_data[-1,])

volatility <- data.frame(matrix(0,1,2))
colnames(volatility) <- c("code", "stdev")

for(i in 1:length(volatility_data)){
  volatility[i,1] <- colnames(volatility_data)[i]
  volatility[i,2] <- sd(volatility_data[,i],na.rm = TRUE)
}

#volatility

volatility_data <- stock_dataset[,-1]

volatility_data <- log(volatility_data[-nrow(volatility_data),]) - log(volatility_data[-1,])

volatility <- data.frame(matrix(0,1,2))

colnames(volatility) <- c("code", "stdev")

for(i in 1:length(volatility_data)){

volatility[i,1] <- colnames(volatility_data)[i]

volatility[i,2] <- sd(volatility_data[,i],na.rm = TRUE)

}

自己回帰による決定係数

自己回帰でフィッティングの良い銘柄はどれなのかを単純に知りたいので、3600銘柄分計算させます。

#AR(1) estimation

ar1_estimation <- stock_dataset[,-1]
ar1_rsuares <- data.frame(matrix(0,1,2))
colnames(ar1_rsuares) <- c("code", "adj_rsquare")

for(i in 1:length(ar1_estimation)){
  #objective variable
  ar1_y <- ar1_estimation[-nrow(ar1_estimation),i]
  #input variable
  ar1_x <- ar1_estimation[-1,i]
  #combine data
  ar1_data <- data.frame(cbind(ar1_y, ar1_x))
  
  #estimation
  if(class(try(lm(data = ar1_data, formula = ar1_y ~. )))=="lm"){
    model <- lm(data = ar1_data, formula = ar1_y ~. )
    model_sum <- summary(model)
    
    ar1_rsuares[i,1] <- colnames(ar1_estimation)[i]
    ar1_rsuares[i,2] <- model_sum$adj.r.squared    
  }
  else{
    ar1_rsuares[i,1] <- colnames(ar1_estimation)[i]
    ar1_rsuares[i,2] <- "NULL"
  }
}

#AR(1) estimation

ar1_estimation <- stock_dataset[,-1]

ar1_rsuares <- data.frame(matrix(0,1,2))

colnames(ar1_rsuares) <- c("code", "adj_rsquare")

for(i in 1:length(ar1_estimation)){

#objective variable

ar1_y <- ar1_estimation[-nrow(ar1_estimation),i]

#input variable

ar1_x <- ar1_estimation[-1,i]

#combine data

ar1_data <- data.frame(cbind(ar1_y, ar1_x))

#estimation

if(class(try(lm(data = ar1_data, formula = ar1_y ~. )))=="lm"){

model <- lm(data = ar1_data, formula = ar1_y ~. )

model_sum <- summary(model)

ar1_rsuares[i,1] <- colnames(ar1_estimation)[i]

ar1_rsuares[i,2] <- model_sum$adj.r.squared

}

else{

ar1_rsuares[i,1] <- colnames(ar1_estimation)[i]

ar1_rsuares[i,2] <- "NULL"

}

以上より、見てみたかった指標の3600銘柄分を一気に計算することができました。

可視化

試しに一番成長率の高かった株価をプロットしてみます。

#visualization
library(ggplot2)
stock_code <- "sdj7612"

df <- data.frame(stock_dataset[,which(colnames(stock_dataset) %in% c("date", stock_code))])

g <- ggplot(df,aes(x= date, y= df[,2]))+
     geom_line(colour = "magenta",linetype = 2,size = 0.5)+
     geom_smooth (method = "lm")+
     xlab("Date")+
     ylab("Stock Price")+
     ggtitle("Stock Price")
plot(g)

#visualization

library(ggplot2)

stock_code <- "sdj7612"

df <- data.frame(stock_dataset[,which(colnames(stock_dataset) %in% c("date", stock_code))])

g <- ggplot(df,aes(x= date, y= df[,2]))+

geom_line(colour = "magenta",linetype = 2,size = 0.5)+

geom_smooth (method = "lm")+

xlab("Date")+

ylab("Stock Price")+

ggtitle("Stock Price")

plot(g)

おわりに

株価データにカジュアルにアクセスできる環境を手に入れることができたので、今後は様々な指標や時系列解析・業界ごとの分析を試してみたいと思います。

顧客生涯価値(CLV)の計算 with R

顧客生涯価値(CLV：Customer Lifetime Value)を計算してくれるRのコード（Calculating Customer Lifetime Value with Recency, Frequency, and Monetary (RFM)）があったので、今更感がありますが取り上げたいと思います。

目次

・顧客生涯価値の数式
・データセット
・関数
・データセットの読み込みと加工
・再購買率とRFMとの関係
・再購買率の推定
・顧客生涯価値の計算
・参考情報

顧客生涯価値の数式

まず、顧客生涯価値の数式は以下の通りです。

t：年や月などの期間
n：顧客が解約するまでの期間合計
r：保持率（1-解約率）
P(t)：t期に顧客から得られる収益
d：割引率

rは数式上では固定ですが、実際にはデモグラ属性（年齢、地理情報、職種など）や行動（RFMなど）や在職中かどうかなどの要因により変わりうるものだと考えられます。参考文献のブログでは、このrのロジスティック回帰による推定がなされています。

データセット

データ名：CDNow
概要：1997年の第一四半期をスタート時点とした顧客の購買行動データ
期間：1997年1月〜1998年6月
顧客数：23570
取引レコード数：69659
変数：顧客ID、購入日、購入金額
入手方法：DatasetsでCDNOW dataset (full dataset)をダウンロード

関数

参考文献にはgetDataFrame関数、getPercentages関数、getCLV関数の三つの関数が出てきますが、CLVの計算に必要なのはgetDataFrame関数、getCLV関数の二つです。getPercentages関数はRecencyなどに応じて細かく分析する際に用います。

getDataFrame関数・・・生のデータセットから、指定した期間に応じたRecencyのデータを作成する関数です。

getDataFrame <- function(df,startDate,endDate,tIDColName="ID",tDateColName="Date",tAmountColName="Amount"){
  
  #日付について降順でデータフレームを並び替える
  df <- df[order(df[,tDateColName],decreasing = TRUE),]
  
  #開始期間の前と終了期間の後のデータを取り除く
  df <- df[df[,tDateColName]>= startDate,]
  df <- df[df[,tDateColName]<= endDate,]
  
  #重複したIDの行を除外して新しいデータフレームを生成する
  newdf <- df[!duplicated(df[,tIDColName]),]
  
  # 終了期間に対してのRecencyを計算する（最も小さい値は最も最近であることを示す）
  Recency <- as.numeric(difftime(endDate,newdf[,tDateColName],units="days"))
  
  # newdfデータフレームに日付の列を加える
  newdf <- cbind(newdf,Recency)
  
  # table関数やtapply関数の返す順番と揃うようにデータフレームをIDごとに並べ替える
  newdf <- newdf[order(newdf[,tIDColName]),]
  
  # 頻度を計算する
  fre <- as.data.frame(table(df[,tIDColName]))
  Frequency <- fre[,2]
  newdf <- cbind(newdf,Frequency)
  
  # 取引ごとの金額を計算する
  m <- as.data.frame(tapply(df[,tAmountColName],df[,tIDColName],sum))
  Monetary <- m[,1]/Frequency
  newdf <- cbind(newdf,Monetary)
  
  return(newdf)
  
}

getDataFrame <- function(df,startDate,endDate,tIDColName="ID",tDateColName="Date",tAmountColName="Amount"){

#日付について降順でデータフレームを並び替える

df <- df[order(df[,tDateColName],decreasing = TRUE),]

#開始期間の前と終了期間の後のデータを取り除く

df <- df[df[,tDateColName]>= startDate,]

df <- df[df[,tDateColName]<= endDate,]

#重複したIDの行を除外して新しいデータフレームを生成する

newdf <- df[!duplicated(df[,tIDColName]),]

# 終了期間に対してのRecencyを計算する（最も小さい値は最も最近であることを示す）

Recency <- as.numeric(difftime(endDate,newdf[,tDateColName],units="days"))

# newdfデータフレームに日付の列を加える

newdf <- cbind(newdf,Recency)

# table関数やtapply関数の返す順番と揃うようにデータフレームをIDごとに並べ替える

newdf <- newdf[order(newdf[,tIDColName]),]

# 頻度を計算する

fre <- as.data.frame(table(df[,tIDColName]))

Frequency <- fre[,2]

newdf <- cbind(newdf,Frequency)

# 取引ごとの金額を計算する

m <- as.data.frame(tapply(df[,tAmountColName],df[,tIDColName],sum))

Monetary <- m[,1]/Frequency

newdf <- cbind(newdf,Monetary)

return(newdf)

}

getPercentages関数・・・Recencyなどの回数に応じて、購入した顧客の割合などを計算するための関数です。

require(plyr)
getPercentages <- function(df,colNames){
  #指定した列の名前とBuyからなるベクトルの作成
  Var <- c(colNames,"Buy")
  #指定した変数とBuyが合致する列だけを抽出
  df <- df[,names(df) %in% Var,drop=F]
  
  #購買した回数に応じたデータ数を返す
  a <- ddply(df,Var,summarize,Number=length(Buy))
  #回数に応じたデータ数に対して、全体の割合を返す
  b <- ddply(a,
             .(),
             .fun=function(x){
               transform(x, Percentage=with(x,round(ave(Number,a[,names(a) %in% Var,drop=F],FUN=sum)/ave(Number,a[,names(a) %in% colNames,drop=F],FUN=sum),2)))
             })
  #1列目を除外する
  b <- b[b$Buy==1,-1]
  
  return(b)
  
}

require(plyr)

getPercentages <- function(df,colNames){

#指定した列の名前とBuyからなるベクトルの作成

Var <- c(colNames,"Buy")

#指定した変数とBuyが合致する列だけを抽出

df <- df[,names(df) %in% Var,drop=F]

#購買した回数に応じたデータ数を返す

a <- ddply(df,Var,summarize,Number=length(Buy))

#回数に応じたデータ数に対して、全体の割合を返す

b <- ddply(a,

.(),

.fun=function(x){

transform(x, Percentage=with(x,round(ave(Number,a[,names(a) %in% Var,drop=F],FUN=sum)/ave(Number,a[,names(a) %in% colNames,drop=F],FUN=sum),2)))

})

#1列目を除外する

b <- b[b$Buy==1,-1]

return(b)

}

getCLV関数・・・Recency、Frequency、Monetary、購入者の数（1人と置く）、コスト（0としている）、期間、割引率、推定したモデルをもとにCLVを計算する関数です。

getCLV<-function(r,f,m,n,cost,periods,dr,pModel){
  
  df<-data.frame(period=c(0),r=c(r),f=c(f),n=c(n),value=c(0))
  
  for(i in 1:periods){
    backstep<-df[df$period==i-1,]
    nrow<-nrow(backstep)
    for(j in 1:nrow){
      r<-backstep[j,]$r
      f<-backstep[j,]$f
      n<-backstep[j,]$n
      p<-predict(pModel,data.frame(Recency=r,Frequency=f),type='response')[1]
      buyers<-n*p
      df<-rbind(df,c(i,0,f+1,buyers,buyers*(m-cost) / (1+dr)^i))
      df<-rbind(df,c(i,r+1,f,n-buyers,(n-buyers)*(-cost)  / (1+dr)^i ))
    }
  }
  
  return(sum(df$value))
  
}

getCLV<-function(r,f,m,n,cost,periods,dr,pModel){

df<-data.frame(period=c(0),r=c(r),f=c(f),n=c(n),value=c(0))

for(i in 1:periods){

backstep<-df[df$period==i-1,]

nrow<-nrow(backstep)

for(j in 1:nrow){

r<-backstep[j,]$r

f<-backstep[j,]$f

n<-backstep[j,]$n

p<-predict(pModel,data.frame(Recency=r,Frequency=f),type='response')[1]

buyers<-n*p

df<-rbind(df,c(i,0,f+1,buyers,buyers*(m-cost) / (1+dr)^i))

df<-rbind(df,c(i,r+1,f,n-buyers,(n-buyers)*(-cost) / (1+dr)^i ))

}

return(sum(df$value))

}

データセットの読み込みと加工

#　CDNOW_SAMPLE.txtを読み込む
df <- read.table(file = "CDNOW_master.txt",header=F)

# 取引ごと且つ顧客ごとに、ID、取引日、購入金額の列からなるデータフレームを生成する
df <- as.data.frame(cbind(df[,1],df[,2],df[,4]))

# 列名の変更
names <- c("ID","Date","Amount")
names(df) <- names

#日付の型の変更
df[,2] <- as.Date(as.character(df[,2]),"%Y%m%d")

#　CDNOW_SAMPLE.txtを読み込む

df <- read.table(file = "CDNOW_master.txt",header=F)

# 取引ごと且つ顧客ごとに、ID、取引日、購入金額の列からなるデータフレームを生成する

df <- as.data.frame(cbind(df[,1],df[,2],df[,4]))

# 列名の変更

names <- c("ID","Date","Amount")

names(df) <- names

#日付の型の変更

df[,2] <- as.Date(as.character(df[,2]),"%Y%m%d")

> head(df)
  ID       Date Amount
1  1 1997-01-01  11.77
2  2 1997-01-12  12.00
3  2 1997-01-12  77.00
4  3 1997-01-02  20.76
5  3 1997-03-30  20.76
6  3 1997-04-02  19.54

> head(df)

ID Date Amount

1 1 1997-01-01 11.77

2 2 1997-01-12 12.00

3 2 1997-01-12 77.00

4 3 1997-01-02 20.76

5 3 1997-03-30 20.76

6 3 1997-04-02 19.54

再購買率とRFMとの関係

まず初めにデータセットを加工します。ロジットの推定における説明変数用のデータとして19970101〜19980228のデータを用い、被説明変数にあたる購入したかどうかのデータを19980301〜19980430のデータを用いて作ります。

# 取引期間の範囲を指定する
startDate_history <- as.Date("19970101","%Y%m%d")
endDate_history <- as.Date("19980228","%Y%m%d")

# 2ヶ月間の購入サイクルにあたる予測期間を指定する
startDate_forcast <- as.Date("19980301","%Y%m%d")
endDate_forcast <- as.Date("19980430","%Y%m%d")

# getDataFrame関数でRFMデータに変換する
history <- getDataFrame(df,startDate_history,endDate_history)
forcast <- getDataFrame(df,startDate_forcast,endDate_forcast)

# 購入サイクルを60日とし、Recencyを離散化する
history$Recency <- history$Recency %/% 60 

# Monetaryを10ドル単位で離散化する
breaks <- seq(0,round(max(history$Monetary)+9),by=10)
history$Monetary <- as.numeric(cut(history$Monetary,breaks,labels=FALSE))

# RFMのデータフレームに購入・非購入データを加える
Buy <- rep(0,nrow(history))
history <- cbind(history,Buy)

# Buyに関して、予測期間（1998/3/1~1998/4/30）の期間で再度購入しているならば1を返すようにする
history[history$ID %in% forcast$ID, ]$Buy<-1
train　<-　history

# 取引期間の範囲を指定する

startDate_history <- as.Date("19970101","%Y%m%d")

endDate_history <- as.Date("19980228","%Y%m%d")

# 2ヶ月間の購入サイクルにあたる予測期間を指定する

startDate_forcast <- as.Date("19980301","%Y%m%d")

endDate_forcast <- as.Date("19980430","%Y%m%d")

# getDataFrame関数でRFMデータに変換する

history <- getDataFrame(df,startDate_history,endDate_history)

forcast <- getDataFrame(df,startDate_forcast,endDate_forcast)

# 購入サイクルを60日とし、Recencyを離散化する

history$Recency <- history$Recency %/% 60

# Monetaryを10ドル単位で離散化する

breaks <- seq(0,round(max(history$Monetary)+9),by=10)

history$Monetary <- as.numeric(cut(history$Monetary,breaks,labels=FALSE))

# RFMのデータフレームに購入・非購入データを加える

Buy <- rep(0,nrow(history))

history <- cbind(history,Buy)

# Buyに関して、予測期間（1998/3/1~1998/4/30）の期間で再度購入しているならば1を返すようにする

history[history$ID %in% forcast$ID, ]$Buy<-1

train　<-　history

データを確認します。

> head(train)
   ID       Date Amount Recency Frequency Monetary Buy
1   1 1997-01-01  11.77       7         1        2   0
2   2 1997-01-12  12.00       6         2        5   0
8   3 1997-11-25  20.96       1         5        3   0
13  4 1997-12-12  26.48       1         4        3   0
24  5 1998-01-03  37.47       0        11        4   0
25  6 1997-01-01  20.99       7         1        3   0

> head(train)

ID Date Amount Recency Frequency Monetary Buy

1 1 1997-01-01 11.77 7 1 2 0

2 2 1997-01-12 12.00 6 2 5 0

8 3 1997-11-25 20.96 1 5 3 0

13 4 1997-12-12 26.48 1 4 3 0

24 5 1998-01-03 37.47 0 11 4 0

25 6 1997-01-01 20.99 7 1 3 0

# getPercentages関数を用いて、Recency変数に基づいて購入者の割合を計算する
colNames <- c("Recency")
p <- getPercentages(train,colNames)

# getPercentages関数を用いて、Recency変数に基づいて購入者の割合を計算する

colNames <- c("Recency")

p <- getPercentages(train,colNames)

> p
   Recency Buy Number Percentage
2        0   1   1180       0.45
4        1   1    581       0.28
6        2   1    279       0.22
8        3   1    198       0.17
10       4   1    163       0.14
12       5   1    249       0.05
14       6   1    316       0.03
16       7   1     13       0.03

> p

Recency Buy Number Percentage

2 0 1 1180 0.45

4 1 1 581 0.28

6 2 1 279 0.22

8 3 1 198 0.17

10 4 1 163 0.14

12 5 1 249 0.05

14 6 1 316 0.03

16 7 1 13 0.03

60日以内に購入した顧客（Recency=0）のうち、45%が再び購入しているようです。

# getPercentages関数を用いて、Monetary変数に基づいて購入者の割合を計算する
colNames <- c("Monetary")
p <- getPercentages(train,colNames)

# getPercentages関数を用いて、Monetary変数に基づいて購入者の割合を計算する

colNames <- c("Monetary")

p <- getPercentages(train,colNames)

> head(p,10)
   Monetary Buy Number Percentage
2         1   1     43       0.05
4         2   1    727       0.09
6         3   1    760       0.15
8         4   1    525       0.17
10        5   1    333       0.16
12        6   1    218       0.17
14        7   1    117       0.16
16        8   1     87       0.16
18        9   1     45       0.13
20       10   1     44       0.19

> head(p,10)

Monetary Buy Number Percentage

2 1 1 43 0.05

4 2 1 727 0.09

6 3 1 760 0.15

8 4 1 525 0.17

10 5 1 333 0.16

12 6 1 218 0.17

14 7 1 117 0.16

16 8 1 87 0.16

18 9 1 45 0.13

20 10 1 44 0.19

100ドル購入した顧客（Monetary=10）のうち、19%が再び購入しているようです。

# getPercentages関数を用いて、Frequency変数に基づいて購入者の割合を計算する
colNames<-c("Frequency")
p<-getPercentages(train,colNames)

# getPercentages関数を用いて、Frequency変数に基づいて購入者の割合を計算する

colNames<-c("Frequency")

p<-getPercentages(train,colNames)

> head(p,10)
   Frequency Buy Number Percentage
2          1   1    399       0.03
4          2   1    455       0.11
6          3   1    411       0.18
8          4   1    299       0.23
10         5   1    293       0.32
12         6   1    212       0.38
14         7   1    164       0.42
16         8   1    123       0.41
18         9   1    108       0.50
20        10   1     70       0.49

> head(p,10)

Frequency Buy Number Percentage

2 1 1 399 0.03

4 2 1 455 0.11

6 3 1 411 0.18

8 4 1 299 0.23

10 5 1 293 0.32

12 6 1 212 0.38

14 7 1 164 0.42

16 8 1 123 0.41

18 9 1 108 0.50

20 10 1 70 0.49

10回購入したことのある顧客（Frequency=10）は49%が再び購入しているようです。

再購買率の推定

RFM（Recency、Frequency、Monetary）のデータに基づいて、再購買率をロジスティック回帰によって推定し、予測確率を用いて顧客生涯価値を計算します。

# 再購買の割合をRecencyでロジスティック回帰
r.glm = glm(Percentage~Recency,family=quasibinomial(link='logit'),data=p)
p_r <- p

# 再購買の割合をRecencyでロジスティック回帰

r.glm = glm(Percentage~Recency,family=quasibinomial(link='logit'),data=p)

p_r <- p

> summary(r.glm)

Call:
glm(formula = Percentage ~ Recency, family = quasibinomial(link = "logit"), 
    data = p)

Deviance Residuals: 
      Min         1Q     Median         3Q        Max  
-0.085865  -0.067551   0.000333   0.041373   0.116927  

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.28538    0.11395  -2.504   0.0462 *  
Recency     -0.47037    0.04035 -11.657  2.4e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasibinomial family taken to be 0.005805806)

    Null deviance: 1.052105  on 7  degrees of freedom
Residual deviance: 0.034439  on 6  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 6

> summary(r.glm)

Call:

glm(formula = Percentage ~ Recency, family = quasibinomial(link = "logit"),

data = p)

Deviance Residuals:

Min 1Q Median 3Q Max

-0.085865 -0.067551 0.000333 0.041373 0.116927

Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) -0.28538 0.11395 -2.504 0.0462 *

Recency -0.47037 0.04035 -11.657 2.4e-05 ***

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasibinomial family taken to be 0.005805806)

Null deviance: 1.052105 on 7 degrees of freedom

Residual deviance: 0.034439 on 6 degrees of freedom

AIC: NA

Number of Fisher Scoring iterations: 6

# 再購買の割合を頻度でロジスティック回帰
f.glm = glm(Percentage~Frequency,family=quasibinomial(link='logit'),data=p)
p_f <- p

# 再購買の割合を頻度でロジスティック回帰

f.glm = glm(Percentage~Frequency,family=quasibinomial(link='logit'),data=p)

p_f <- p

> summary(f.glm)

Call:
glm(formula = Percentage ~ Frequency, family = quasibinomial(link = "logit"), 
    data = p)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-0.79432  -0.03744   0.10447   0.19982   0.52309  

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -1.48217    0.21490  -6.897 6.05e-09 ***
Frequency    0.12601    0.01139  11.067 1.67e-15 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasibinomial family taken to be 0.08078336)

    Null deviance: 25.652  on 55  degrees of freedom
Residual deviance:  4.053  on 54  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 7

> summary(f.glm)

Call:

glm(formula = Percentage ~ Frequency, family = quasibinomial(link = "logit"),

data = p)

Deviance Residuals:

Min 1Q Median 3Q Max

-0.79432 -0.03744 0.10447 0.19982 0.52309

Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) -1.48217 0.21490 -6.897 6.05e-09 ***

Frequency 0.12601 0.01139 11.067 1.67e-15 ***

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasibinomial family taken to be 0.08078336)

Null deviance: 25.652 on 55 degrees of freedom

Residual deviance: 4.053 on 54 degrees of freedom

AIC: NA

Number of Fisher Scoring iterations: 7

# 再購買の割合をMonetaryでロジスティック回帰
m.glm = glm(Percentage~Monetary,family=quasibinomial(link='logit'),data=p)
p_m <- p

# 再購買の割合をMonetaryでロジスティック回帰

m.glm = glm(Percentage~Monetary,family=quasibinomial(link='logit'),data=p)

p_m <- p

> summary(m.glm)

Call:
glm(formula = Percentage ~ Monetary, family = quasibinomial(link = "logit"), 
    data = p)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-0.40487  -0.16209   0.00298   0.15201   0.96271  

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -2.61381    0.27229  -9.599 2.53e-09 ***
Monetary     0.07482    0.01266   5.908 6.04e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasibinomial family taken to be 0.06566381)

    Null deviance: 4.2966  on 23  degrees of freedom
Residual deviance: 1.8027  on 22  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 5

> summary(m.glm)

Call:

glm(formula = Percentage ~ Monetary, family = quasibinomial(link = "logit"),

data = p)

Deviance Residuals:

Min 1Q Median 3Q Max

-0.40487 -0.16209 0.00298 0.15201 0.96271

Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) -2.61381 0.27229 -9.599 2.53e-09 ***

Monetary 0.07482 0.01266 5.908 6.04e-06 ***

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasibinomial family taken to be 0.06566381)

Null deviance: 4.2966 on 23 degrees of freedom

Residual deviance: 1.8027 on 22 degrees of freedom

AIC: NA

Number of Fisher Scoring iterations: 5

RecencyとFrequencyによるロジスティック回帰

model <- glm(Buy~Recency+Frequency,family=quasibinomial(link='logit'),data=train)
pred <- predict(model,data.frame(Recency=c(0),Frequency=c(1)),type='response')

1 2	model <- glm(Buy~Recency+Frequency,family=quasibinomial(link='logit'),data=train) pred <- predict(model,data.frame(Recency=c(0),Frequency=c(1)),type='response')

> pred
        1 
0.2579282

> pred

0.2579282

顧客生涯価値の計算

推定したロジットを用いて、生涯価値を計算します。

# 3期間に関する顧客生涯価値の計算（Recency=0、Frequency=1、平均収益=100ドル、割引率=2%）
v<-getCLV(0,1,100,1,0,3,0.02,model)

1 2	# 3期間に関する顧客生涯価値の計算（Recency=0、Frequency=1、平均収益=100ドル、割引率=2%） v<-getCLV(0,1,100,1,0,3,0.02,model)

> v
[1] 63.91906

1 2	> v [1] 63.91906

それではさっそく、1998年5月〜6月のデータを用いて、今回推定した顧客生涯価値が妥当なのかどうかを確かめたいと思います。

#トレーニングデータをもとにCLVを計算する
train_pred <- cbind(train,lifetimevalue=matrix(,nrow(train),1))

pb <- txtProgressBar(min = 1, max = nrow(train_pred), style = 3)

for(i in 1:nrow(train_pred)){
  train_pred$lifetimevalue[i] <- getCLV(train_pred$Recency[i],train_pred$Frequency[i],100,1,0,3,0.02,model)
  setTxtProgressBar(pb, i) 
}

#1998年5月1日〜6月30日の間の取引データからIDと購入金額を集計する。
startDate_test <- as.Date("19980501","%Y%m%d")
endDate_test <- as.Date("19980630","%Y%m%d")
test <- getDataFrame(df,startDate_test,endDate_test)
test_monetary <- dplyr::select(test,ID,Monetary)

#予測したCLVと実際の取引金額のデータをマージする。
train_pred_2 <- data.frame(train_pred[,colnames(train_pred) %in% c("ID","lifetimevalue")])
train_pred_2 <- merge(train_pred,test_monetary,by = "ID",all=T)
train_pred_2$Monetary.y[is.na(train_pred_2$Monetary.y)] <- 0

#トレーニングデータをもとにCLVを計算する

train_pred <- cbind(train,lifetimevalue=matrix(,nrow(train),1))

pb <- txtProgressBar(min = 1, max = nrow(train_pred), style = 3)

for(i in 1:nrow(train_pred)){

train_pred$lifetimevalue[i] <- getCLV(train_pred$Recency[i],train_pred$Frequency[i],100,1,0,3,0.02,model)

setTxtProgressBar(pb, i)

}

#1998年5月1日〜6月30日の間の取引データからIDと購入金額を集計する。

startDate_test <- as.Date("19980501","%Y%m%d")

endDate_test <- as.Date("19980630","%Y%m%d")

test <- getDataFrame(df,startDate_test,endDate_test)

test_monetary <- dplyr::select(test,ID,Monetary)

#予測したCLVと実際の取引金額のデータをマージする。

train_pred_2 <- data.frame(train_pred[,colnames(train_pred) %in% c("ID","lifetimevalue")])

train_pred_2 <- merge(train_pred,test_monetary,by = "ID",all=T)

train_pred_2$Monetary.y[is.na(train_pred_2$Monetary.y)] <- 0

予測したCLVと実際の取引金額データを散布図で描き、回帰線を引く。

library(ggplot2)
g <- ggplot(train_pred_2,aes (x = lifetimevalue,y = Monetary.y))
g <- g +  geom_point(shape = 20,size = 0.8,na.rm = TRUE)
g <- g + geom_smooth(method = glm) + xlab("lifetimevalue") + ylab("Monetary")
plot(g)

library(ggplot2)

g <- ggplot(train_pred_2,aes (x = lifetimevalue,y = Monetary.y))

g <- g + geom_point(shape = 20,size = 0.8,na.rm = TRUE)

g <- g + geom_smooth(method = glm) + xlab("lifetimevalue") + ylab("Monetary")

plot(g)

CLVが上がれば、1998年5月1日〜6月30日の間（未来）の取引金額が増すような傾向が出ています。

参考情報

RFM Customer Analysis with R Language
Calculating Customer Lifetime Value with Recency, Frequency, and Monetary (RFM)

OpenCV&Pythonで画像の類似度を計算させる〜イケメンの顔比較

・動機
・やりたいこと
・準備
・類似度の計算
・実行コード
・実行結果
・おまけ

動機

画像系の技術にあまり関心が無かったのですが、とある知人が福士蒼汰のような雰囲気の男性が好みであると発言されたことを発端に、福士蒼汰に最も顔の近い知人を見つけるというプライベートなミッションを仰せつかりました。
そこで、まずは最も楽だろうと思われる、画像間の類似度を計算する方法について調べてみました。顔のパーツを検知して、パーツ同士で比較するなどのレベルでは無いことをご了承下さい。ちなみに、比較画像は国民的アイドルである嵐のメンバーの画像としました。今回の実践で、与えられた画像の中で、福士蒼汰に最も近そうな嵐のメンバーの写真がわかることになります。

やりたいこと

画像間の類似度の計算

準備

・Python2.7
・OpenCV

OSX環境における準備にあたっては以下の情報が参考になりました。
PythonでOpenCVを使う
 Mac OS X で OpenCV 3 + Python 2/3 の開発環境を整備する方法

類似度の計算

こちらのブログにある計算手法とコードを使いました。
How-To: Python Compare Two Images

紹介されているコードで以下の評価指標が計算できます。
・Mean Squared Error (MSE)
・Structural Similarity Index (SSIM)・・・0〜1の値を取ります。

指標について、詳しくはこちらの論文に書かれています。
Image Quality Assessment: From Error Visibility to Structural Similarity

実行コード

紹介されていたコードは画像サイズが同じでないと計算ができなかったので、まずは画像のサイズを整えるためのコードが以下のようになります。コードはこちらのものを使いました。
Python – アンチエイリアスで写真をキレイに縮小

# -*- coding: utf-8 -*-
import commands
from PIL import Image
import re

# 縮小する際の画像の高さピクセル
PHOTO_HEIGHT = 400

# 画像があるフォルダのフルパス
BASE_DIR = "/任意のフォルダのパス/"

# 画像の正規表現名
PHOTO_REGEX = r"pic00.*.[jpg|jpeg]"

# リサイズ後の画像の接頭語
PHOTO_RESIZE_PREFIX = "r_"

def main():
    # 画像フルパスを取得
    _cmd = "cd {} && ls".format(BASE_DIR)
    l = commands.getoutput(_cmd)
    l = l.split("\n")
    l = [_l for _l in l if re.match(PHOTO_REGEX, _l)]

    # 出力用のフォルダを生成
    commands.getoutput("mkdir {}/output".format(BASE_DIR))

    # 既存ファイルを readモードで読み込み
    for _l in l:
        before_path = '{}/{}'.format(BASE_DIR, _l)
        filename = '{}{}'.format(PHOTO_RESIZE_PREFIX, _l)
        after_path = '{}/output/{}'.format(BASE_DIR, filename)
        resize(before_path, after_path, filename=_l)  # 縮小


def resize(before, after, height=PHOTO_HEIGHT, filename="", aa_enable=True):
    """
    画像をリサイズする
    :param str before: 元画像ファイルパス
    :param str after: リサイズ後の画像ファイルパス
    :param int height: リサイズ後の画像の高さ
    :param bool aa_enable: アンチエイリアスを有効にするか
    :return:
    """
    # 画像をreadonlyで開く
    img = Image.open(before, 'r')
    # リサイズ後の画像ピクセルを計算
    before_x, before_y = img.size[0], img.size[1]
    x = int(round(float(height / float(before_y) * float(before_x))))
    y = height
    resize_img = img
    if aa_enable:
        # アンチエイリアスありで縮小
        resize_img.thumbnail((x, y), Image.ANTIALIAS)
    else:
        # アンチエイリアスなしで縮小
        resize_img = resize_img.resize((x, y))

    # リサイズ後の画像を保存
    resize_img.save(after, 'jpeg', quality=100)
    print "RESIZED!:{}[{}x{}] --> {}x{}".format(filename, before_x, before_y, x, y)


# 実行
main()

# -*- coding: utf-8 -*-

import commands

from PIL import Image

import re

# 縮小する際の画像の高さピクセル

PHOTO_HEIGHT = 400

# 画像があるフォルダのフルパス

BASE_DIR = "/任意のフォルダのパス/"

# 画像の正規表現名

PHOTO_REGEX = r"pic00.*.[jpg|jpeg]"

# リサイズ後の画像の接頭語

PHOTO_RESIZE_PREFIX = "r_"

def main():

# 画像フルパスを取得

_cmd = "cd {} && ls".format(BASE_DIR)

l = commands.getoutput(_cmd)

l = l.split("\n")

l = [_l for _l in l if re.match(PHOTO_REGEX, _l)]

# 出力用のフォルダを生成

commands.getoutput("mkdir {}/output".format(BASE_DIR))

# 既存ファイルを readモードで読み込み

for _l in l:

before_path = '{}/{}'.format(BASE_DIR, _l)

filename = '{}{}'.format(PHOTO_RESIZE_PREFIX, _l)

after_path = '{}/output/{}'.format(BASE_DIR, filename)

resize(before_path, after_path, filename=_l) # 縮小

def resize(before, after, height=PHOTO_HEIGHT, filename="", aa_enable=True):

"""

画像をリサイズする

:param str before: 元画像ファイルパス

:param str after: リサイズ後の画像ファイルパス

:param int height: リサイズ後の画像の高さ

:param bool aa_enable: アンチエイリアスを有効にするか

:return:

"""

# 画像をreadonlyで開く

img = Image.open(before, 'r')

# リサイズ後の画像ピクセルを計算

before_x, before_y = img.size[0], img.size[1]

x = int(round(float(height / float(before_y) * float(before_x))))

y = height

resize_img = img

if aa_enable:

# アンチエイリアスありで縮小

resize_img.thumbnail((x, y), Image.ANTIALIAS)

else:

# アンチエイリアスなしで縮小

resize_img = resize_img.resize((x, y))

# リサイズ後の画像を保存

resize_img.save(after, 'jpeg', quality=100)

print "RESIZED!:{}[{}x{}] --> {}x{}".format(filename, before_x, before_y, x, y)

# 実行

main()

これだけではサイズが同一にならなかったので、追加で以下のコードを実行しました。

from PIL import Image
img = Image.open('output/r_pic005.jpg') # image extension *.png,*.jpg
new_width  = 300
new_height = 400
img = img.resize((new_width, new_height), Image.ANTIALIAS)
img.save('resize/re_pic005.jpg')

from PIL import Image

img = Image.open('output/r_pic005.jpg') # image extension *.png,*.jpg

new_width = 300

new_height = 400

img = img.resize((new_width, new_height), Image.ANTIALIAS)

img.save('resize/re_pic005.jpg')

以下は、類似度計算の実行用コードとなります。

# import the necessary packages
from skimage.measure import structural_similarity as ssim
import matplotlib.pyplot as plt
import numpy as np
import cv2
%matplotlib inline

# import the necessary packages

from skimage.measure import structural_similarity as ssim

import matplotlib.pyplot as plt

import numpy as np

import cv2

%matplotlib inline

def mse(imageA, imageB):
    # the 'Mean Squared Error' between the two images is the
    # sum of the squared difference between the two images;
    # NOTE: the two images must have the same dimension
    err = np.sum((imageA.astype("float") - imageB.astype("float")) ** 2)
    err /= float(imageA.shape[0] * imageA.shape[1])
    
    # return the MSE, the lower the error, the more "similar"
    # the two images are
    return err
 
def compare_images(imageA, imageB, title):
    # compute the mean squared error and structural similarity
    # index for the images
    m = mse(imageA, imageB)
    s = ssim(imageA, imageB)
 
    # setup the figure
    fig = plt.figure(title)
    plt.suptitle("MSE: %.2f, SSIM: %.2f" % (m, s))
 
    # show first image
    ax = fig.add_subplot(1, 2, 1)
    plt.imshow(imageA, cmap = plt.cm.gray)
    plt.axis("off")
 
    # show the second image
    ax = fig.add_subplot(1, 2, 2)
    plt.imshow(imageB, cmap = plt.cm.gray)
    plt.axis("off")
 
    # show the images
    plt.show()

def mse(imageA, imageB):

# the 'Mean Squared Error' between the two images is the

# sum of the squared difference between the two images;

# NOTE: the two images must have the same dimension

err = np.sum((imageA.astype("float") - imageB.astype("float")) ** 2)

err /= float(imageA.shape[0] * imageA.shape[1])

# return the MSE, the lower the error, the more "similar"

# the two images are

return err

def compare_images(imageA, imageB, title):

# compute the mean squared error and structural similarity

# index for the images

m = mse(imageA, imageB)

s = ssim(imageA, imageB)

# setup the figure

fig = plt.figure(title)

plt.suptitle("MSE: %.2f, SSIM: %.2f" % (m, s))

# show first image

ax = fig.add_subplot(1, 2, 1)

plt.imshow(imageA, cmap = plt.cm.gray)

plt.axis("off")

# show the second image

ax = fig.add_subplot(1, 2, 2)

plt.imshow(imageB, cmap = plt.cm.gray)

plt.axis("off")

# show the images

plt.show()

# load the images -- the original, the original + contrast,
# and the original + photoshop
original = cv2.imread("resize/re_pic006.jpg")
contrast = cv2.imread("resize/re_pic005.jpg")
shopped = cv2.imread("resize/re_pic003.jpg")
    
# convert the images to grayscale
original = cv2.cvtColor(original, cv2.COLOR_BGR2GRAY)
contrast = cv2.cvtColor(contrast, cv2.COLOR_BGR2GRAY)
shopped = cv2.cvtColor(shopped, cv2.COLOR_BGR2GRAY)

# load the images -- the original, the original + contrast,

# and the original + photoshop

original = cv2.imread("resize/re_pic006.jpg")

contrast = cv2.imread("resize/re_pic005.jpg")

shopped = cv2.imread("resize/re_pic003.jpg")

# convert the images to grayscale

original = cv2.cvtColor(original, cv2.COLOR_BGR2GRAY)

contrast = cv2.cvtColor(contrast, cv2.COLOR_BGR2GRAY)

shopped = cv2.cvtColor(shopped, cv2.COLOR_BGR2GRAY)

# initialize the figure
fig = plt.figure("Images")
images = ("Original", original), ("Contrast", contrast), ("Photoshopped", shopped)
 
# loop over the images
for (i, (name, image)) in enumerate(images):
    # show the image
    ax = fig.add_subplot(1, 3, i + 1)
    ax.set_title(name)
    plt.imshow(image, cmap = plt.cm.gray)
    plt.axis("off")
 
# show the figure
plt.show()
 
# compare the images
compare_images(original, original, "Original vs. Original")
compare_images(original, contrast, "Original vs. Contrast")
compare_images(original, shopped, "Original vs. Photoshopped")