Empty

2013年4月から社会人になりました。

気分転換に競技プログラミング: Google Code Jam Qualification Round 2012: Problem C. Recycled Numbers

概要

今日免許の更新に行ってきた。
過去2回の違反(通行禁止違反、進路変更禁止違反)により違反運転者講習に区別され、2時間に渡って大変わかりやすい講習を受けてきた。
講習中余ったリソースを用いて、Google Code Jam Qualification Round 2012 Problem C. Recycled Numbers について考えていたわけだが、その時のアイディアを帰宅後にコード化したものを晒す。まだlargeは解けていないが...

solve.py

#!/usr/bin/env python
# coding: utf-8

"""
Google Code Jam
Qualification Round 2012
Problem C. Recycled Numbers
http://code.google.com/codejam/contest/1460488/dashboard#s=p2
"""

import sys

def read():
    f = sys.stdin
    num = int(f.next().strip())
    return [f.next().strip() for i in range(num)]

def solve(case):
    line = case.split(" ")
    A = int(line.pop(0))
    B = int(line.pop(0))

    def reverse_patterns(num):
        str_num = str(num)
        patterns = {}
        for i in range(1, len(str_num)):
            s1 = str_num[:i]
            s2 = str_num[i:]
            if s1.startswith("0") or s2.startswith("0"):
                continue
            r_num = int(s2 + s1)
            if num == r_num:
                continue
            elif A <= r_num <= B:
                patterns.setdefault(r_num, 0)
        return patterns.keys()

    pairs = {}
    num_set = [n for n in range(A, B+1)]
    for num in num_set:
        for r_num in reverse_patterns(num):
            pair = tuple(sorted([num, r_num]))
            pairs.setdefault(pair, 0)
    return len(pairs.keys())

def main():
    for i, case in enumerate(read()):
        answer = solve(case)
        print "Case #%d: %d" % (i+1, answer)

if __name__ == "__main__":
    main()

気分転換に競技プログラミング: Google Code Jam Qualification Round 2012: Problem B. Dancing With the Googlers 改

前回の続き。
前回はとりあえず全探索的な感じで書いたので、アルゴリズム効率が悪くLargeなケースだと解けずにいたが、先日一緒にこの問題について悩んでいた後輩がその翌日会った矢先に「解けました!」と嬉しそうにPerlのコードを見せてきた。そのコードをPythonで書き直したものを下に晒しておく。(きっと需要はないだろうが)

solve.py

#!/usr/bin/env python
# coding: utf-8

"""
Google Code Jam
Qualification Round 2012
Problem B. Dancing With the Googlers
http://code.google.com/codejam/contest/1460488/dashboard#s=p1
"""

import sys

def read():
    f = sys.stdin
    num = int(f.next().strip())
    return [f.next().strip() for i in range(num)]

def solve(case):
    line = case.split(" ")
    N = int(line.pop(0))
    S = int(line.pop(0))
    P = int(line.pop(0))
    T = map(int, line)
    count = other = 0
    if P == 0:
        return len(T)
    elif P == 1:
        for t in T:
            if t >0: count += 1
        return count
    else:
        for t in T:
            if t >= 3*P-2:
                count += 1
            elif (t == 3*P-3 or t == 3*P-4):
                other += 1
        if other >= S:
            return S + count
        else:
            return other + count

def main():
    for i, case in enumerate(read()):
        answer = solve(case)
        print "Case #%d: %d" % (i+1, answer)

if __name__ == "__main__":
    main()

それにしても、その後輩が直接何かに熱中した姿を見たことがなかったのもあり、競技プログラミングを通してそんな無邪気な一面を見せてくれてこちらも嬉しかった。修了あとわずかだが(無事に修了できる気はしないが)、研究室で定期的に競技プログラミングの問題に取り組めるといいなと思っている今日この頃である。そして今日も修論から現実逃避をしたのであった。

気分転換に競技プログラミング: Google Code Jam Qualification Round 2012: Problem B. Dancing With the Googlers

また研究からの現実逃避をした

ゼミ後に後輩と一緒に Code Jam の問題で楽しんだので気分転換したので、コードを晒してみる。
やった問題は Problem B. Dancing With the Googlers

しかし、アルゴリズムがカスでLargeのデータセットは解けずにいる。
また気が向いた時に改良する予定。

#!/usr/bin/env python
# coding: utf-8

"""
Google Code Jam
Qualification Round 2012
Problem B. Dancing With the Googlers
http://code.google.com/codejam/contest/1460488/dashboard#s=p1
"""

import sys
import itertools

def read():
    f = sys.stdin
    num = int(f.next().strip())
    return [f.next().strip() for i in range(num)]

