Create a gist now

Instantly share code, notes, and snippets.

Embed
What would you like to do?
寒い冬の過ごし方。 #vgadvent2013 〜アイスクリーム統計学をSQLで書いてみた〜

寒い冬の過ごし方。

この記事はVOYAGE GROUP エンジニアブログ : Advent Calendar 2013 の5日目の記事になります!

すっかり寒くなりましたね。 こんな寒い冬はこたつに入ってアイスクリームを食べるのが一番です。

12月っぽいアイス(お菓子?)の家

アイスクリームと言えば、**アイスクリーム統計学**ですね。 店長さんが最高気温による予測客数を知りたいとのことなので回帰係数を求めてみます。

アイスクリーム統計学 4. 単回帰 ※外部サイトです。

このウェブサイトではExcelを使って計算しているのですが、 あいにくExcelのインストールされたPCを持っていなかったのでSQLで求めてみます。

ちなみに弊社では非エンジニアもSQLを書いてDBに問い合わせ、分析を行っているので 仮に回帰分析の説明を求められてもExcelでやるよりずっとスムーズに伝わるはずです!?

計測データのテーブル作成

CREATE TABLE visitorlog 
(
 id int auto_increment primary key, 
 temp int not null, 
 visits int not null
);
INSERT INTO visitorlog
(temp, visits)
VALUES
(33,382),
(33,324),
(34,338),
(34,317),
(35,341),
(35,360),
(34,339),
(32,329),
(28,218),
(35,402),
(33,342),
(28,205),
(32,368),
(25,196),
(28,304),
(30,294),
(29,275),
(32,336),
(34,384),
(35,368);

相関係数(cor)

ページには

相関係数=偏差積の平均/(xの標準偏差*yの標準偏差)

とあります。※ xの標準偏差 = √(x-xの平均)^2/xデータ数

標準偏差(stdev)

SQLで書くとこんな感じになります。

select 
  sqrt(sum(pow(temp-avg_temp,2))/cnt) stdev_temp,
  sqrt(sum(pow(visits-avg_visits,2))/cnt) stdev_visits 
from visitorlog v ,(
  select avg(temp) avg_temp, avg(visits) avg_visits,count(*)cnt from visitorlog
)a;

偏差積の平均

こうなりますね。

select avg(b.dev_prod) avg_dev_prod from(
  select 
    (temp-avg_temp) * (visits-avg_visits) dev_prod
  from visitorlog v ,(
    select avg(temp) avg_temp, avg(visits) avg_visits from visitorlog
  )a
)b

なので、それらを組み合わせれば、

select avg_dev_prod/(stdev_temp*stdev_visits) cor from (
select avg(b.dev_prod) avg_dev_prod from(
  select 
    (temp-avg_temp) * (visits-avg_visits) dev_prod
  from visitorlog v ,(
    select avg(temp) avg_temp, avg(visits) avg_visits from visitorlog
  )a
)b
)c,(
  select 
  sqrt(sum(pow(temp-avg_temp,2))/cnt) stdev_temp,
  sqrt(sum(pow(visits-avg_visits,2))/cnt) stdev_visits 
from visitorlog v ,(
  select avg(temp) avg_temp, avg(visits) avg_visits, count(*) cnt 
  from visitorlog
)a
)d

やったね。ちなみにこの結果は0.87。最高気温と客数には強い相関があるといって良さそうです。

切片(intersect)、傾き(slope)

相関係数が出てしまえば、あとは簡単。

回帰直線の傾き=相関係数*((yの標準偏差)/(xの標準偏差))

y切片=yの平均-(傾き*xの平均)

これを計算するには、

select slope, avg_visits-slope*avg_temp intercept from(
select
  cor*(stdev_visits/stdev_temp) slope,
  avg_temp,
  avg_visits    
from(
  select 
    avg_dev_prod/(stdev_temp*stdev_visits) cor,
    stdev_visits,
    stdev_temp,
    avg_temp,
    avg_visits
  from (
    select avg(b.dev_prod) avg_dev_prod from(
      select 
        (temp-avg_temp) * (visits-avg_visits) dev_prod,
        avg_temp,
        avg_visits
      from visitorlog v ,(
        select 
          avg(temp) avg_temp, 
          avg(visits) avg_visits 
        from visitorlog
      )a
    )b
)c,
(
select 
  sqrt(sum(pow(temp-avg_temp,2))/cnt) stdev_temp,
  sqrt(sum(pow(visits-avg_visits,2))/cnt) stdev_visits,
  avg_temp,
  avg_visits
from visitorlog v ,(
  select 
    avg(temp) avg_temp, 
    avg(visits) avg_visits, 
    count(*) cnt 
  from visitorlog
)a
)d
)e
)f

傾きと切片がでましたね。RとかSPSSとか難しいツールを使わなくても ちょっと考えれば使い慣れているSQLでなんとかなりました。

今回のSQLはSQL Fiddleを利用して作成しました。初めて使いましたが便利です。 http://sqlfiddle.com/#!2/db3cf/1

所感

普段便利ツールばっかり使ってて意識してなかったことを ちょっと手間かけて考えてみたのですが、だいぶ理解が深まった気がしました。

来年度から社会人3年目になりますが、 まだまだ思うように仕事をこなせることばかりではなく悔しいことが多々あります。 しかし弊社の先輩や同期、後輩は優秀な人達ばかりで 仕事の仕方、着眼点、思考のプロセス、意識していることなど学べることばかりです。 時間をかけてでも自分も高めていきたいと思います。

明日は、最近お仕事や社内勉強会でご一緒させてもらっている @akiyah さんです。お楽しみに!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment