2015-05-16

ARC037復習 (A,B,C)

開催は4/18だったので一ヶ月近く前のやつですね。

<a href="http://arc037.contest.atcoder.jp/">Welcome to AtCoder Regular Contest 037 - AtCoder Regular Contest 037 | AtCoder</a>arc037.contest.atcoder.jp

このときもなかなかにボロボロでありました。解説聞いたけどちゃんと復習終わってないからやることに。ただD問題は参加時にそもそも問題を読むことすら出来なかったのでまたそのうちに。

A - 全優

未来視した点数が80点を下回るならその分だけ勉強時間を割いてあげましょうって問題。

まあAはいいですよねAは。

N = int(raw_input())
m = map(int,raw_input().split(" "))
print sum([0 if m[i] >= 80 else 80 - m[i] for i in range(N)])

lambda式使えばさらにまとめられるなーと思いましたけど特に意味はなかった。

N = int(raw_input())
print sum(map(lambda x: 0 if int(x) >= 80 else 80 - int(x), raw_input().split(" ")))

B - バウムテスト

与えられた無向グラフにおいて閉路を持たない連結成分(木)の個数を求める問題。

一応通したものの、見返すと40行ぐらいのなかなかひどいコード書いてました。pythonは可読性が売りなのにこれは無い。

方針としては、木を節点の集合(set())で表すことにして、

辺の両端の二節点のどちらも既存の木に属してない → 新しく木を生成
辺の両端の二節点の片方だけが既存の木に属している → 属してない方をその木に属させる
辺の両端の二節点の両方が既存の木に属している→
- 属している木が同じ → 閉路があるので木ではない
- 属している木が異なる → 和集合をとって一つの木にする

それで最後に木の数と一度も使われてない節点の数の合計を返すって方針でした。

N, M = map(int,raw_input().split(" "))
tree_num = 0
tree_set = []
unused = set(range(1,N+1))
for i in range(M):
    u, v = map(int,raw_input().split(" "))
    if u in unused:
        unused.remove(u)
        if v in unused:
            unused.remove(v)
            tree_set.append([1,set([u,v])])
            tree_num += 1
        else:
            for j in range(tree_num):
                if v in tree_set[j][1]:
                    tree_set[j][1].add(u)
                    break
    else:
        if v in unused:
            unused.remove(v)
            for j in range(tree_num):
                if u in tree_set[j][1]:
                    tree_set[j][1].add(v)
                    break
        else:
            for j in range(tree_num):
                if u in tree_set[j][1]:
                    for k in range(tree_num):
                        if v in tree_set[k][1]:
                            if j == k:
                                tree_set[k][0] = 0
                            else:
                                tree_set[j][0] &= tree_set[k][0]
                                tree_set[j][1] |= tree_set[k][1]
                                tree_set.pop(k)
                                tree_num -= 1
                            break
                    break
print sum([tree_set[i][0] for i in range(tree_num)]) + len(unused)

うん、クソですわ。対比の為に置いたけどこれはひどいw 同じ方針でももっと綺麗に書けると思うんですけどめんどくさくなっちゃいました。

で、解説で説明された方針は幅優先探索(dfs)を行うってものでした。探索済みのものにチェックして行って、探索中に探索済みにものに出会ったらそれは閉路だろって方針。実装してみて分かったんですが、単純に全節点を始点としてdfsを行い、一度も探索済みの節点に出会わなければ1を返すってだけで良かったんですね。

def dfs(here, before):
    is_tree = True
    if searched[here]:
        is_tree = False
    else:
        searched[here] = True
        for next in link[here]:
            if next == before:
                continue
            is_tree &= dfs(next, here)
    return is_tree

N, M = map(int,raw_input().split(" "))
tree_num = 0
searched = [0 for i in range(N)]
link = [[] for i in range(N)]
for i in range(M):
    u, v = map(int,raw_input().split(" "))
    u -= 1; v -= 1
    link[u].append(v)
    link[v].append(u)
print sum(dfs(i,-1) for i in range(N))

