Haskellで数独を解いてみた

関数プログラミング実践入門を読んだのでHaskell数独を解いてみました。実装は以下のようになりました。

module Try.Sudoku where

import Data.List
import Data.Function

type Cell = (Int, Int)
type Board = [(Cell, Int)]


solve :: Board -> [Board]
solve board | length board == 81 = [board]
solve board = do
    concat $ map solve [ (cell, n) : board
                        | let remains = concat $ (\x -> if x `elem` map (\x -> fst x) board then [] else [x]) <$> cells
                        , let cell = maximumBy (compare `on` length . used board) remains
                        , n <- concat $ (\x -> if x `elem` used board cell then [] else [x]) <$> [1..9]]


cells :: [Cell]
cells = [(x, y) | x <- [0..8], y <- [0..8]]

area :: Cell -> Int
area (x, y) = y `div` 3 * 3 + x `div` 3

used :: Board -> Cell -> [Int]
used board cell = nub [ n
                      | (cell', n) <- board
                      , any (\f -> f cell == f cell') [snd, fst, area]
                      ]

sudokuAns :: IO ()
sudokuAns = case solve problem of
    answer : _ -> mapM_ print $ format answer
    []         -> putStrLn "invalid problem"

format :: Board -> [[Int]]
format = map (map snd) . transpose . groupBy ((==) `on`(fst . fst)) . sort

problem :: Board
problem = [ ((3, 0), 8),
            ((5, 0), 1),
            ((6, 1), 4),
            ((7, 1), 3),
            ((0, 2), 5),
            ((4, 3), 6),
            ((6, 3), 8),
            ((6, 4), 1),
            ((1, 5), 2),
            ((4, 5), 3),
            ((0, 6), 6),
            ((7, 6), 7),
            ((8, 6), 5),
            ((2, 7), 3),
            ((3, 7), 4),
            ((3, 8), 2),
            ((6, 8), 6)
            ]  

main :: IO ()
main = do
    sudokuAns

各セルとボードの情報を表すための型シノニムを定義します。

type Cell = (Int, Int)
type Board = [(Cell, Int)]

ボード上には9 × 9のセルが並び、各セルを初期化する場合は以下のようになります。

cells :: [Cell]
cells = [(x, y) | x <- [0..8], y <- [0..8]]

問題となる数独のボードは以下のように表せられます。

problem :: Board
problem = [ ((3, 0), 8),
            ((5, 0), 1),
            ((6, 1), 4),
            ((7, 1), 3),
            ((0, 2), 5),
            ((4, 3), 6),
            ((6, 3), 8),
            ((6, 4), 1),
            ((1, 5), 2),
            ((4, 5), 3),
            ((0, 6), 6),
            ((7, 6), 7),
            ((8, 6), 5),
            ((2, 7), 3),
            ((3, 7), 4),
            ((3, 8), 2),
            ((6, 8), 6)
            ]  

((3, 0), 8) は四行一列目が8で埋められていることを表し、リストを使ってすでに埋められている数値を表しています。

数独では3×3の区間で数値が被らないようにする必要があるので、対象のセルが3×3のどの区間に存在しているか判定するために以下の関数を定義します。

area :: Cell -> Int
area (x, y) = y `div` 3 * 3 + x `div` 3

ボードとセルを指定して、対象のセルに入る値の候補を取得できるようにするため以下の関数ではすでに使われている数値を取得します。

used :: Board -> Cell -> [Int]
used board cell = nub [ n
                      | (cell', n) <- board
                      , any (\f -> f cell == f cell') [snd, fst, area]
                      ]

リスト内包表記を使って入る値の候補を取得しますが重複を除去するためにnubを使っています。 nubを使うと以下のようにリスト内の重複を除去してくれます。

ghci>[n | n <- [0,1,2,2,3,4,5]]
[0,1,2,2,3,4,5]
ghci>nub [n | n <- [0,1,2,2,3,4,5]]
[0,1,2,3,4,5]

(cell', n) <- boardで引数で渡したボードを扱いやすくするようにセルと値に分けています。 any (\f -> f cell == f cell') [snd, fst, area]では指定したセルの縦、横、3×3の区間ですでに使用されているものでフィルターをかけています。anyは関数とリストを受け取り、リスト内の各要素に関数適応後に一つでもTrueを返すものがあれば結果がtrueとなります。

ghci>any (\a -> a == 1) [1, 2, 3]
True

フィルターをかけた後にすでに使われている数値の重複を除去するためnubを使っています。

問題を解いているメインの関数は以下になります。

solve :: Board -> [Board]
solve board | length board == 81 = [board]
solve board = do
    concat $ map solve [ (cell, n) : board
                        | let remains = concat $ (\x -> if x `elem` map (\x -> fst x) board then [] else [x]) <$> cells
                        , let cell = maximumBy (compare `on` length . used board) remains
                        , n <- concat $ (\x -> if x `elem` used board cell then [] else [x]) <$> [1..9]]

数独には複数の回答がありえるため関数の型はsolve :: Board -> [Board]となっております。処理の大まかな流れはボード上で埋められる数の候補が一番多いセルを取得ごにそれぞれのセル埋めた場合で処理を分岐させるのを繰り返し、ボード上の全てのセルが埋まったら答えの一つとして返すというようになっています。 セルが全て埋まっているかどうかは以下で判定しています。

solve board | length board == 81 = [board]

それから以下の関数で再帰的にセルを埋めています。

solve board = do
    concat $ map solve [ (cell, n) : board
                        | let remains = concat $ (\x -> if x `elem` map (\x -> fst x) board then [] else [x]) <$> cells
                        , let cell = maximumBy (compare `on` length . used board) remains
                        , n <- concat $ (\x -> if x `elem` used board cell then [] else [x]) <$> [1..9]]

リスト内包表記では引数で与えられたボードの先頭に一つ埋めたセルを追加するようにしています。[ (cell, n) : board | ,,,]それからmap solve を適用することで再帰的に処理を呼び出すようにしています。最後にconcatを使うことでboard型から[board]型に変換しています。

ghci>concat $ [[0..3],[4..5]]
[0,1,2,3,4,5]

リスト内包表記の中身についてみていきます。

let remains = concat $ (\x -> if x `elem` map (\x -> fst x) board then [] else [x]) <$> cells

についてはremainsにまだ埋められていないcellを代入しています。

let cell = maximumBy (compare `on` length . used board) remains

まだ数値が埋められていないセルの中で入る値の候補が一番多いものを取得しています。

n <- concat $ (\x -> if x `elem` used board cell then [] else [x]) <$> [1..9]]

それから上記の処理により、入る値の候補のみを取得しています。よってこのリスト内包表記ではまだ埋められていないセルのうち入る値の候補が一番多いセルの各値の候補をボードの先頭に追加したボードのリストを返しています。

回答結果の表示は以下の処理になります。

sudokuAns :: IO ()
sudokuAns = case solve problem of
    answer : _ -> mapM_ print $ format answer
    []         -> putStrLn "invalid problem"

format :: Board -> [[Int]]
format = map (map snd) . transpose . groupBy ((==) `on`(fst . fst)) . sort

sudokuAns関数のcase式によりリストのうちの最初の回答のみを評価して表示するようにしています。

実行すると以下のように回答結果が表示されます。

[3,6,9,8,4,1,7,5,2]
[2,7,8,6,9,5,4,3,1]
[5,1,4,3,2,7,9,8,6]
[7,3,1,5,6,4,8,2,9]
[4,9,5,7,8,2,1,6,3]
[8,2,6,1,3,9,5,4,7]
[6,4,2,9,1,8,3,7,5]
[1,5,3,4,7,6,2,9,8]
[9,8,7,2,5,3,6,1,4]

また、この回答方法は複数回答ある場合は先頭一つを取り出すようになっているので高速に結果を返していますが、時間をかけて全ての件数を取得することもできます。

print $ length $ solve problem
19283

Pythonでjanomeを使って形態素解析してみた

pipでjanomeを入れて形態素解析してみた

今までPythonで日本語の形態素解析を入れる時はまずmecabをインストールしてからPythonから呼び出せるように共有ライブラリの依存関係を更新するなどして少し手間がかかっていたのですが、Pythonで作られた形態素解析のライブラリであるjanomeを使えば楽に形態素解析を行えるとのことなので試してみたいと思います。

janomeインストール

pipで一発で入ります。

pip install janome

インストール後の動作確認

使い方についてはgithubのドキュメントを読むので十分そうです。 自分はこんな感じで動作確認してみました。 f:id:steavevaivai:20180616195053p:plain

一般化線形モデルの訓練

一般化線形モデルの訓練についてまとめてみました。
ちゃんと確認が取れているわけではなく、間違っている部分があるかもしれませんがご了承ください。

誤差の削減

損失関数によりモデルの当てはまりの良さの計算が行える。 モデルの訓練では損失の計算、パラメータの更新、損失の計算を繰り返しにより誤差が少なくなるようにします。 モデルの訓練を進める中で損失が変化しなくなったり、パラメータの更新が非常にゆっくりになります、この状態をモデルの収束と呼びます。 モデルの訓練では収束する前までに最良のモデルを見つけることが目的となります。

モデルの訓練の流れ

  1. 教師データから訓練に使うサンプルデータを取り出す
     - フルバッチ: 全ての教師データを使う
     - 確率的勾配降下法(SGD): 1つのサンプルデータを使う
     - ミニバッチ: 複数のサンプルデータを使う
  2. 取り出したサンプルデータを使って損失、勾配を計算する
  3. 勾配を元と学習率を元に現在のパラメータを更新する
  4. 更新後のパラメータで教師データ、テストデータの損失を計算する

サンプルデータの取り出し

パラメータを更新するためには訓練データからデータを取り出して損失を計算する必要があるのですが、データの取り出しかたにはいくつか種類があります。 フルバッチでは全ての訓練データを使います。確率的勾配降下法では訓練データからランダムに一つのデータを取り出します。ミニバッチでは適当な数件のデータを取り出して使用します、件数を増やした方が学習の精度は上がるはずですがその分計算コストがかかるはずなので動かしながらちょうど良いサイズを決めるなどなのでしょうか。

確率的勾配効果法では損失を以下のように求めます。dは訓練データの値とします。
f:id:steavevaivai:20180616095058p:plain
ここでは絶対値を2乗していますが、損失は大きさの単位なので負の数にならないために絶対値を使った方が内容的に正しいのかと思います。

フルバッチでは全ての訓練データを使用するので以下のように表せられます。
f:id:steavevaivai:20180616095129p:plain

ミニバッチでは数件のデータを取り出して一件ごとは確率的勾配降下法と同様に損失を計算するのですが、その後平均をとるようにしています。
f:id:steavevaivai:20180616095142p:plain

勾配の計算

パラメータの更新では損失が小さくなる方向にパラメータを変えていくのですが、そのためにはパラメータと損失の勾配
f:id:steavevaivai:20180616095247p:plain を求める必要があります。

パラメータwと損失Eを以下のグラフで表した時、求めたい勾配はこのようになります。 f:id:steavevaivai:20180616095306p:plain

確率的勾配降下法について、パラメータが1つのシンプルなモデルで勾配を求めてみたいと思います。
f:id:steavevaivai:20180616095318p:plain

w, bの初期値が1でサンプルデータとしてx=3, y=11のデータががあったとします。

求める勾配を以下のようにして、hを小さい数0.001としたら
f:id:steavevaivai:20180616095331p:plain
のようになりE(w+0.001)とE(w-0.001)をそれぞれ求める

E(w+0.001) = 1/2 || y - (w+0.001) * x + b || ^2    
  = 1/2 || 11 - (1+0.001) * 3 + 1|| ^2    
  = 40.4730045    

同様に

E(w-0.001) = 1/2 || y - (w-0.001) * x + b|| ^2    
= 1/2 || 11 - (1-0.001) * 3 + 1|| ^2    
= 40.5270045    

となる。よって

∂E/∂w ≡ (40.4730045 - 40.5270045) / 0.002 = -27    

と勾配は求められる。

勾配の計算を一般化について、例えばパラメータが一つのモデルであれば f:id:steavevaivai:20180616095346p:plain
ここで||d - wx - b ||をuと置いたら E=1/2 ||u||^2となるので
f:id:steavevaivai:20180616095357p:plain
と表すことができる。よって求めたい∂E/∂wは
f:id:steavevaivai:20180616095411p:plainとなる
ここで u = ||d - wx - b ||としているがd - wx - b > 0の時は∂u/∂w = -x
d - wx - b < 0の時は∂u/∂w = xと考えることができるので
f:id:steavevaivai:20180616104610p:plain
と簡単に表すことができる これを使えば先ほどの勾配計算は以下のようになる

∂E/∂w = - (d - y(x;w)x = - (11 - 1*3 -1) * 1
      = - (11 - 1*3 - 1) * 1 = -7

パラメータの更新

勾配を求めたら次のパラメータの更新を行う。現在のパラメータから先ほど求めた勾配に学習率を掛けたものを引いたものが新しいパラメータとなる。
y=wx+bのモデルについて現在のパラメータwが1で勾配が -7, 学習率εが0.01とした場合、アタ新しいパラメータは以下のように求められる。
w_new = w - ε ∂E/∂w = 1 - 0.01 * (-7) = 1.07
この時の学習率が学習の進み方に影響を与え、小さすぎるとゆっくりにしか進まず時間がかかり、大きいと一気に学習が進むが最良のモデルを見つける前に収束するなどします。

学習率の決め方

学習率の決め方については手動で試行錯誤に選んだものを固定で使用するなどありますが、学習の回数に応じて学習率を小さくする(ε = ε0 - αt)方法などもあります。 他の有名な方法としてAdaGradがあり、これは複数のパラメータがあった場合に、パラメータ毎の学習の進み方により学習率を変動させるイメージでパラメータ毎の勾配∂E/∂w_tをg_tとして以下の式で表せられる。
f:id:steavevaivai:20180616095450p:plain
これはパラメータ毎の勾配の各学習毎の合計を使って学習率を変動させている。

モメンタム

モデルの収束性能を向上させる方法としてはモメンタムがある。大きい学習率で学数を進めるとパラメータがジグザグに進むがモメンタムを利用して以下のように更新すると通常の学習に比べて大きくジグザグすることは少なくなる。
f:id:steavevaivai:20180616095502p:plain

Swiftで100マス計算を実装してみた

まず、実装したソースは以下のようになりました。

let col = [1,2,3,4,5,6,7,8,9,10]
let row = [1,2,3,4,5,6,7,8,9,10]

class Answer {
    let col: [Int]
    let row: [Int]
    let ans: [[Int]]
    init(col: [Int], row: [Int], ans: [[Int]]) {
        self.col = col;
        self.row = row;
        self.ans = ans;
    }
}

func hyakumasu(_ a:((_ x: Int,_ y: Int) -> Int), _ col: [Int], _ row: [Int]) -> Answer {
    return Answer(col: col, row: row, ans: col.map { (i) -> [Int] in row.map{ (j) -> Int in a(i, j) }})
}

func showAns(_ ans: Answer) {
    let colAns = ans.ans.map{col -> String in
        col.map{ (a) -> String in ansFormat(a)}.reduce("", + )
    }
    print(ansFormat(0) + ans.col.map{(a) -> String in ansFormat(a)}.reduce("", + ))
    for (a1, a2) in Zip2Sequence(_sequence1: ans.row, _sequence2: colAns) {
        print(ansFormat(a1) + a2)
    }
}

func ansFormat(_ ans: Int) -> String {
    return spaceSeq(4 - ans.description.count) + ans.description
}
func spaceSeq(_ x:Int) -> String {
    return [String](repeating: " ", count: x).reduce("", +)
}

showAns(hyakumasu(+, col, row))
print("")
showAns(hyakumasu(*, col, row))

各処理について

計算結果を保持するためのクラスは以下のように定義しています。

class Answer {
    let col: [Int]
    let row: [Int]
    let ans: [[Int]]
    init(col: [Int], row: [Int], ans: [[Int]]) {
        self.col = col;
        self.row = row;
        self.ans = ans;
    }
}

計算は以下の処理で行なっています。

func hyakumasu(_ a:((_ x: Int,_ y: Int) -> Int), _ col: [Int], _ row: [Int]) -> Answer {
    return Answer(col: col, row: row, ans: col.map { (i) -> [Int] in row.map{ (j) -> Int in a(i, j) }})
}

[Int]型の各行と列の要素同士に対して引数で受け取った(( x: Int, y: Int) -> Int)型の関数を適用するのは以下の処理になります。(( x: Int, y: Int) -> Int)型はInt型の引数を2つ受け取りInt型を返す関数になります。

col.map { (i) -> [Int] in row.map{ (j) -> Int in a(i, j) }})

それから計算結果をAnswer型にして返しています。

計算結果の各項目について、4文字詰でフォーマットして表示するため以下の関数を定義して4文字に足りない分はスペースを詰めて表示が崩れないようにしています。

func ansFormat(_ ans: Int) -> String {
    return spaceSeq(4 - ans.description.count) + ans.description
}
func spaceSeq(_ x:Int) -> String {
    return [String](repeating: " ", count: x).reduce("", +)
}

計算結果の表示には以下の関数を使用しています。

func showAns(_ ans: Answer) {
    let colAns = ans.ans.map{col -> String in
        col.map{ (a) -> String in ansFormat(a)}.reduce("", + )
    }
    print(ansFormat(0) + ans.col.map{(a) -> String in ansFormat(a)}.reduce("", + ))
    for (a1, a2) in Zip2Sequence(_sequence1: ans.row, _sequence2: colAns) {
        print(ansFormat(a1) + a2)
    }
}

colAnsには計算結果のInt型の配列に対してフォーマットをかけて行毎で連結させて[String]にしたものをセットしています。

let colAns = ans.ans.map{col -> String in
    col.map{ (a) -> String in ansFormat(a)}.reduce("", + )
}

print(ansFormat(0) + ans.col.map{(a) -> String in ansFormat(a)}.reduce("", + )) の部分では入力行の配列にフォーマットをかけておいて先頭に0を付与したものを表示しています。

for (a1, a2) in Zip2Sequence(_sequence1: ans.row, _sequence2: colAns) {
    print(ansFormat(a1) + a2)
}

の部分では各列と列毎の計算結果をフォーマットしたものを連結して表示するようにしています。

計算をするときは以下のように適用する関数と行、列の関数を引数として渡しています。

showAns(hyakumasu(+, col, row))
print("")
showAns(hyakumasu(*, col, row))

結果として以下のように計算結果が表示されるのが確認できます。

   0   1   2   3   4   5   6   7   8   9  10
   1   2   3   4   5   6   7   8   9  10  11
   2   3   4   5   6   7   8   9  10  11  12
   3   4   5   6   7   8   9  10  11  12  13
   4   5   6   7   8   9  10  11  12  13  14
   5   6   7   8   9  10  11  12  13  14  15
   6   7   8   9  10  11  12  13  14  15  16
   7   8   9  10  11  12  13  14  15  16  17
   8   9  10  11  12  13  14  15  16  17  18
   9  10  11  12  13  14  15  16  17  18  19
  10  11  12  13  14  15  16  17  18  19  20

   0   1   2   3   4   5   6   7   8   9  10
   1   1   2   3   4   5   6   7   8   9  10
   2   2   4   6   8  10  12  14  16  18  20
   3   3   6   9  12  15  18  21  24  27  30
   4   4   8  12  16  20  24  28  32  36  40
   5   5  10  15  20  25  30  35  40  45  50
   6   6  12  18  24  30  36  42  48  54  60
   7   7  14  21  28  35  42  49  56  63  70
   8   8  16  24  32  40  48  56  64  72  80
   9   9  18  27  36  45  54  63  72  81  90
  10  10  20  30  40  50  60  70  80  90 100

言語処理のための機械学習入門を読んでみた

自然言語処理のための機械学習入門を読んだので載っていた内容を簡単にまとめたいと思います。

1.必要な数学的知識

最適化問題、確率論、情報理論について述べられている

最適化問題

制約付き条件で最適化を行う方法としてラグランジュの未定乗数法についてがある。ラグランジュの未定乗数法は3章クラスタリングEMアルゴリズムと4章分類のサポートベクトルマシンで使用されている。

確率論

確率分布に従い文章内の単語から分類する。パラメータを最適化するための最尤推定、MAP推定などがある。

情報理論

文章内の単語の出現率が独立に同一の確率分布に従う(independently, identically distributed: i.i.d.)という前提の元に最尤推定などのパラメータ最適を行うが、実際のデータは乱雑さがありその尺度であるエントロピーの求める KLダイバージェンスとJSダイバージェンスは文章ないに出現する単語の確率分布の差を求めることが出来る

2.文章および単語の数学的表現

nグラム(n-gram)

nグラムとは隣り合って出現したn単語のことを言う

文書、文のベクトル表現

bag-of-wordsという文書のベクトル表現方法について説明 例えば"nature or nurture? nurture passes nature."という文章は以下のように表現できる xd = (n("nature", d"), n("nurture", d), n("or", d), n("passes",d)) = (2,2,1,1)

文章に対する前処理とデータスパースネス問題

文章内の単語を全てベクトル化するのはあまり有益ではないので、不要な単語は前処理として取り除くようにしておく、取り除く単語のことをストップワードと呼ぶ。 50000語の辞書があった時に100語からなる新聞をベクトル化するとほとんどの要素は0になる、このように0の要素が大きい場合データは疎であるといえデータを処理するために必要な統計量が十分に得られない問題が発生する、このような問題をデータスパースネス問題という

3.クラスタリング

似ているもの同士をグループにするクラスタリングについて述べられている。クラスタリングではラベル付けされた教師データを学習させるのではなく、似ているもの同士でグループを作るようにする。

k-平均法

最初に適当に分けてしまって、それからよりうまく別れるように調整していくことによってクラスタリングを行う方法。

EMアルゴリズム

入力データに含まれない隠れ変数まで考慮してクラスタリングを行う方法

クラスタリングにおける問題点や注意点

クラスタリングではどんなクラスタが出来上がるか事前に知ることができなく、最終的なクラスタ数は何にすれば良いのかという問題がある。 k-平均法やEMアルゴリズムなど繰り返しに基づくクラスタリング手法では結果が初期値に依存するので、いくつかの初期値で計算することが必要になってくる。

4.分類

ナイーブベイズ分類器

ベイズの定理に基づいて各クラスに分類される確率を計算し、最大となったものを分類先として取得する。

SVM(サポートベクターマシン)

SVMは直線を引いて分類を行う、カーネル法と組み合わせることで非線形な分類も可能になる。

本の中では分類手法としてナイーブベイズ分類器とSVMについて説明したが、他にはロジスティック回帰や決定木、それから再帰ニューラルネットワーク と畳み込みニューラルネットワークも使われているようです

5.系列ラベリング

各単語の品詞をタグ付けする方法のことを系列ラベリングと呼んでいる。

隠れマルコフモデル

隠れマルコフモデルは系列ラベリングに用いることが出来るシンプルなモデル。

階層ベイズモデルについての読書メモ

岩波データサイエンス vol.1の読書メモ

階層ベイズ最初の一歩

階層ベイズモデル(階層事前分布を使った統計モデル)からベイズ統計モデルの良さを解説

階層ベイズモデルは似たようなパラメータに制約を与える特徴がある
 → それによりパラメータが多い場合であっても統計モデルをデータにうまく当てはめられる

例題

12才の男子学生(のべ1082人)の身長を測定し、1年後に2回目の測定をしたデータを使う。最初に10調査県で50人ぐらいの男子学生の身長を測定、それからランダムに選ばれた5県に対して普通の給食(甲タイプ)とは異なる乙タイプの給食を食べさせる。1年後に身長を測定し乙タイプの給食が身長の伸びに与える効果を調査する。 正解として乙タイプ給食の効果はゼロで設定してテストデータを生成した。つまり身長増減と給食タイプは無関係ということが導出できれば良い。

まず"直線当てはめ"を使ってみる

まず直線当てはめができる条件として以下を満たしていなければいけない
f:id:steavevaivai:20180506184855p:plain:h18
それぞれ以下を表す Y{i,j}: 測定i回目における県jの平均値 Y{i,j} ~ N(μ{i,j}, σ2): 測定i回目における県jの平均値は平均μ{i,j}で標準偏差がσの正規分布に従う

測定i回目における県jでの身長の平均μ{i,j}は以下になっているものとする。
(測定1の甲・乙) μ
{1,j} = β1 (β1は1回目の測定の平均)
(測定2の甲) μ{2,j} = β1 + β2 (β2は身長の変化量)
(測定2の乙) μ{2,j} = β1 + β2 + β3(β3は乙タイプとの比較量)
以下の説明変数を用意する A_i: 第i回目の測定 X_j: 県毎の異なる給食タイプ(甲タイプが0で乙タイプが1)
f:id:steavevaivai:20180506184910p:plain:h18
この各項についている(β
1, β2, β3)を係数と呼ぶ

μ_{i,j}をデータに当てはめる

動作確認用のデータ(県、標本サイズ、平均身長、身長の標準偏差、測定回、給食タイプが含まれる)をダウンロードしてきてRのglm関数を使ってμ{i,j}のモデルへの当てはまりの良さ(尤度)を最大化するように(β1, β2, β3)を計算してみる。 そしたらβ2 = 5.647, β3 = -1.720(甲タイプの身長増加の平均は5.65cm、乙タイプの身長増加の平均は5.65 - 1.72 = 3.93cm)となった。この時p=0.045という値で有意水準5%に収まっているので乙タイプの給食では身長の伸びが悪くなるという結果が出てしまった。

モデルの不備を見直す

μ_{i,j}へのデータの当てはめを105回繰り返してみる、p < 0.05が成り立つため5000回未満になるはずだがシュミレーションでは10000回以上 p < 0.05となった。このように結果が一致していなのはモデルに不備があるため。

(階層ではない)ベイズモデルにしてみる

以下のシンプルなモデルをベイズモデルにするためにはパラメータ(β1, β2, β_3)の事前分布の指定が必要になる。
f:id:steavevaivai:20180506184855p:plain:h18
f:id:steavevaivai:20180506184910p:plain:h18

事前分布とは
→ パラメータ毎に推定したい範囲を確率分布として指定したもの。どんな値でも良いとしてする場合はむ情報事前分布と呼ばれる。データのばらつきを表すσの事前分布を検討する。

以下の事前分布でベイズモデルに当てはめる - 係数{βk}と標準偏差σある値をとる時、身長平均がY_i,jとなる確率は正規分布N(μ{i,j}, σ2)から求めれる - 係数{βk}は無情報事前分布に従う - 標準偏差σは無情報事前分布に従う これで計算したらβ3は95%の区間(95%ベイズ信頼区間)が-2.411 ~ -0.0065となった

階層ベイズモデルにしてみる

より良い統計モデルを目指してみる。平均ち{Y_i,j}の不確かさ(標準誤差)を(標準偏差) / √{標本サイズ} で近似してみる。測定回i,県j毎の標準誤差をS{i,j}とすることで身長の平均値Y{i,j}を以下のように表せられる。
f:id:steavevaivai:20180506185012p:plain:h18

それから測定i回目における県jでの身長の平均μ{i,j}を決める線形予測子の構成を考える
(測定1の甲・乙) μ
{1,j} = β1 + r{1,j}
(測定2の甲) μ{2,j} = β1 + r{1,j} + β2 + r{2,j}
(測定2の乙) μ
{2,j} = β1 + r{1,j} + β2 + r{2,j} + β3
r
{1,j}とr{2,j}はそれぞれ"1回目の測定における県毎の身長差"、"2回目の測定における県毎の身長増加の差"を表している。これらを合わせると以下の式で表せられる。
f:id:steavevaivai:20180506185022p:plain:h18
事前分布について 係数{b_k}についてはベイズモデルと同様に無情報事前分布を指定 {r_i,j}の事前分布としてはr_i,j ~ N(0, σ2_i)を指定する(県毎の差(r_i,j)は平均0, 標準偏差σ
iの正規分布に従う)

これで計算したら95%ベイズ信頼区間が0に入っている → 甲タイプと乙タイプで差が特にないと言える

階層ベイズモデルの利点

県毎の差を表すパラメータを扱えるようになった。それにより県毎の影響を受けない正確な推定が行えるようになっている。

Hadoopエコシステム周辺について

Hadoopは大規模データの分散処理を支えるオープンソースのソフトウェアフレームワークとして知られていますが、それだけで完結するという訳ではなく他のソフトウェアと連携して利用する必要があったりして初めてHadoopに触れる人などは把握しづらいと思われますので、わかる範囲でHadoopエコシステムとその周辺で使われている技術についてまとめてみたいと思います。

Hadoop

大規模分散処理の基盤として使われているソフトウェアです。Hadoopの機能は分散ファイルシステム(Hadoop Distributed File System)と分散処理基盤(YARN)の2つの機能に分けることができます。

Hadoop Distributed File System(HDFS)

HDFS分散ファイルシステムになりましてブロック単位でファイルを管理しデフォルトでは3ブロックのレプリケーションを作成します。データを保存するサーバではDataNodesのプロセスを動かす必要がありまして、ブロックの保存先を参照するにはNameNodeのプロセスを動かす必要があります。だいたいMasterサーバでNameNodeを動かしてSlaveサーバでDataNodeのプロセスを動かすことが多いかと思います。

Yet Another Resource Negotiator(YARN)

Hadoopの2.2以降で使われているリソース管理、ジョブスケジューリングです。2.2までのHadoopではリソース管理、ジョブスケジューリングとしてMapReduce Version 1が使われていたのですが、ノードとタスクの数が増えるとスケールアウトしないという問題があったのでYARNではそれが修正されています。YARNとMapReduce Version1でどこが違うのかと言いますと、まずMapReduce Version1ではJobTrackerと呼ばれる単一のマスター・プロセスが、クラスターのリソースを管理し、全てのジョブの実行順序の調整まで行います。リソースの管理とジョブの実行順序管理の2つをJobTrackerでまとめて実行する場合スケールアウトできなくなるのでYARNではリソース管理とジョブの実行順序管理を分けて実行するようにしています。YARNではMasterサーバでResourceManagerを動かしSlaveサーバでNodeManagerを動かします、このResourceManagerがHadoopクラスタのリソース管理それからNodeManagerがジョブの実行順序管理を行いましてこのように分けることでMapReduce Version1で発生したスケールアウトしない問題に対策されています。
*ここで出ているジョブとタスクは分散処理の単位になっていましてSparkアプリケーションであればMapメソッドなどの分散処理は複数のタスクで構成され、タスクはジョブで構成されます。ジョブは実際に各Slaveで実行するものになります。sparkアプリケーションは実行時にコアの数、コアごとのスレッドの数を指定するのですがコアの数5のスレッド数1で指定した状態でmap処理を呼び出した場合、データを5分割して処理を実行します。この5分割したものが実際にslaveで実行される単位でジョブになります。

Apache Tez

YARNはリソース管理、ジョブスケジューリングですがApache Tezはその上で動くアプリケーションになっていまして特徴としては処理を有向非循環グラフ(DAG)で繋げることでメモリを有効に活用して処理を行うことができます。詳しくは以下のリンクに詳しく載っています。 https://techblog.yahoo.co.jp/architecture/2015_06_ditributed_system/

Apache Hive

Apache HiveHadoopクラスタ上のRDBです。TEZ上の分散処理エンジンで実行することができます。HiveQLといったクエリを投げることでデータの操作が行えます。HDFS上に保存されたファイルがデータの実体でCSV形式やORC形式で保存されています。大量の時系列データを扱う場合は日付でパーティションを切ってファイルの実体の保存先自体を分けたり、一つのパーティションに大量のファイルがある場合は1つにまとめるなどして読み込みを早くするといった対策はあります。

Apache Spark

Hiveと同じようにDAGを構成する分散処理のエンジンです。バッチ集計やstreaming処理のプログラムを書くことができます。Scalaで作られているのでScala, Javaはもちろん使えてPythonも使えるので扱いやすいのかと思います。

Apache Sqoop

Apache Sqoopは既存のOracleMySQLPostgreSQLなどのRDBからHiveにデータを取り込むためのソフトウェアです。

Apache ZooKeeper

Apache ZooKeeperとは分散処理の環境で同期を行うためのソフトウェアでネーミングベースでロックを取りに行くことで分散処理環境で同期を取れるようにしています。 http://oss.infoscience.co.jp/hadoop/zookeeper/docs/r3.3.1/zookeeperProgrammers.html

presto

Hiveに比べてレスポンスがすぐに帰ってくるようです https://www.slideshare.net/frsyuki/presto-hadoop-conference-japan-2014?ref=https://tug.red/entry/2014/07/10/150250/ https://blog.treasuredata.com/blog/2015/03/20/presto-versus-hive/