def solve(case):
    line = case.split(" ")
    N = int(line.pop(0))
    S = int(line.pop(0))
    P = int(line.pop(0))
    T = map(int, line)

    def check_triple(num, S_flg):
        a = int(num / 3)
        r = num % 3 # remainder
        def get_s_triple():
            if r == 0: return (a-1, a, a+1)
            elif r == 1: return (a-1, a+1, a+1)
            elif r == 2: return (a, a, a+2)
        def get_not_s_triple():
            if r == 0: return (a, a, a)
            elif r ==1: return (a, a, a+1)
            elif r ==2: return (a, a+1, a+1)
        if S_flg:
            triple =  get_s_triple()
        else:
            triple = get_not_s_triple()
        for _p in triple:
            if _p < 0: return 0
            if _p >= P: return 1
        return 0

    answer_set = []
    for S_indexes in itertools.combinations(range(N), S):
        pattern = []
        for i in range(N):
            if i in S_indexes:
                S_flg = True
            else:
                S_flg = False
            pattern.append(check_triple(T[i], S_flg))
        answer_set.append(sum(pattern))
    return max(answer_set)

def main():
    dataset = read()
    for i, case in enumerate(dataset):
        answer = solve(case)
        print "Case #%d: %d" % (i+1, answer)

if __name__ == "__main__":
    main()


英語の問題、問題を理解するのに一苦労。
続けていけば、きっと明るい未来がまってるはず。
また近いうちに現実逃避します。

気分転換に競技プログラミング: Google Code Jam Qualification Round 2012: Problem A. Speaking in Tongues

Google Code Jam 2012 の Qualification Round Problem A. Speaking in Tonguesを解いたので、書いたコードを晒してみる。

solve.py

#!/usr/bin/env python
# coding: utf-8

"""
Google Code Jam 2012: Problem A. Speaking in Tongues
http://code.google.com/codejam/contest/1460488/dashboard

$ solve.py < input.in
"""

import sys
import string

sample_set = [
    "ejp mysljylc kd kxveddknmc re jsicpdrysi",
    "rbcpc ypc rtcsra dkh wyfrepkym veddknkmkrkcd",
    "de kr kd eoya kw aej tysr re ujdr lkgc jv",
]

answer_set = [
    "our language is impossible to understand",
    "there are twenty six factorial possibilities",
    "so it is okay if you want to just give up",
]

def understand_sample():
    replace_set = {}
    for i, text in enumerate(sample_set):
        answer_text = list(answer_set[i])
        if not len(text) == len(answer_text): sys.exit("Error sample")
        for j, char in enumerate(list(text)):
            replace_set[char] = answer_text[j]
    if len(replace_set) - 1 == len(string.lowercase):
        pass
    elif len(replace_set) == len(string.lowercase):
        for char in list(string.lowercase):
            if not char in replace_set.keys(): key_char = char
            if not char in replace_set.values(): value_char = char
        replace_set[key_char] = value_char
    elif len(replace_set) + 1 == len(string.lowercase):
        key_chars = []
        value_chars = []
        for char in list(string.lowercase):
            if not char in replace_set.keys():
                key_chars.append(char)
            if not char in replace_set.values():
                value_chars.append(char)
        value_chars.sort()
        key_chars.sort()
        if value_chars == key_chars:
            replace_set[key_chars[0]] = value_chars[1]
            replace_set[key_chars[1]] = value_chars[0]
    return replace_set

def read():
    f = sys.stdin
    num = int(f.next().strip())
    return [f.next().strip() for i in range(num)]

def main(replace_set):
    for i, text in enumerate(read()):
        try:
            replaced_text = "".join([replace_set.get(char)\
                                    for char in list(text)])
        except TypeError:
            raise
        print "Case #%d: %s" % ( i+1, replaced_text )

if __name__ == "__main__":
    replace_set = understand_sample()
    main(replace_set)

A-small-practice.in