さっきより大分すっきりしました。シンプルにかける方が筋が良いって言葉が耳に痛い。

C - 億マス計算

N²マス計算を行ったとき小さい方からK番目の数を求める問題。

1≤N≤3.0*10⁴ であるため当然ながら全部を計算することは不可能ですよね。オンタイムのときは行成分と列成分をとりあえずソートすると何か出来そうな気がする、ぐらいしか思いつかなくてもう何も出来ませんでした。

「小さい方からK番目の値がXである」を「X-1以下の数はK個未満だがX以下の数はK個以上ある」って感じで言い換えるの身に付いてないっすねえ。「X以下の数がK個以上あるような最小のX」を見つければ良いっていうのは言われてみればなるほどなんですけど、その場で思いつくのはなかなか出来ない。早く身につけたい発想ですね。

探索は上の方針で二分探索を行えば良いので、あとは「X以下の数がK個以上ある」判定をどう行うかって話だったんですけど、高橋ジグザグ法も聞いてなるほどーてなりました。答え聞いたら当たり前のことではあるんですけど、ああいうのを思いつける人になりたい。

def check(X):
    count = 0 
    j = N - 1
    for i in range(N):
        while j >= 0 and a[i] * b[j] > X:
            j -= 1
        count += j + 1
    return count >= K

N, K = map(int,raw_input().split(" "))
a = map(int,raw_input().split(" "))
b = map(int,raw_input().split(" "))
a.sort()
b.sort()
left = 0
right = 10**18
while right - left > 1:
    mid = (left + right) / 2
    if check(mid):
        right = mid
    else:
        left = mid
print right

この前やったABC002やその前のABC023でも二分探索の復習したんで、次はオンタイムで解けるようにしたいです。今日のARCも頑張ります！

2015-05-15

ABC002

今日のラボプロコンでABC002をやりました。

<a href="http://abc002.contest.atcoder.jp/">Welcome to AtCoder Beginner Contest 002 - AtCoder Beginner Contest 002 | AtCoder</a>abc002.contest.atcoder.jp

最近の奴に比べると難易度が本当に教育的だって思いましたね。この前のABC023の制約とかもうね(;_;) 定期的に記事書かないと習慣がつかないと思ったので上げることに。D問題が結構ためになりました。

A-正直者

入力二つの大きい方を出力する問題。一行で書きたいなら皆こうなるであろう(pythonなら)。

print max(map(int,raw_input().split(" ")))

B-罠

入力された文字列の母音を取り除いて出力する問題。ちょっとfilter使いたかったので使ってみた。

ban_word = set(["a","i","u","e","o"])
print "".join(filter(lambda w: w not in ban_word,list(raw_input())))

C-直訴

与えられた3点がなす三角形の面積を出す問題。ヒント見ちゃったので乗っかっちゃいましたw (0,0), (a,b), (c,d) で構成される三角形の面積は|ad−bc|/2だっての高校でやった覚えあるなあ。

xa, ya, xb, yb, xc, yc = map(int,raw_input().split(" "))
print abs((xa-xc)*(yb-yc)-(xb-xc)*(ya-yc)) / 2.

D-派閥

N人の議員についてのM個の人間関係から、全員が知り合いであるような最大のグループの人数を求める問題。文章にすると分かりにくいですが、いわゆる最大クリーク*1問題です。

最大クリーク問題 - Wikipedia

「あっこれグラフ理論でやった奴だ」ってちょっと思いましたね笑とは言ってもそもそも最大クリーク問題の最適な解き方とかよく知らなかったんで地道に解くことにしました。

初めにこのD問題が一番ためになったと言ったのですが、なぜかというとプロコンでしばしば見られる問題の言い換えを意識することが出来たからです。「グラフの最大クリークの節点数Kを求める」と言われても僕は困ってしまうのですが、「グラフにK個の節点からなるクリークが存在するようなKの最大値を求める」と考えると話が分かりやすくなります。

