【MATLAB】当ブログの検索流入件数の傾向と予測

はてなブログでProにアップグレードして独自ドメインに移行後約250日が経ちました

検索流入もだんだんと上がってきたので，今後の伸び方を大雑把に予想してみます

単純に線形近似してみる

独自ドメインへ移行後の日数，累計記事数，検索流入件数は図1のようになりました

また，検索流入件数を線形近似した結果とその残差も載せました

図1．線形近似した結果

問題点

検索流入件数を見ると，130日目あたりまでは変化が小さく徐々に増加していた一方で，130日目より後では急激に増加し大きく変化しながら全体的に増加しているように見えます

これはブログ界でよく知られている現象で，約半年後にGoogleに信頼できるサイトであると判断されることによるものだと思われます

したがって，Googleに様子見されている129日目までと，130日目以降の検索流入件数は分けて考えることにしました

130日目より前と後で別々に線形近似

130日目より前と後で別々に線形近似した結果．図2のようになりました

それぞれの残差プロットを図3と図4に示しました

図2．線形近似

図3．残差プロット（129日目まで）

図4．線形近似（130日目以降）

モデル式$ y=p1*x+p2 $について，フィッティングの結果下表のようになりました

	p1（95% の信頼限界）	p2 （95% の信頼限界）	R二乗値	修正済みR二乗値	自由度	RMSE
129日目まで	0.1411 (0.1094, 0.1729)	21.65 (19.27, 24.03)	0.3786	0.3737	127	6.7856
130日目以降	0.2306 (0.1222, 0.339)	0.2306 (0.1222, 0.339)	0.1625	0.1534	92	14.3587

129日目まではグラフを見ても値を見てもまぁまぁ当てはまってるなぁと思います

130日目以降の線形近似の残差プロットで約50日周期の変化が見られますが，今回は大雑把な傾向を知りたいので今回はノイズとして扱うことにします（つまり，約50日周期の残差の変動は気にしない．）

しかしこの変動が大きく日数が足りてないため，このフィッティングが果たして適切なのかイマイチ分かりません

数式的な考察

検索流入件数を数式的に表してみます

なお，以下の「○○日目」とは，全て独自ドメイン移行後の日数を意味します

まず，ある記事$ i $を$ n_i $日目に公開したときの，$ n $日目の検索流入件数$ f_i(n) $を考えます

検索流入は公開してすぐに増えるわけではなく，ある程度日数が経ってから検索上位にヒットするようになるので，$ n_i+\Delta n_i $日目以降に一日あたり平均$ K_i $の検索流入があるとすると，$ f_i(n) $は式(1)と表せます

現実には検索流入件数は日によってかなり変動しますが，ここで考えている記事数は150記事以上を考えているため，各記事の検索流入の変動を合わせれば全体の変動は小さくなります

そのため一つの記事の検索流入が一定であるとして考えても良いだろうということです

なお，式中の$ u(x) $はステップ関数で，$ x\geq 0 $のとき$ 1 $, $ x<0 $のとき$ 0 $をとります