30
ejp mysljylc kd kxveddknmc re jsicpdrysi
rbcpc ypc rtcsra dkh wyfrepkym veddknkmkrkcd
de kr kd eoya kw aej tysr re ujdr lkgc jv
hello i am the google code jam test data
how are you
aynny iynny aynny iynny aynny iynny aynny iynny aynny iynny aynny iynny aynny iynny aynny ieeeeeeeee
y n f i c w l b k u o m x s e v z p d r j g a t h a q set k oset xa ynfd
schr rkxc tesr aej dksl tkrb xc
wep rbedc tbe dvcyo ks y resljc ie ser dvcyo re erbcp vcevmc
seneia jsicpdrysid rbcx dksfc rbca ypc dvcyoksl xadrcpkcd ks rbc dvkpkr
rbkd kd de chfkrksl k bygc re le rbc nyrbpeex
kr tyd rbc ncdr ew rkxcd kr tyd rbc nmjpdr ew rkxcd
mcr mkvd ie tbyr bysid ie
rbkd bcpc kd ljsveticp yfrkgyrci rtcsra dcgcs fymkncp wjmm yjre se okfonyfo sykmrbpetksl xyabcx
k bygc ncdrci wpjkr dvkoc ysi xees set k dbymm ncdr aej rbc lja
eb byk kx ks jp fexvjrcp cyrksl aejp fbccqnjplcpd ysi leelmcpcdksl aejp rchrq
ys cac wep ys cac ysi y vklces wep y vklces
ymm aejp nydc ypc ncmesl re cppep rbc dveesa nypi
aej vkddci eww rbc fbkfocs myia
set kd rbc djxxcp ew ejp myfo ew ikdfesrcsr
na rbc vpkfoksl ew xa rbjxnd dexcrbksl tkfoci rbkd tya fexcd
ks y tepmi ew ikpctemgcd ysi mkesd dexcrkxcd rbc pypcdr fpcyrjpc kd y wpkcsi
lpccrksld fbccdc vevdkfmc rbc sjxncp aej bygc ikymci kd fjppcsrma ejr ew vepofbevd
tba ie vpelpyxxcpd ymtyad xkh jv bymmetccs ysi fbpkdrxyd
kx fexxysicp dbcvypi ysi rbkd kd xa wygepkrc vpenmcx es rbc leelmc feic uyx
w ew rte czjymd w ew esc czjymd esc
wep k ncrtccs rbpcc ysi s w ew k czjymd w ew k xksjd esc vmjd w ew k xksjd rte
bet ypc aej bemiksl jv ncfyjdc kx y veryre
ip qykjd ip qykjd ip qykjd ip qykjd eeeeeeeeeeeeb ip qykjd
tbeeeeeeeeeeeeeeeeeeeyyyyyyyyy k oset f vmjd vmjd

output

Case #1: our language is impossible to understand
Case #2: there are twenty six factorial possibilities
Case #3: so it is okay if you want to just give up
Case #4: xoggk d yl wxo vkkvgo ekso uyl wonw sywy
Case #5: xkf yto akj
Case #6: yabba dabba yabba dabba yabba dabba yabba dabba yabba dabba yabba dabba yabba dabba yabba dooooooooo
Case #7: a b c d e f g h i j k l m n o p q r s t u v y w x y z now i know my abcs
Case #8: next time wont you sing with me
Case #9: for those who speak in a tongue do not speak to other people
Case #10: nobody understands them since they are speaking mysteries in the spirit
Case #11: this is so exciting i have to go the bathroom
Case #12: it was the best of times it was the blurst of times
Case #13: let lips do what hands do
Case #14: this here is gunpowder activated twenty seven caliber full auto no kickback nailthrowing mayhem
Case #15: i have bested fruit spike and moon now i shall best you the guy
Case #16: oh hai im in ur computer eating your cheezburgers and googleresing your textz
Case #17: an eye for an eye and a pigeon for a pigeon
Case #18: all your base are belong to error the spoony bard
Case #19: you pissed off the chicken lady
Case #20: now is the summer of our lack of discontent
Case #21: by the pricking of my thumbs something wicked this way comes
Case #22: in a world of direwolves and lions sometimes the rarest creature is a friend
Case #23: greetings cheese popsicle the number you have dialed is currently out of porkchops
Case #24: why do programmers always mix up halloween and christmas
Case #25: im commander shepard and this is my favorite problem on the google code jam
Case #26: f of two equals f of one equals one
Case #27: for i between three and n f of i equals f of i minus one plus f of i minus two
Case #28: how are you holding up because im a potato
Case #29: dr zaius dr zaius dr zaius dr zaius ooooooooooooh dr zaius
Case #30: whoooooooooooooooooooaaaaaaaaa i know c plus plus

ながーいコードになってしまい悔しいです。久々に競技プログラミングすると楽しいですね。また気分転換にやりたいと思います。

最近

概要

最近のことについてだらだらと綴るエントリー.思考の整理になればいいなと.

研究

熱心に指導して頂いている先生のおかげで、残り一年を全力を注ぐべき研究テーマが決まり、やる気満々な状況。
だが、最近論文なりゼミ資料を作っていて、自分の情報をまとめる力・文章力の欠如を感じている。情報を調べて良い感じのタイミングでまとめるというのが苦手みたい。
その点を意識して地道に克服していきたい。今年は国際学会2本出す! 夏と来春。(英語苦手だけど・・・)

勉強会

