オブジェクトムービーのお引越し  オブジェクトムービー

このブログで作成&紹介していたオブジェクトムービーは、新しいブログ『ぐりぐり写真を作ろう』にお引越しいたします。
新しい作品は『ぐりぐり写真を作ろう』をみてくださいね。古い作品も順次移動していきます。
3

2015/12/15

スターバックスのコーヒーを1000杯飲んで気がついたこと 【ポアンカレのパン(1)】 #vgadvent2015  プログラミング

※この記事はVOYAGE GROUP Advent Canlendarの14日目の記事です。

※※かつ、この記事は【ポアンカレのパン】シリーズの第1弾です。(アドベントカレンダーなのに別のシリーズのスタート地点にしてしまいました)

みなさんはポアンカレのパンの話をご存知でしょうか。

数学者ポアンカレがパン屋で1kgのパンを毎日買っていたが、1年間毎日計測していたら実はパンの平均は950gで、パン屋のズルを見抜いたと言う統計学の本に出てくる逸話です。

これを知ったとき私はいくつかのことを思いました。
・統計ってすごい!ポアンカレすごい!
・自分でも同じようなことやってみたい!
・1kgのパンって大きすぎじゃない?いまだと売っている食パンは1斤340gくらいだよ?毎日食パン3斤食べるか?
・他にも細かいこといくつか

この【ポアンカレのパン】シリーズでこれらの気持ちを解決して行くとして、まず最初の今回は『自分でも同じようなことやってみたい!』についてです。

パンはそんなに買わないので、スターバックスのコーヒーでやってみることにしました。思い立ったのが2013年の12月で、平日の朝と夕方にコーヒーを買って、その重さを計り続けました。2年かかって2015年12月に目標の1000杯を計測終わりました。長かった。

クリックすると元のサイズで表示します
こんな風に、タンブラーに入れてもらったコーヒーをタンブラーごと測って、タンブラーの重さ189gを引いてコーヒーの重さとします。

トールサイズのタンブラーにショートサイズを入れてもらうという方針にしています。タンブラーにいっぱいに入れる量だと日々の散らばりが少ないだろうと思ったので、ゆとりのあるタンブラーに対して途中まで入れるような状況を作ったのです。
紙コップでショートサイズコーヒーを買うと、なみなみ入っている状態で約200gでしたので、200gを基準にしてどんな風な分布をしているのか、調べてみましょう。

だんだんサンプルサイズが大きくなって行くようなヒストグラムを作ってみました。
クリックすると元のサイズで表示します

ですが結構分かりづらいので、素直にサンプルサイズ1000の時の静止画も用意しました。
クリックすると元のサイズで表示します
赤い線は、サンプルから出した平均と分散を元に正規分布曲線を書いたものです。

ポアンカレの時とは逆で、200gよりも多めに入れてもらっていることがわかります!平均247.2gですね!
1000gのパンを5%少なく950gで売っていたポアンカレのパン屋に比べて、200gのコーヒーを20%以上多く240g入れてくれていた私のスターバックスはなんと素晴らしいのでしょう!!

ところで、スターバックスの公称のショートサイズコーヒーの量はどれくらいなのでしょうか。
http://www.starbucks.co.jp/howto/store/order.html
ここに書いてありました。

クリックすると元のサイズで表示します
ショート
Short(240ml)
各サイズの容量は、目安です。


。。。
なんだ、公称が240か。
まあでも、平均で7.2gもおおく入れてくれていたので、感謝です!!

【ポアンカレのパン】シリーズは、気が向いたら続きを書きます。

明日のVOYAGE GROUP Advent Canlendarは、@yukimineです。お楽しみに!

2

2015/4/6

ABテストの失敗例:平均を安易に使う  

私が経験したABテストの失敗例を書いていきます。

【失敗例1】
平均を安易に使う

【設定】宝くじ
全部で100000枚
1等 100000円が10枚
2等 10000円が100枚
3等 1000円が1000枚
4等 100円が10000枚
はずれ (100000 - 11110)枚

A: 赤いパンツを履いてもらった人10000人に1人1枚ずつ購入してもらう
B: パンツの色の指定なしに10000人に1人1枚ずつ購入してもらう
として、ABの平均当たり金額を比べてみようとします。

