<jupyter_start><jupyter_text>a gentle exploration tour through the Monad-Bayes libraryhttps://github.com/adscib/monad-bayes<jupyter_code>:e OverloadedStrings

import qualified Graphics.Vega.VegaLite as VL
import Graphics.Vega.VegaLite hiding (filter)
import IHaskell.Display.Hvega (vlShow)
import Data.Aeson (ToJSON(toJSON), Value)
import Data.Aeson (encode)
import qualified Data.ByteString.Lazy.Char8 as BL

import Control.Monad (liftM2, replicateM)
import Control.Monad.IO.Class (liftIO)
import Control.Monad.Bayes.Class
import Control.Monad.Bayes.Sampler
import Control.Monad.Bayes.Inference.RMSMC as RMSMC
import Control.Monad.Bayes.Sequential
import Control.Monad.Bayes.Population
import Control.Monad.Bayes.Traced.Static (Traced)
import Control.Monad.Bayes.Inference.SMC
import Data.Vector (fromList, toList)
import Statistics.Sample.KernelDensity (kde)


import Numeric.Log
import Control.Monad.Bayes.Class


barPlot :: [Double] -> VegaLite
barPlot dataSamples = 
  let desc = description "Probability Games"
      dat = dataFromColumns []
          . dataColumn "Result" (Numbers dataSamples)
      trans = transform
        . binAs [ MaxBins 20 ] "Result" "binnedData"
      conf = configure
        . configuration (Axis [ DomainWidth 1 ])
        . configuration (BarStyle [ MContinuousBandSize 20 ])
        . configuration (SelectionStyle [ ( Single, [ On "dblclick" ] ) ])
      barOpts = [MOpacity 1.0, MColor "blue"]
      enc = encoding
            . position X [PName "Result", PmType Quantitative, PAxis [AxGrid True, AxTitle "Toss Value"]]
            . position Y [PName "binnedData", PmType Quantitative, PAxis [AxGrid False]]
  in toVegaLite [desc, width 600, height 200, dat [], trans [], conf [], mark Bar barOpts, enc []]
  
linePlot :: [Double] -> [Double] -> VegaLite
linePlot xSamples ySamples = 
  let desc = description "Logistic Regression"
      dat = dataFromColumns []
          . dataColumn "X" (Numbers xSamples)
          . dataColumn "Y" (Numbers ySamples)
      conf = configure
        . configuration (Axis [ DomainWidth 1 ])
        . configuration (SelectionStyle [ ( Single, [ On "dblclick" ] ) ])
      enc = encoding
            . position X [PName "X", PmType Quantitative, PAxis [AxGrid True, AxTitle "Toss Value"]]
            . position Y [PName "Y", PmType Quantitative, PAxis [AxGrid False]]
  in toVegaLite [desc, width 600, height 200, dat [], conf [], mark Line [], enc []]

pointPlot :: [Double] -> [Double] -> VegaLite
pointPlot xSamples ySamples = 
  let desc = description "Logistic Regression"
      dat = dataFromColumns []
          . dataColumn "X" (Numbers xSamples)
          . dataColumn "Y" (Numbers ySamples)
      conf = configure
        . configuration (Axis [ DomainWidth 1 ])
        . configuration (SelectionStyle [ ( Single, [ On "dblclick" ] ) ])
      enc = encoding
            . position X [PName "X", PmType Quantitative, PAxis [AxGrid True, AxTitle "Toss Value"]]
            . position Y [PName "Y", PmType Quantitative, PAxis [AxGrid False]]
  in toVegaLite [desc, width 600, height 200, dat [], conf [], mark Point [MOpacity 0.2], enc []]
  
density2DPlot :: [Double] -> [Double] -> VegaLite
density2DPlot xSamples ySamples = 
  let desc = description "Logistic Regression"
      dat = dataFromColumns []
          . dataColumn "X" (Numbers xSamples)
          . dataColumn "Y" (Numbers ySamples)
      conf = configure
        . configuration (Axis [ DomainWidth 1 ])
        . configuration (SelectionStyle [ ( Single, [ On "dblclick" ] ) ])
        . configuration (View [Stroke (Just "transparent")])
      enc = encoding
            . position X [PName "X", PBin [MaxBins 20], PmType Quantitative, PAxis [AxGrid True, AxTitle "Toss Value"]]
            . position Y [PName "Y", PBin [MaxBins 20], PmType Quantitative, PAxis [AxGrid False]]
            . color [ MAggregate Count, MName "col", MmType Quantitative, MScale [SReverse False, SScheme "blues"[0.0, 1.0]]]
  in toVegaLite [desc, width 300, height 300, dat [], conf [], mark Rect [], enc []]<jupyter_output><empty_output><jupyter_text>A Gaussian Toy Model<jupyter_code>gaussian :: MonadInfer m => [Double] -> m (Double, Double)
gaussian os = do
  c <- uniform (-0.5) 0.5
  s <- uniform 0.5 1.0
  mapM_ (score . normalPdf c s) os
  return (c, s)
  
generateDataGaussian :: MonadSample m
      => Int  -- ^ Number of points
      -> m [Double] -- ^ List of latent and observable states from t=1