3ヶ月全10回、運営してきた勉強会が先週の木曜に本の内容が終わると共に終焉となった。
主催者として至らない点が多々あり、終始、参加者の皆さんに助けてもらいながらの運営だった。個人的に実りの多い活動だったと思っている。

参加者の皆さん、一緒に運営してくれた @haman29 に感謝しています。ありがとうございました m(_ _)m

収穫としては、技術的な知識はもちろんだが、それよりも、魅力的な参加者の皆さんとの交流できたこと。一番は、自分が以前より人を巻き込んで物事を進めていくことが出来るように変われたことが最大の収穫じゃないかなと思っている。
その他の感想としては、いつも成果発表時間に成果物示せられなかったのが悔しいので、決められた時間内に最大のアウトプットを出すことを意識していろいろと精進したい。

参加者の皆さんから勉強会への感想として頂いた意見の中で、

  • 毎週必ず行われるところが、勉強する習慣を作ってくれてよかった。

という意見を複数人から頂き、参加者の一人として、なかなか予習する時間を確保するのに大変さを感じ、週一の開催はキツいな~と思っていたので、意外な発想だった。
何かをやるための仕組みを作るって以前よんだ本: 最少の時間と労力で最大の成果を出す「仕組み」仕事術 にも書いてあったような・・・。時間の使い方が下手な自分には無理くりに、勉強会という、技術的にインプット・アウトプットする時間を定期的に作るのは向いている仕組みかもしれない。

今後の方針としては、インプット・アウトプットがバランスが丁度いい感じの勉強会スタイルを模索している。テーマは以前、集合知関連をに興味があるのでそっち系(Data Mining, Machine Learning)の何かをやる予定。
今回の入門ソーシャルデータは、各種Social APIをたたいて、多様な視覚化ツールなどを試したりとか、実践よりな印象を持っているが、個人的な嗜好では、応用の部分に重点を置くよりも、基礎理論から応用部分まで体系的な理解をしたいなと思っているので、
機械学習のアルゴリズムの理解して、○○な課題に自分でアルゴリズムを実装して試してみるとか系。まさに集合知プログラミングはそんな感じだった。入門ソーシャルデータも手軽に横断的に効率的に知識収集する用途には最適な本で勉強になったんだけど、個人的には基礎理論を抑えた上で応用してみるといった教科書的な流れでやっていけたらいいなーと思ってる。

また違った方向性としては、コンテスト型の勉強会。一つの課題を出して、それを参加者各自で最適だと思うアプローチを考えて、実装し適用。で発表してもらって参加者間のアイディアをシェアする形式。
(例えば、A,Bに分類したいデータがあって、どうしたら精度よく分類できるか考えて実装・・・)

まだまだほとんど白紙に近いので、何かいいアイディアがあれば、お気軽に教えてください。
4月下旬もしくはGW明けにできればいいなと考えている。


私生活

最近は、ほぼ毎日研究室で起きて、研究、技術的な勉強、サーバ管理者の仕事、たまに就活、といったことをして、眠くなったら研究室のソファーで寝るという、ある意味とても充実した日々を送っている。
現在の借りている家賃が、比較的高いのもあって、もったいないので、激安な物件に引越をすることにした。来月の中旬~下旬に引越する予定。荷造りや各種手続きが億劫。


就活

5社受けていて、お祈り経験はまだなく順調なのかなと。しかし、まだまだ自己分析が甘い部分があるので、4月以降の選考開始までに煮詰めないといけない。
自分の中で、どんな人になりたいのか・・・とてもフワフワしている状況。携わった会社の社員さん達にとてもお世話になっているので誠意をもって就職活動をしていく所存。


サーバ管理者

ちょうど1日前、ファイルサーバのHDD障害でRaid5のシステムが壊れるという大きめの障害が発生した。(HDD2本がイっちゃったみたい。)
急遽、助教の先生と一緒に新たなサーバを立てることになり、おかげさまで、この2日間はFreeBSDさん、バックアップ用のHDDさん達、助教の先生と長い時間を過ごした。

でも、奇跡的に4~5日前に運用開始後7年目にして初のフルバックアップデータを取っていたので助かった。(今までは共有ディレクトリだけの日時バックアップという運用体制。)

作業内容としては、新たなファイルサーバ/NFSサーバの構築と最新のバックアップデータ/ユーザアカウント情報の移行というシンプルな内容。
でもスムーズにいかなかった。Samba2系->3系への移行だったり、FreeBSD 6系 -> 9系への以降だったりと、そのままの移行だけでは動かないこと多く時間がかかった。
。30時間くらいかけてどうにか復旧。

運用しているシステムの障害から復旧という経験が初めてだったので、いろいろと学ぶことが多くいい経験になった。ちゃんと復旧できたのも自信になった。

しかし、中間発表のResume提出を間近に控える時期に起きたのは、正直勘弁してほしかったなーと思う点もあったが、実践の中で知識を得て・すぐに試して・結果復旧となり研究室に貢献できて、とてもやりがいを感じられてよかった。
障害対応がインフラエンジニアを強くする~という話を身をもって体験できた。

今のアーキテクチャや運用体制はイケてないので、在籍している間に耐障害性を考慮した単一障害点のない・管理しやすい、かっこいいアーキテクチャーにしたいな~と妄想している。まずはFreeBSDさんからDebianさんに乗り換えたい。


総括

充実しているが、プライベート皆無。
それが空しく感じることがないといえば嘘になるので、自分の幸せな状態ってどういう状態なのかを考えて、答えを出して本当の意味のリア充になりたいものです。

ホワイトボードプログラミング

概要

さきほど、就活の選考の一環で、ホワイトボードプログラミングしてきたので、感想とか。

問題

2つのストリング(検索対象文字列/検索文字列)があり、検索文字列が含まれているかどうかを判定する関数を書いてください。

書いた

Cか、C++がいいんだけど・・・pythonでもいいからということで、時間いっぱいに使って書いたのがコレ。

def search (text, keyword):
    k_len = len(keyword)
    for i in range(len(text)):
        _i = i
        for j, _k in enumerate(keyword):
            if text[_i] == _k:
                if (j+1) == k_len: return True
                _i += 1; continue
            else: break
    return False

研究室に帰って確認したら、一応探索出来ているけど、両者の文字列が似通っていて長かったりするとすごい時間がかかる気がする。

感想

こんな基本的なことがすぐにコード化できなくて悔しかった。いい経験させてもらいました。
面接官の方、ホワイトボードに向かってぶつぶつ言っている気持ち悪い男と、20分くらい同室でお付き合いしていただきありがとうございました。

いい刺激と、来週頭に届くお祈り、ありがたく頂戴致します。
アルゴリズムの勉強再開します。

マンガでわかる統計学を読んだ

概要

最近は就活の時期ということもあり、人と会うことが多くなり、どんなことをしたいのか? と聞かれるたびに、「データマイニング機械学習の技術を使って世の中をより良くしたい。」といったことを良く言及している。
しかし、いざ専門的な知識を持っているのかと聞かれると自信を持ってあるとは言えない自分がいる、かつ基礎的な知識についても怪しい節を感じている始末。
で、これではイカンっ!と思って、前々から読みたかったけど積んでいるままであった 「マンガでわかる統計学」シリーズを手に取って再勉強したので、そのメモなどをエントリーする。

読んだ本

今回は上の本のまとめ。回帰分析編はまた後にまとめる予定。

目次

第1章 データの種類をたしかめよう!
第2章 データ全体の雰囲気をつかもう! (数量データ編)
第3章 データ全体の雰囲気をつかもう! (カテゴリーデータ編)
第4章 基準値と偏差値
第5章 確立を求めよう!
第6章 2変数の関連を調べよう!
第7章 独立性の検定をマスターしよう!

カテゴリーデータと数量データ (第1章あたり)

カテゴリーデータ
測れないデータ
選択肢の幅が等間隔でないもの。

e.g.

  • Q. ご感想
    • ["おもしろかった", "ややおもしろかった", "ふつう", "ややつまらなかった", "とてもつまらなかった"]
  • Q. 英検の級
    • ["1級", "準1級", "2級", ... ]
      • 級間の対象語数や難易度の差は等しくないのでカテゴリーデータとして扱うのが一般的。
数量データ
測れるデータ
ある目盛りと隣の目盛りが等間隔であるデータ。標準化されているもの?(上手く定義できない...)

e.g.

  • 身長、体重、実行時間、とか
カテゴリーデータを数値データとして扱う (第1章あたり)

扱いやすいという理由から、カテゴリーデータを数量データとしてみなす場合も多い。
e.g.

とてもおもしろかった 5 points
ややおもしろかった 4 points
どちらともいえない 3 points
ややつまらなかった 2 points
とてもつまらなかった 1 points

度数分布表 & ヒストグラム (第2章~4)