「グラフにK個の節点からなるクリークが存在するか？」という問題はK-クリーク問題と呼ばれます。これを解くのは(計算量の話は今回は置いとくとして)簡単です。グラフのN個の節点の内K個の節点を選んでしえば, それらがクリークを成しているかは任意の二節点が辺で結ばれているかを調べてしまえば確かめられます。これをN個の節点から取れるK個の節点の組み合わせ, $_NC_K$ 通り全て調べればK-クリーク問題は解けるわけです。

よって今回の問題は, KをNから1まで動かして行って各KでK-クリーク問題を解き, TRUEだったKの最大値はいくらだったかを返せば良いことになります。こういう解き方を初めて自力でできたので簡単な問題とはいえ結構嬉しかったです(^^)

from itertools import combinations as comb

def k_clique_solve(relation, max_members, k):
    for members in comb(max_members,k):
        for r in comb(members,2):
            if r in relation:
                continue
            break
        else:
            return True
    else:
        return False

N, M = map(int,raw_input().split(" "))
max_members = range(N)
relation = set()

for i in range(M):
    x,y = map(int,raw_input().split(" "))
    relation.add((x-1,y-1))
 
result = 1
for k in range(N,0,-1):
    if k_clique_solve(relation,max_members, k):
        result = k
        break
print result

今回は1≤N≤12なので全然大丈夫ですが、 ARCや最近のABCだと制約が厳しいため二分探索を求めてきます。せっかくなのでやってみました。

from itertools import combinations as comb

def k_clique_solve(relation, max_members, k):
    for members in comb(max_members,k):
        for r in comb(members,2):
            if r in relation:
                continue
            break
        else:
            return True
    else:
        return False

N, M = map(int,raw_input().split(" "))
max_members = range(N)
relation = set()

for i in range(M):
    x,y = map(int,raw_input().split(" "))
    relation.add((x-1,y-1))
 
left = 1
right = N + 1
while right - left > 1:
    mid = (left + right) / 2
    if k_clique_solve(relation, max_members, mid):
        left = mid
    else:
        right = mid
print left

まあNがちっちゃいので2ms速くなっただけでしたね笑今回のは良い練習になったと思うので、次は本番でこういう解き方をしたいなあと思います。

*1:任意の二節点が辺で結ばれているグラフをクリークという

2015-05-10

ABC023復習

　昨日AtCoder Beginner Contest 023 に参加しました。結果はボロボロ..orzせめてC問題はちゃんと解きたかったなあ。

<a href="http://abc023.contest.atcoder.jp/">Welcome to AtCoder Beginner Contest 023 - AtCoder Beginner Contest 023 | AtCoder</a>abc023.contest.atcoder.jp

珍しくちゃんと全部復習したので上げることにしました。ブログを書くことで復習へのモチベが上がると嬉しいですね。

A：加算王

　まあワンライナーでいきますよね

print sum(map(int, list(raw_input())))

B：手芸王

　解説では素直にシミュレーションしても十分間に合うってことだったんですが、pythonでstringの連結はあんまりやりたくなかったのでこのような形にしました。真ん中から見て行って、不正なブロックをくっつけている場合はその時点で-1を返します。文字列作ってないだけでやってることは一緒ですね。

N = int(raw_input())
S = list(raw_input())
block = ["a", "b", "c"]
if N % 2 == 0:
    print -1
else:
    center = N / 2
    for i in range(center+1):
        if S[center-i] == block[(1-i) % 3] and S[center+i] == block[(1+i) % 3]:
            continue
        print -1
        break
    else:
        print center

C：収集王

　全マスをチェックするのは計算量的に絶対無理って分かってたんですが、良い手が思いつかず部分点だけ取りに行きました。最後の二重カウントを省くのがN回回すだけで済むとなぜ気付かなかったのか..

from collections import defaultdict as ddict
R, C, K = map(int,raw_input().split(" "))
N = int(raw_input())
rc = [0 for i in range(R)]
cc = [0 for i in range(C)]
candy = [None for i in range(N)]
for i in range(N):
    r, c = map(int,raw_input().split(" "))
    candy[i] = (r-1,c-1)
    rc[r-1] += 1
    cc[c-1] += 1
rn = ddict(int)
cn = ddict(int)
for i in range(R):
    rn[rc[i]] += 1
for i in range(C):
    cn[cc[i]] += 1
count = 0
for key in rn:
    count += rn[key] * cn[K-key]
for r, c in candy:
    if rc[r] + cc[c] == K:
        count -= 1
    elif rc[r] + cc[c] == K + 1:
        count += 1
print count

D：撃墜王

　これもオーダーがきつかったですね。最初から部分点狙いでbit-DPしてみたけどそれでも間に合わないっていう(^^;;) chokudaiさんも仰ってましたけど、二分探索をちゃんと使いこなせるようになりたいです。

def check(X, H, S, N):
    time_limit = [0 for i in range(N)]
    for i in range(N):
        tmp = (X - H[i]) / S[i]
        if tmp < 0:
            return False
        if tmp < N:
            time_limit[tmp] += 1
    count = 0
    for i in range(N):
        count += time_limit[i]
        if count > i + 1:
            return False
    else:
        return True

N = int(raw_input())
H = [0 for i in range(N)]
S = [0 for i in range(N)]
for i in range(N):
    H[i], S[i] = map(int,raw_input().split(" "))
left = max(H) - 1
right = max(H[i] + S[i]*(N-1) for i in range(N))
while right - left > 1:
    mid = (left + right) / 2
    if check(mid, H, S, N):
        right = mid
    else:
        left = mid
print right

いつになったらビギナーを脱せるんですかねえ。まあ頑張ります。

2015-05-09

言語処理の為の機械学習入門をつまみ食い(その1)：パラメータの最尤推定

　最近MLaPP(Machine Learning: a Probabilistic Perspective)の輪講が我がラボで始まったのですが、話が大分抽象的+英語で早くもヤバい空気が流れ始めました(^^;)

そこで、ちょっと機械学習で使う確率系の話を復習する為に、

言語処理のための機械学習入門 (自然言語処理シリーズ)

を読み直そうと思います。この本は二年前から所持しており分からなくなったときにちょくちょく読んでいるのですが、そろそろまとめ的なのを作ろうかと思いここに書くことにしました。

あくまで入門なのですが、簡潔に書いてある＋日本語（超重要！）であるところが素晴らしい。もちろんさらに深層に挑むなら、MLaPPもそうだけどPRML読めって話ですけどね。

そんなこんなでちょびちょびやって行こうと思います。

具体例：コイントスはフェアなのか？

　コイントスに使うコインは表と裏が出る確率が同じだからこそ意味がありますが、ふとイカサマを疑うときがあります。そこで、実際に10回コイントスして調べてみましょう。

コインA : 表, 裏, 表, 裏, 表, 表, 表, 裏, 裏, 裏
コインB : 表, 裏, 裏, 裏, 裏, 表, 裏, 裏, 裏, 裏

すごく大雑把ですが、コインAは表と裏が出る確率はそれぞれ0.5だと考えられ、フェアと言えそうです。逆にコインBは表と裏が出る確率はそれぞれ0.2と0.8、コインへの細工を疑いますね笑

このように、実際に試行したデータを用いることで我々は直感的にコインの表の出る確率（と裏の出る確率）を調べる事が出来ます。今回のテーマである最尤推定法は、すごくシンプルに言えばこのコインの表の出る確率を求める手法です。

コイントスの確率分布の記述(読み飛ばし可)

　コイントスは取り得る値(確率変数X)が表か裏かの二つしかなく、確率θで表が、確率(1-θ)で裏が出ると書けます*1(θは確率分布のパラメータ。 0 < θ < 1)。すなわち確率変数Xが値x(表か裏)をとる確率 $p(X=x;\theta)$ は以下のように書けます。

$p(X=x;\theta) = \delta(x,表)\theta + \delta(x,裏)(1-\theta) \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ = \delta(x,表)\theta + (1-\delta(x,表))(1-\theta) \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ = \theta^{\delta(x,表)}(1-\theta)^{1-\delta(x,表)}$

ここで $\delta(x,a)$ はデルタ関数と呼ばれ、x=aのときに1を、そうでないときに0を返す関数。

※要するに $p(X=表;\theta) = \theta$ 、 $p(X=裏;\theta) = 1-\theta$ になるってことです。これさえ分かってれば問題ないです。

尤度とは

　さて先ほど $p(X=x;\theta)$ を記述しましたが、これにより実データが生成される確率を記述することが出来ます。ここで、コイントスは一回ごとに結果が独立であり、かつ全てが同一の確率分布に従っている(それぞれの試行の確率を $p(X=x;\theta)$ で書ける)*2ので、実データ $D = [x_1, x_2, ... , x_n$ ]の生成確率 $P(D)$ は

$P(D)=\prod_{i}^{n}p(x_i)$

と書けます。それぞれのコイントスの確率の積になってますね。この実データDの生成確率を尤度と呼びます。通常は積の形だと扱いにくいので、対数を取った値、

$\log P(D) = \log \prod_{i}^{n}p(x_i) \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ = \sum_{i}^{n} \log p(x_i)$

を用いることが多く、これを対数尤度と呼びます。

最尤推定

　やっと本題に入れます。最尤推定法とは、尤度が最も大きくなるようにパラメータ(今回はθ)を決定する方法です。対数尤度の方が扱いやすいため、通常は対数尤度を最大化することが多いです。尤度（実際のデータが生成される（起こる）確率）が最も大きくなるようにパラメータを決定することで、パラメータを現実のデータにフィットさせている訳です。

さてコイントスの話に戻ると、コインAのデータ $D_A = [表, 裏, 表, 裏, 表, 表, 表, 裏, 裏, 裏$ ] の対数尤度は

$\log P(D_A) = \sum_{i=1}^{10} \log p(x_i) = 5\log p(表) + 5\log p(裏) = 5(\log\theta + \log(1-\theta))$

これを微分すると

$\displaystyle \frac{d log P(D_A)}{d \theta} = 5( \frac{1}{\theta} - \frac{1}{1-\theta} ) = \frac{5(1-2\theta)}{\theta(1-\theta)}$

となります。直感の通り、このデータにとって最も尤もらしい(※駄洒落ではない)表の出る確率(θ)は0.5となります。

同様にコインBのデータ $D_B =[表, 裏, 裏, 裏, 裏, 表, 裏, 裏, 裏, 裏$ ] についても

$\log P(D_B) = \sum_{i=1}^{10} \log p(x_i) = 2\log p(表) + 8\log p(裏) = 2(\log\theta + 4\log(1-\theta))$

$\displaystyle \frac{d log P(D_B)}{d \theta} = 2( \frac{1}{\theta} - \frac{4}{1-\theta} ) = \frac{2(1-5\theta)}{\theta(1-\theta)}$

となり、これもまた直感の通りθ=0.2となっていることがわかります。非常にシンプルな例ではありますが、最尤推定法でパラメータを求めることが出来ました！

終わりに

　繰り返しになりますが、最尤推定法は実際の観測データにフィットするようにパラメータを決定する方法です。今回の例は非常にシンプルですが、もっと複雑な問題を扱う場合も根本の考えは一緒だと思います。

　ところで、たとえイカサマをしていないコインであっても、10回の試行中偶然にも裏を9回出した場合は最尤推定法でθ=0.1が求まります。実データからはこう求まっても、コインは平等だと思っている我々に取っては少し納得がいきません。θが大体0.5位になりそうだと分かっている場合に、それをθに反映させたいと思うときがあるわけです。こんな場合に用いるのが最大事後確率推定(MAP推定)なのですが、これについてはまたの機会に。

*1:このような確率分布をベルヌーイ分布(Bernoulli distribution)と言う

*2:独立に同一の確率分布に従う(independently, identically distributed)。 i.i.d.としばしば略記される

俵言

しがない社会人が書く、勉強とかのこと。最近は機械学習や kaggle 関連がメイン。