\begin{aligned} f_i(n)&= \left\{ \begin{aligned} &K_i& (n\geq n_i+\Delta n_i )\\ &0&(n< n_i+\Delta n_i ) \end{aligned} \right.\\ &=K_iu(n-(n_i+\Delta n_i)) \end{aligned}\tag{1}

さらに$ \Delta n_i,\,K_i $も同様に，全記事について$ \Delta n_i = \Delta n,\,K_i=K $ と考えしまえば，

\begin{aligned} f_i(n)&=Ku(n-(n_i+\Delta n)) \end{aligned}\tag{2}

すると，$ n $日目の全検索流入件数$ \displaystyle f(n)=\sum_i f_i(n) $ は，

\begin{aligned} f(n)&=K\sum_i u(n-(n_i+\Delta n))\\ &=K\sum_i u((n-\Delta n)-n_i) \end{aligned}\tag{3}

式中の$ u((n-\Delta n)-n_i) $ は，$ n-\Delta n\geq n_i $のときに$ 1 $，$ n-\Delta n< n_i $のときに$ 0 $となることから，「$ \Delta n $日前に記事$ i $が存在していたかどうか」を表します

したがって，$ \displaystyle \sum_i u((n-\Delta n)-n_i) $は「$ \Delta n $日前の時点の累積記事数」を表しており，$ n $日目の検索流入件数は$ \Delta n $日前の時点での累積記事数に比例すると言えます

一日あたり平均$ p $記事のペース書いていて0日目で$ q $記事既に存在していた時，$ n $日目の累計記事数は$ pn+q $と表されますから，$ n $日目の検索流入件数$ f(n) $は式 (3) より，

\begin{aligned} f(n)&=K\sum_i \left( p(n-\Delta n)+q\right) \\ &= K\left(\frac{1}{2}pn(n+1)-pn\Delta n+qn\right) \\ &=K\left(\frac{1}{2}pn ^2+\left(\frac{1}{2}p-p\Delta n+q\right)n\right) \\ \end{aligned}\tag{4}

上の結果から，あるペースで記事を新たに公開していると，検索流入件数には日数の二乗に比例する項が現れることが分かりました

二次近似

130日目以降のデータについて，今度は二次近似をした結果，近似曲線は図4のようになりました

約50日周期のノイズが大きすぎてそれを拾ってしまったように見えます

図4. 二次近似の結果

	p1（95% の信頼限界）	p2（95% の信頼限界）	p3 （95% の信頼限界）	R二乗値	修正済みR二乗値	自由度	RMSE
130日目以降	0.003805 (-0.0006172, 0.008227)	-1.113 (-2.677, 0.4521)	147.6 (11.74, 283.5)	0.1885	0.1707	91	14.2110

式 (4) が成り立つとして，現在の一日あたりの公開ペースの平均$ p\approx 0.5 $, 一記事あたりの一日の平均検索流入数$ K\approx 0.4 $なので，二次の項の係数$ \frac{pK}{2} \approx 0.1 $ですから，二次の成分を見出すには日数が足りなすぎるのかもしれません

数式モデルへのデータ適合問題として解く

今度は式 (4) の式の形で関数を作成し，$ p, q, \Delta n, K $を決定させます

この方法なら，各パラメータの範囲を指定できるので，良い結果が得られるかもしれません

ただ，式 (4) は独自ドメイン移行直後から成り立つものではありません

さらに，130日目より前ではGoogleに様子見されていて，130日目より後と完全に傾向が異なるため，130日目以降の傾向が130日目より前にもあったと仮定したときの0日目の検索流入は0から大きくずれるでしょう

したがって，モデル式を式 (5) のように切片$ a $を追加して再定義した上でデータ適合を行いました．（以下，$ \frac{1}{2}p \ll q $より$ \frac{1}{2}p $の項は無視します．）

\begin{aligned} f(n)&=K\left(\frac{1}{2}pn ^ 2+\left(-p\Delta n+q\right)n\right)+a \end{aligned}\tag{5}

（パラメータの数が多くて解が1つに定まらない気がしますが，今回はそれっぽい可能解が見つかれば良いことにします）

データ適合の結果

実行結果

検索流入件数（130日目以降）を線形近似，データ適合した結果

パラメータの結果は（$ K, p, \Delta n, q, a $の順）

    0.1467    0.0710  604.0735   32.1189  185.9286

1000日目までのグラフを見ると，データ適合はちょっと上がりすぎな気はしますが，線形近似よりかは現実味があるように見えます（筆者の欲に強く影響された主観ですが）

300日目までのグラフを見ると，線形近似は良さそうですが，データ適合の方は130日目〜140日目の一時的な減少に影響されて，最小値が140日目付近にあるように適合されてしまいました

モデル式の定義では，各記事の一日あたりの検索流入件数は一定であり，公開する記事が減る（つまり，記事を削除する）ことも想定していないため，nが正であれば単調増加するはずです

改善

次に，nが正のときに単調増加するようにモデル式に制約条件を加えます

しかし，MALABではデータ適合のパラメータに制約式を与える方法が見つからないので，少し工夫をしました

モデル式の二次関数について，最小値をとる$ n $を$ n^* $ とすると，

\begin{aligned} n^*=\frac{p\Delta n-q}{p}\tag{5} \end{aligned}

ここで，式 (6) のように同値変形します

\begin{aligned} n\text{が非正のとき最小値をとる} &\Leftrightarrow n^*\text{が非正}&(\because\text{ 二次の係数が正なので下に凸})\\ &\Leftrightarrow \frac{p\Delta n-q}{p}\text{が非正}\\ &\Leftrightarrow p\Delta n-q\text{が非正}&(\because p>0)\\ \end{aligned}\tag{6}

以上から，$ T=p\Delta n-q $とおいて$ q=p\Delta n-T $を代入し，パラメータ$ T $の最大値を0に指定してデータ適合を行えば良いことが分かりました

式にすると，

\begin{aligned} f(n)&=K\left(\frac{1}{2}pn ^2-pn\Delta n+(p\Delta n-T)n\right)+a\\ &=K\left(\frac{1}{2}pn ^2-Tn\right)+a\\ &(K,p,\,\Delta n,\,a>0,\,T\leq0) \end{aligned}\tag{6}

この結果，下図のようになりました

検索流入件数（130日目以降）を線形近似，データ適合（改善後）した結果

各パラメータの値は（$ K, p, T, a $の順）

    0.0504    0.0306   -0.0045   48.4653

改善前や線形近似よりかなり現実的な曲線となりました

どのパラメータも制約条件の端の値にはならず，制約条件の中の局所的最適解に落ち着いているようです

$ K, p, \Delta n $は現在の当ブログの値（$ K\approx 0.4,\,p\approx 0.5 $）とオーダーが一つずれていますが，まだ日数が少ないですし，これ以上は問わないことにします

さて，上のパラメータの値から$ q $の値は$ q = p\Delta n-T $で出るかと思いきや，式 (6) で$ \Delta n $は消えてしまったため$ \Delta n $の値が定まらず，$ q $と$ \Delta n $の組み合わせは無限に存在することになってしまいました

指定したパラメータの数が多すぎたみたいです

全てのパラメータの値が定まらなくてもそれっぽい曲線は導けたので，今回はこれで充分ってことにします

モデリングあたりの知識が皆無なので近いうち勉強したいな...

MATLABのコードの完成版

Google Analyticsでデータをcsv形式でエクスポートした後，日数の列を追加し，MATLABにaccessというtable型の変数にインポートして使います

130日目以降について，線形近似とデータ適合の２つの解析結果を，300日目までを拡大して表示するグラフと1000日目まで俯瞰するグラフの2つを出力します

さらに，線形近似の統計データと，データ適合によるパラメータの適合結果をコマンドウィンドウに表示します

コードを見る（クリックして下さい）

検索流入の解析（130日目以降，線形近似とデータ適合）

おわりに

今回は線形回帰とデータ適合の2つの手法を用いて検索流入件数の大まかな傾向を把握しましたが，未来の検索流入件数を予想するにはまだ不十分でした

今後は定期的に新しいデータを追加して解析してみて様子見することにします

（本ブログのaboutページに掲載しています）

また，毎月の成果報告でも解析結果の報告をする予定です

二乗に比例してくれると良いなぁ...

本記事執筆中，試行錯誤がなかなか面白いと感じる一方で，統計の知識が圧倒的に足りてないことに気づきました

統計検定2級は持っているものの，それは分布の特徴が分かる程度の基礎知識で，モデリングあたりについては全く知りません

そのあたりの問題が準一級以上で出るそうなので，近いうち受けようかと思います