用語
階級
ある範囲で区切ったときのそれぞれのカテゴリー e.g. 300円以上400円未満の商品
階級値
ある階級の値域の中央値 e.g. 350(円)
度数
ある階級に所属している数 e.g. 15(個)
相対度数
全体を1としたときの割合 e.g. 全体:100(個) ある階級の度数:15(個) その階級の相対度数: 0.15
ヒストグラム
底辺を階級の幅に、面積が度数を表すように高さを定めた長方形を密着して並べて作った統計グラフの一種。なんとなーくデータ全体の雰囲気を直観的に掴むことができる.
平均(arithmetric mean)
相加平均・算術平均を指すことが多い。全データの合計/データ数 で求まる。 (その他に、幾何平均、調和平均といった違う概念を指す場合もある。)
中央値(median)
有限個のデータを小さい順位並べたとき、中央に位置する値。データ数が偶数個の場合は、真ん中2つの算術平均をとる。
標準偏差(standard deviation)
散らばりの程度を表す。(1データあたりの平均からのズレ) 0以上。母集団の標準偏差と、標本の標準偏差を求める式は若干異なるが、ほとんどの場合は母集団全体を調べられないので、標本の標準偏差を用いる場合が多い。
階級の幅を決め方
  • スタージェスの公式を用いて幅を決めるやり方が紹介されていた。
  • そのほかの方法、[入門解析法 永田先生 p.9]にも載ってた.
コードに起こす (理解度チェック)

本の内容ではないけど、理解しているならコードに起こせるよね?という考えで組み込み関数を極力使わずに、
本の内容を計算するためのスクリプトを書いた。

#!/usr/bin/env python
# coding: utf-8

import math
import random
from prettytable import PrettyTable


def struges_fomula(dataset):
        """スタージェスの公式"""
        category_num = 1 + math.log(len(dataset), 10) / math.log(2, 10)
        max_v, min_v = max(dataset), min(dataset)
        return float(max_v - min_v) / category_num


def original_fomula(dataset):
        """教科書に載ってた方法
        参考: 入門統計解析法 日科技連出版社 p.9
        """
        m = 1 # 測定範囲
        n = len(dataset) # データ数
        max_v, min_v = max(dataset), min(dataset)
        R = max_v - min_v # 範囲
        h = math.sqrt(n) # 仮の区間数
        return int( (R/h) / m) * m # 区間の幅 測定範囲の整数倍に丸め込む


class statistics():

        def __init__(self, n=100, max=100):
                self.make_dataset(n, max)

        def make_dataset(self, n, max):
                """無作為にデータセットを作る"""
                self.dataset =  [random.random() * 100 for i in range(n)]

        def set_dataset(self, dataset):
                """指定のデータセットを用いる"""
                self.dataset = dataset

        def calculate(self):
                """統計量を算出する"""
                self._calculate_average()
                self._calculate_standard_deviation()
                self._calculate_T_score()

        def _calculate_average(self):
                """平均を算出する"""
                sum = 0.0
                for data in self.dataset:
                        sum += data
                self.average = sum / len(self.dataset)

        def _calculate_standard_deviation(self):
                """標準偏差を算出する"""
                sum = 0.0
                for data in self.dataset:
                        sum += pow( data - self.average, 2)
                self.standard_deviation = math.sqrt(sum / len(self.dataset))

        def _calculate_T_score(self):
                """偏差値を算出する"""
                self.T_scores = []
                for data in self.dataset:
                        standard_value = float(data - self.average) / self.standard_deviation
                        self.T_scores.append(standard_value * 10 + 50.0)

        def dump_data(self):
                """算出した統計量を表示する"""
                # display each score and each T-score
                fields = ["id", "score", "T-score"]
                pt = PrettyTable(fields=fields)
                for i, data in enumerate(self.dataset):
                        pt.add_row([i+1, data, "%2.2f" % self.T_scores[i]])
                pt.printt()

                # display average and standard deviation
                static_fields = ["average", "standard_deviation"]
                static_pt = PrettyTable(fields=static_fields)
                static_pt.add_row(["%2.2f" % self.average, "%2.2f" % self.standard_deviation])
                static_pt.printt()


        def histogram(self, func=struges_fomula):
                """ヒストグラムを作成する"""
                category_distance = func(self.dataset)
                sorted_dataset = sorted(self.dataset)
                category_start_v = min(self.dataset)
                category_end_v = category_start_v + category_distance
                category_freq = []
                _count = 0
                for val in sorted_dataset:
                        while True:
                                _key = "%s ~ %s" % tuple(map(int,[category_start_v, category_end_v]))
                                if category_start_v <= val < category_end_v:
                                        _count += 1
                                        break
                                else:
                                        category_start_v = category_end_v
                                        category_end_v += category_distance
                                        category_freq.append({"key":_key, "count":_count})
                                        _count = 0
                # display histogram
                fields = ["range", "visualization", "freqquency", "conposition ratio (%)"]

                pt = PrettyTable(fields=fields)
                pt.set_field_align("visualization" , "l")
                for info in category_freq:
                        pt.add_row([info["key"], "#" * info["count"], info["count"], "%3.2f" % (float(info["count"]) * 100 / len(self.dataset))])
                pt.add_row(["sum", "", len(self.dataset), "%3.2f" % 100])
                pt.printt()


