Haskellでドメイン固有言語

• こちら(GADTとDSL)とこちら(確率的プログラミングとGADT/DSL)とHaskellのEDSL
• GADTs(Generalized Algebraic Data Types)は、Haskellのデータ型である代数型データタイプを一般化したもの
• これを使って、AST(abstract syntax tree)を作ることができる。このASTは構文木から意味に関係ある部分だけを取り出した木構造である。これを作ることで、この代数型データタイプの処理に特化した構文の構成が決まり、それは、(あたかも)(Haskellの中に)、一つの言語(特定の用途に特化した言語〜ドメインと固有言語)が構築されたことと同じである。この埋め込まれた言語を Embedded Domain Specific Language(EDSL)と言う
• 確率的プログラミングをHaskellで実装していくと、Haskell上にProbabilisitic programmingに特化したEDSLを作ることになる、とそういう話
• DSLは抽象的なものなので、意味を持たせるためにはインタープリタの実装も必要なのだそうだ…
• さてやってみよう
  • こちら(確率的プログラミングとGADT/DSL)をなぞってみる

stack new GADTtest

• として、環境を作り、cabal, yamlを修正する。言語拡張とかパッケージとかの要求、それのバージョン指定とか、stack buildでメッセージが出るたびに対応していくだけで、最終的にはコンパイルできるようになる
  • cabal

executable GADT-exe
  hs-source-dirs:      app
  main-is:             Main.hs
  ghc-options:         -threaded -rtsopts -with-rtsopts=-N
  build-depends:       base
                     , GADT
                     , rvar
                     , random-fu
                     , plots
                     , diagrams-lib
                     , diagrams-rasterific
  default-language:    Haskell2010

• • yaml

extra-deps:
- erf-native-1.0.0.1
- polynomial-0.7.3
- plots-0.1.0.2
- diagrams-lib-1.4.1.2

• 次のファイルをMain.hsとして、stackで作った領域のapp/　以下に置き、stack buildする

{-# LANGUAGE GADTs #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE FlexibleContexts #-}

import Data.Random
import qualified Data.Random as Random
import Plots
import Diagrams.Prelude hiding (scale, (<>), Vector)
import Diagrams.Backend.Rasterific.CmdLine
import qualified Data.Random.Distribution.Bernoulli as Bernoulli
import Data.List

type Prob = Double

data Dist a where
  Pure :: a -> Dist a
  Bind :: Dist b -> (b -> Dist a) -> Dist a
  Primitive :: RVar a -> Dist a
  Conditional :: (a -> Prob) -> Dist a -> Dist a

condition :: (a -> Prob) -> Dist a -> Dist a
condition = Conditional

instance Functor Dist where
  fmap f d = d >>= (pure . f)

instance Applicative Dist where
  pure = Pure
  f <*> a = a >>= \x -> fmap ($x) f

instance Monad Dist where
  (>>=) = Bind

instance Sampleable Dist IO a where
  sampleFrom _ (Pure x) = pure x
  sampleFrom s (Bind d f) = do
    x <- sampleFrom s d
    sampleFrom s (f x)
  sampleFrom s (Primitive d) = sampleFrom s d
  sampleFrom _ (Conditional _ _) = undefined


bernoulli a = Primitive $ Bernoulli.bernoulli a

prior :: Dist a -> Dist (a, Prob)
prior (Conditional c d) = do
  (x,s) <- prior d
  return (x, s * c x)
prior (Bind d f) = do
  (x,s) <- prior d
  y <- f x
  return (y,s)
prior d = do
  x <- d
  return (x,1)

mh :: Int      -- burn-in 期間
   -> Int      -- サンプリング数
   -> Dist a   -- `Conditional`を含む事前分布
   -> Dist [a] -- 事後分布からサンプリングした点列上の確率分布
mh dc tc d = fmap (map fst . drop dc) $ proposal >>= iterate (dc + tc)
  where
    proposal = prior d
    iterate 0 _ = pure []
    iterate n (x,s) = do
      (y,r) <- proposal
      accept <- bernoulli $ min 1 (r / s)
      let next = if accept then (y,r) else (x,s)
      rest <- iterate (n-1) next
      return $ next:rest

normalry a b = Primitive $ Random.normal a b

type Point = (Double, Double)
type Param = (Double, Double)

linear :: Dist Param
linear = do
  a <- normalry 0 5
  b <- normalry 0 5
  return (a,b)

point :: Main.Point -> Dist Param -> Dist Param
point (x,y) = condition (\(a,b) -> Data.Random.pdf (Random.Normal (a*x + b) 1) y)

points :: [Main.Point] -> Dist Param -> Dist Param
points ps d = foldl (flip point) d ps

trainingData :: IO [(Double, Double)]
trainingData =
  let points = [-5..5]
      equation = \x -> 2 * x + 1
      addNoise x = Random.sample $ Random.Normal x 1
   in sequence . map (\x -> (,) x <$> addNoise (equation x)) $ points

main :: IO ()
main = do
  td <- trainingData

  putStrLn "Sampling..."
  xs <- Random.sample $ mh 100 100000 (points td linear)

  putStrLn "Drawing..."
  r2AxisMain $ myaxis $ nub xs
    where
      myaxis :: [(Double, Double)] -> Axis B V2 Double
      myaxis mydata1 = r2Axis &~ do
        scatterPlot mydata1 $ key "data 1"

• 実行

stack exec -- GADTtest-exe -o out.png