もちろん、赤いパンツを履いても当選確率は変わらないという前提ですよ。

Rでシミュレーションをしてみます。


平均あたり金額は、Aが40.86円、Bが30.96円になりました。
Bに対してAは30%以上差がついているし、10000人という十分過ぎる人数でテストしているので、これはAが有意に差があるにちがいない!と思ってしまいます。

でもこれは間違いなんですね。
シミュレーションを10000回繰り替えしてみます。



クリックすると元のサイズで表示します

ヒストグラムを見ると25円くらいから60円くらいまで分布が広がっているように見えます。
計算で5パーセンタイルと95パーセンタイルを出してみると、26.65円と59.10円なので、シミュレーション1回毎に90%の確率でこの範囲に入ることがわかります。

ということは、さっきのABの範囲はこの中に収まるわけで、ぜんぜん有意な差ではないんですね。
Bに対してAは30%以上差がついているのに有意な差ではないというのは不思議ですが、それは元の分布の標準偏差が大きいからです。

今回は説明しやすさのために宝くじを例に出しましたが、低額なものも高額なものも扱っているサイトでの「購入金額」でも同様に標準偏差が大きくなって、ABテストがやりにくくなります。
こういう場合は「平均当たり金額」ではなくて、「当たった人数」や「3等以上が当たった人数」などのような、普通のABテストに近い設定にしたほうがいいです。
0

2014/12/15

ジバニャン方程式の作り方 #vgadvent2014  



※この記事はVOYAGE GROUP エンジニアブログ : Advent Calendar 2014の16日目の記事として書かれています。

11月にジバニャン方程式発見のツイートをして9000以上のRTをいただきました。




このエントリーではこのジバニャン方程式をどうやって作ったのかを紹介します。
実際には試行錯誤をしながら道具や環境を整えて行ったのでこの順番に作ったわけでは無いのですが、説明しやすい順番で書きますね。

0. 自分のモチベーションと目的を確認する。



仕事じゃないんだからモチベーションとか目的とか言わなくていいでしょ、と思うでしょう。私もそう思います!
でも、自分の趣味とはいえ何日も何週間もやっているとそうも言っていられなくなるんですよね。
ジバニャン方程式は、構想を含めなければ10月22日から手をつけたようで、発表が11月12日なので、3週間くらいやっていたひとりプロジェクトでした。
3週間!!!
ふりかえって見ると長いですね。

途中で何度か目的を見失って、このままだと完成しないなーと感じることがありました。
そういうときはたいてい、こだわりすぎて作業が前進しないときなので、何をしたかったのか自分に問いかけてこだわりを割り切ったりしました。

【やりたかったこと】
・ジバニャンを方程式で書く(以前やったアンパンマン方程式の続編)
・その方程式はTeXでかっこいい数式にする
・方程式は1個にする。(別で定義した式を使わない)

【割り切って捨てたこと】
・かっこいいDSL
・正確なジバニャンの描画(口とか、目のまわりとか)
・maxとminは使いたくなかったけど、使うことにした。(元々は四則演算とn乗と絶対値だけのつもりでした)

そして今思えば、TeXにするのは余分なこだわりだったかなと思っています。

1. 方程式の見える化の環境を作る。



方程式は文字列なので、目に見えるグラフのようなものにしておかないとお絵描きすることはできません。今回は言語としてRを選択したので等高線を書くcontourを使いました。
例えば

$$ 0 = 1 - \left(\frac{x-119}{103} \right)^2 - \left(\frac{y-56}{86} \right)^2 $$

この楕円の方程式であれば、

クリックすると元のサイズで表示します

こんな風に表示できます。

$x$,$y$の2変数の方程式なので、 $z=f(x,y)$ の形の3次元空間内の曲面とみなすこともできて、その方が考えやすいときもあります。
rglを使って3次元の状態を見られるようにもしておきましょう。
ここに貼っているのは画像ですが、Rで実行するとマウスでぐりぐり動かせる3Dモデルになります。

クリックすると元のサイズで表示します