if __name__ == "__main__":
        # 日本史 得点
        history_dataset = [73, 61, 14, 41, 49, 87, 69, 65, 36, 7, 53, 100, 57, 45, 56, 34, 37, 70]
        # 生物 得点
        biology_dataset = [59, 73, 47, 38, 63, 56, 15, 53, 80, 50, 41, 62, 44, 26, 91, 35, 53, 68]
        # p.78
        practice_dataset = [16.3, 22.4, 18.5, 18.7, 20.1]


        sta = statistics()

        # 基本統計量 算出 about 日本史
        print("----- 日本史 -----")
        sta.set_dataset(history_dataset)
        sta.calculate()
        sta.dump_data()
        sta.histogram()

        # 基本統計量 算出 about 生物
        print("----- 生物 -----")
        sta.set_dataset(biology_dataset)
        sta.calculate()
        sta.dump_data()
        sta.histogram()
        #sta.histogram(original_fomula)

        sta.set_dataset(practice_dataset)
        sta.calculate()
結果(一例)
----- 日本史 -----
+----+-------+---------+
| id | score | T-score |
+----+-------+---------+
| 1  |   73  |  58.79  |
| 2  |   61  |  53.52  |
| 3  |   14  |  32.85  |
| 4  |   41  |  44.72  |
| 5  |   49  |  48.24  |
| 6  |   87  |  64.95  |
| 7  |   69  |  57.04  |
| 8  |   65  |  55.28  |
| 9  |   36  |  42.53  |
| 10 |   7   |  29.77  |
| 11 |   53  |  50.00  |
| 12 |  100  |  70.67  |
| 13 |   57  |  51.76  |
| 14 |   45  |  46.48  |
| 15 |   56  |  51.32  |
| 16 |   34  |  41.65  |
| 17 |   37  |  42.96  |
| 18 |   70  |  57.47  |
+----+-------+---------+
+---------+--------------------+
| average | standard_deviation |
+---------+--------------------+
|  53.00  |       22.74        |
+---------+--------------------+
+---------+---------------+------------+-----------------------+
|  range  | visualization | freqquency | conposition ratio (%) |
+---------+---------------+------------+-----------------------+
|  7 ~ 24 | ##            |     2      |         11.11         |
| 24 ~ 42 | ####          |     4      |         22.22         |
| 42 ~ 60 | #####         |     5      |         27.78         |
| 60 ~ 78 | #####         |     5      |         27.78         |
| 78 ~ 96 | #             |     1      |          5.56         |
|   sum   |               |     18     |         100.00        |
+---------+---------------+------------+-----------------------+
...

5章 確立密度関数

ヒストグラムにおける階級の幅を極限まで狭めた曲式が確率密度関数であるという話で、
主要な確立分布を紹介していた。

  • 正規分布
  • 標準正規分布
  • カイ2乗分布
  • t分布
  • F分布

分布表の見方などを説明していた。x軸と確率密度関数の囲む面積は1である。とかの話。

2変数の関連を調べよう! (第6章)

数量データ と 数量データ
単相関係数
数量データ と カテゴリーデータ
相関比
カテゴリーデータ と カテゴリーデータ
クラメールの連関係数

負の相関, 無相関, 正の相関
○○以上なら相関があるといったことが言える統計学的な基準はない.

これらの指標はあくまでも、直線的な相関が見られる場合のみの判断ということに注意。

コードに起こす (理解度チェック)

また、少しコードを書いて理解を深めた。

  • 数量データ と 数量データ: 単相関係数
# numerable_data * numerable_data
dataset = [
        (3000, 7000),
        (5000, 8000),
        (12000, 25000),
        (2000, 5000),
        (7000, 12000),
        (15000, 30000),
        (5000, 10000),
        (6000,15000),
        (8000, 20000),
        (10000, 18000),]

def cal_cc(dataset):
        """calculate simple correlation coefficient"""
        N = len(dataset)
        sum_x, sum_y = 0.0, 0.0
        Sxx, Syy, Sxy = 0.0, 0.0, 0.0
        for x, y in dataset:
                sum_x += x
                sum_y += y
        ave_x = sum_x / N
        ave_y = sum_y / N

        for x, y in dataset:
                Sxx += pow(x - ave_x, 2)
                Syy += pow(y - ave_y, 2)
                Sxy += (x - ave_x) * (y - ave_y)

        return Sxy / math.sqrt(Sxx * Syy)
  • 数量データ と カテゴリーデータ: 相関比