generateDataGaussian n = replicateM n $ normal 0 0.1
runModel :: IO [((Double, Double), Log Double)]
runModel = sampleIOfixed $ do
  ys <- generateDataGaussian 50
  res <- runPopulation $ (RMSMC.rmsmc 100 10 10) (gaussian ys)
  return res
  
runModel<jupyter_output><empty_output><jupyter_text>Modelling Dice<jupyter_code>-- A toy model for dice rolling from http://dl.acm.org/citation.cfm?id=2804317
-- Exact results can be obtained using Dist monad

-- | A toss of a six-sided die.
die :: MonadSample m => m Int
die = uniformD [1..6]

-- | A sum of outcomes of n independent tosses of six-sided dice.
dice :: MonadSample m => Int -> m Int
dice 1 = die
dice n = liftM2 (+) die (dice (n-1))
-- | Toss of two dice where the output is greater than 4.
dice_hard :: MonadInfer m => m Int
dice_hard = do
  result <- dice 2
  condition (result > 8)
  return result
takeFromModel :: Int -> Sequential (Population SamplerIO) a -> IO [a]
takeFromModel nsamples inModel = sampleIOfixed $ do
    res <- runPopulation $ smcMultinomial 1 nsamples inModel
    return $ map fst res
samples <- takeFromModel 100 (dice 3)
(bins, counts) = kde 10 (fromList (fromIntegral <$> samples))
vlShow $ linePlot (toList bins) (toList counts)
samples <- takeFromModel 100 dice_hard
(bins, counts) = kde 50 (fromList (fromIntegral <$> samples))
vlShow $ linePlot (toList bins) (toList counts)<jupyter_output><empty_output><jupyter_text>Two alternating distributions<jupyter_code>prior :: MonadSample m => m Double
prior = return 0.2

modelFalse :: MonadSample m => m Int
modelFalse = uniformD ([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12] :: [Int])

modelTrue :: MonadSample m => m Int
modelTrue = do
  x <- bernoulli 0.5
  y <- uniformD (if x then [1, 2] else [11, 12] :: [Int])
  return y
samples <- takeFromModel 20 modelTrue
vlShow $ barPlot $ fromIntegral <$> samples<jupyter_output><empty_output><jupyter_text>Logistic Regression<jupyter_code>xs = [-10, -5, 2, 6, 10]
xsDouble = fromInteger <$> xs
labels = [False, False, True, True, True]
ysDouble = fromIntegral <$> fromEnum <$> labels
vlShow $ linePlot xsDouble ysDouble<jupyter_output><empty_output><jupyter_text>both normal and gamma are functions that take two parameters and that return distributions with the _MonadSample_ interface:<jupyter_code>:t normal
:t gamma<jupyter_output><empty_output><jupyter_text>we now build a sigmoid function from those. First we build a scaled and shifted normal distribution that is then passed through the exponential.<jupyter_code>sigmoid :: MonadSample m => Double -> m Double
sigmoid x = do
    m <- normal 0 1
    b <- normal 0 1
    sigma <- gamma 1 1
    n <- normal (m * x + b) sigma
    return (1 / (1 + exp (- n)))<jupyter_output><empty_output><jupyter_text>we can take 100 samples of the final distribution at a few x values and plot the resulting model distribution ...<jupyter_code>nsamples = 100
xValues = [-10..10]
yValues <- sequence (takeFromModel nsamples <$> (sigmoid <$> xValues))
xPoints = xValues >>= replicate nsamples
yPoints = mconcat yValues
vlShow $ pointPlot xPoints yPoints<jupyter_output><empty_output><jupyter_text>here is the distribution as a density plot<jupyter_code>vlShow $ density2DPlot xPoints yPoints<jupyter_output><empty_output><jupyter_text>the sigmoidal structure of the distribution becomes quite clear. The goal of the inference algorithm will be to find parameters that correspond to the observed values. We start with the definition of a helper function `obs` that returns an empty instance of the MonadInfer class. This class can sample and score different execution paths.<jupyter_code>obs :: (MonadInfer m) => Double -> Bool -> m ()
obs x label = do
      p <- sigmoid x -- get the y Value (True=1, False=0) of the sigmoid at point x
      score $ (Exp . log) $ if label then p else 1 - p
logisticRegression :: (MonadInfer m) => [(Double, Bool)] -> m Double
logisticRegression dat = do
  mapM_ (uncurry obs) dat
  sigmoid 10

model = logisticRegression (zip xs labels)
samples <- takeFromModel 100 model
(bins, counts) = kde 10 (fromList samples)
vlShow $ linePlot (toList bins) (toList counts)
xValues = xs
yValues <- takeFromModel nsamples model
vlShow $ pointPlot xValues yValues<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10863798869795871, 'fraction_numerical': 0.01833698535347711, 'mean_word_length': 3.7233115468409586, 'pattern_counts': {'":': 0, '<': 68, '<?xml version=': 0, '>': 79, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '90362', 'license': 'MIT', 'n_tokens_mistral': 2976, 'n_tokens_neox': 2728, 'n_words': 937, 'path': 'example/Haskell/bayesMonad/monad_bayes_intro.ipynb', 'repo_name': 'marcus7070/jupyterWith', 'size': '8671'}