実際には $z=f(x,y)$ ではなく $z=-log(-min(f(x,y),0)+1)$ という変換をしています。$z=0$ と言う平面との交わりを見るためと、$z$方向が大きすぎたり小さすぎたりすると、縦横に比べて高さが大きくなりすぎて、見にくくなってしまうので、それを防ぐための工夫です。
(3次元空間で表現するのは@kohskeさんのアイデアをいただきました。)

function_viewer.r

2. DSLを作る。



いきなりDSL作るのか!と思うかも知れませんが、何度も方程式を書き換えていると面倒な作業が多くなってきて、それを軽減するために作るという感じです。
DSL(ドメイン固有言語:domain-specific language)と言っても、楕円の方定式を作る、ミラー、ユニオン、インターセクション、補集合、などの7つ程度です。

方程式は最終的には文字列にしたいので、上記のDSLも方程式文字列を受け取って新しい方程式文字列を返すという関数の形をしています。

function_dsl.r

いくつか解説しますね。

【楕円を書く】

たとえば

ellipse(119,56,103,86)

と書くと楕円を書きます。端が欠けているのは領域をはみ出しているからです。

クリックすると元のサイズで表示します

【インターセクション】

二項演算子 %I% を使うとインターセクションを書いてくれます。
曲線の内側の領域に注目してインターセクションと呼んでいます。
下の画像は、2つの楕円のインターセクションです。

クリックすると元のサイズで表示します

【ミラー】

ジバニャンはほぼ左右対称なので、ミラー関数mirrorを作りました。

クリックすると元のサイズで表示します

【ユニオン】

二つの領域をくっつけるのがユニオンです。union関数にしました。
最初は %U% という二項演算子だったのですが、2個だけじゃなくて複数のユニオンを取ることが多かったので可変長引数の関数にしてあります。

クリックすると元のサイズで表示します

【カット】

余分な部分を消すのがカットです。二項演算子の %C% です。
インターセクションと似ていますね。実は %C% の実装には %I% を使っています。
%I% を使うよりわかりやすいので作りました。 ( %I% を使った方がわかりやすい時もあるので両方残してあります。)

この画像は、上の例の下半分をカットしています。

クリックすると元のサイズで表示します

こんな風にDSLを活用してジバニャンを作り上げていくのです。

3. ひたすら書く。



まず輪郭をつくっておきます。

クリックすると元のサイズで表示します

つぎに顔の内側を作ります。

クリックすると元のサイズで表示します

おでこのぐにゃぐにゃした模様は、楕円の組み合わせでは難しそうだったので、特別に関数を作りました。
登別の温泉に入りながらひらめいた関数です。

$$ 1 - \left|\frac{x}{51}+\frac{10}{51} sin \left( \frac{7}{2} \pi \left|\frac{y}{1.2*51}\right|^{1.2} \right) \right|^\frac{2}{3} - \left|\frac{y}{1.2*51}\right|^{\frac{2}{3}} $$


クリックすると元のサイズで表示します

これを顔の内側の部分からカットします。

クリックすると元のサイズで表示します

顔の内側ができたので、輪郭と合わせます。

クリックすると元のサイズで表示します

目、鼻、口をつけて出来上がり。

クリックすると元のサイズで表示します

クリックすると元のサイズで表示します

これで完成です。
この後、DSLで構成した方程式が余計な記述があったりするので(2重のカッコとか、2重のマイナスとか、絶対値を2乗しているので絶対値不要とか)、手作業で直しておしまいです。

jibanyan_edit.r

まとめ



こうやってみるとドローツールでジバニャンを描いているだけにも見えますね。
おでこの模様は楕円じゃなくて、実は内心自慢のポイントです。

妖怪ウォッチが大好きな娘に見せたら、口が変、と言われました。
口、最後に作ったし、本気で作ると細かいパーツが多くて大変だったから、さぼったのだけど、一瞬でバレました。さすがのドメイン知識です。

明日のVOYAGE GROUP エンジニアブログ : Advent Calendar 2014は、@co3kです!
0

2014/12/15

MathJaxの練習  



When $a \ne 0$, there are two solutions to \(ax^2 + bx + c = 0\) and they are
$$x = {-b \pm \sqrt{b^2-4ac} \over 2a}.$$
0



teacup.ブログ “AutoPage”
AutoPage最新お知らせ