# numerable_data * category_data
dataset2 = [
        (27, "Hermes"),
        (33, "CHANEL"),
        (16, "BURBERRY"),
        (29, "BURBERRY"),
        (32, "CHANEL"),
        (23, "Hermes"),
        (25, "CHANEL"),
        (28, "Hermes"),
        (22, "BURBERRY"),
        (18, "BURBERRY"),
        (26, "CHANEL"),
        (26, "Hermes"),
        (15, "BURBERRY"),
        (29, "CHANEL"),
        (26, "BURBERRY"),]

def cal_cr(dataset):
        """calculate correlation ratio"""
        N = len(dataset)

        label_count = {}
        sum_all = 0.0
        sum_each_label = {}
        sum_squre = {}
        # 級間変動 aboung class variation
        among_class_variation = 0.0
        # 級内変動 within class variation
        within_class_variation = 0.0

        for age, label in dataset:
                sum_each_label.setdefault(label, 0.0)
                label_count.setdefault(label, 0)
                sum_each_label[label] += age
                sum_all += age
                label_count[label] += 1

        for age, label in dataset:
                sum_squre.setdefault(label, 0.0)
                sum_squre[label] += pow(age - sum_each_label[label]/label_count[label], 2)

        within_class_variation = sum(sum_squre.values())

        for label in sum_each_label.keys():
                among_class_variation += label_count[label] * \
                                pow( (sum_each_label[label] / label_count[label]) - (sum_all / N), 2)
        return among_class_variation / (within_class_variation + among_class_variation)
  • カテゴリーデータ と カテゴリーデータ: クラメールの連関係数
# category_data * category_data
dataset3 = {
        "WOMAN":{"PHONE":34, "MAIL":61, "FACE":53},
        "MAN": {"PHONE":38, "MAIL":40, "FACE":74},}

# category_data * category_data
practice_set = {
        "Japanese":{"COFFEE":43, "TEA":33},
        "Westerm":{"COFFEE":51, "TEA":53},
        "Chinese":{"COFFEE":29, "TEA":41}}

def cal_cramer(dataset):
        """calculate Cramer's coefficient of association"""
        sum_each_sex = {}
        sum_each_method = {}
        row_num = len(dataset.keys())
        column_num = len(dataset[dataset.keys()[0]].keys())

        for sex in dataset.keys():
                sum_each_sex[sex] = sum(dataset[sex].values())
                for method in dataset[sex].keys():
                        sum_each_method.setdefault(method, 0.0)
                        sum_each_method[method] += dataset[sex][method]

        N = sum(sum_each_sex.values())

        # カイ二乗統計量
        chi_squared_statistic = 0.0

        for sex in dataset.keys():
                for method in dataset[sex].keys():
                        # 期待度数の計算
                        expected_frequency = float(sum_each_sex[sex] * sum_each_method[method]) / N
                        chi_squared_statistic += pow(dataset[sex][method] - expected_frequency, 2) / expected_frequency

        return math.sqrt(chi_squared_statistic / \
                                (N * ( min([row_num, column_num]) -1 )) )

独立性の検定をマスターしよう! (第7章)

「検定」とは母集団について分析者が立てた仮説が正しいかどうかを標本のデータから推測する分析手法のこと
正確には「統計的仮説検定」という。

独立性の検定, 相関比の検定, 無相関の検定, 母平均の差の検定, 母比率の差の検定など検定の種類はいろいろとあるが、分析手順は基本的に同じ。

  1. 母集団の定義する.
  2. 帰無仮説と対立仮説を立てる.
  3. どの「検定」を行うかを選択する.
  4. 有意水準を決定する.
  5. 標本のデータから検定統計量の値を求める.
  6. 求めた検定統計量の値が棄却域に入っているかどうかを調べる.
  7. 棄却域に入っていたならば、「対立仮説は正しいとの結論」そうでなければ、「帰無仮説は誤っているとはいえない」との結論を下す.
疑問

理解はできたつもりだが、実務で「検定」をどのようにして利用すればいいのかが思いつかなかった。
無相関であることを示したり、この因子は目的変数に効果を与えていないことを示すなど例だけでは、現実世界のどこで使われているのかがイメージできなかった。
また、教科書などで検定をする際に統計量を求めるときには、予め、「自由度2のカイ2乗分布に従う」、「平均が○○で標準偏差が××の正規分布とみなして」...といった前提が与えられているから、問題は解けるが、実際のデータでは、どの分布に従っているかなんてわからんし、そこをどうすればいいのかな~なんてモヤモヤしていた。

まだまだ勉強なり続けます。近いうちに回帰分析編の内容もまとめられたらいいな。