// Dependencies
#r "nuget: XPlot.Plotly"
#r "nuget: MathNet.Numerics"
#r "nuget: Newtonsoft.Json"

// Imports 
open XPlot.Plotly
open MathNet.Numerics
open Newtonsoft.Json

type ParsedRandomVariable = 
    { Name             : string
      Conditionals     : string list
      Distribution     : string
      Parameters       : string list
      Observed         : string option }

type ParsedBayesianModel = ParsedRandomVariable list

open System

// Format: RVName [|Conditionals] ~ Distribution(Parameters) [: observed] 
// [] -> optional
// NOTE: There can't be any spaces in the distribution parameters.
let parseLineOfModel (lineInModel : string) : ParsedRandomVariable = 
 
    // Helper fn to split the string based on a variety of type of delimiters.
    // Resultant type is a list of strings to feed in for the pattern matching.
    let splitToList (toSplit : string) (delimiters : obj) : string list = 
        let split = 
            match delimiters with
            | :? string        as s   -> toSplit.Split(s, StringSplitOptions.RemoveEmptyEntries) 
            | :? array<string> as arr -> toSplit.Split(arr, StringSplitOptions.RemoveEmptyEntries) 
            | :? array<char>   as arr -> toSplit.Split(arr, StringSplitOptions.RemoveEmptyEntries) 
            | _ -> failwithf "Splitting based on delimiters failed as it is neither a string nor an array of strings: Of Type: %A - %A" (delimiters.GetType()) toSplit
        
        Array.toList split

    match splitToList lineInModel " " with
    | nameAndConditionals :: "~" :: distributionParametersObserved ->
        // Get the name and conditionals.
        let splitNameAndConditionals = splitToList nameAndConditionals "|"
        let name = splitNameAndConditionals.[0]
        let conditionals = 
            match splitNameAndConditionals with 
            | name :: conditionals -> 
                if conditionals.Length > 0 then splitToList conditionals.[0] ","
                else []
            | _ -> failwithf "Pattern not found for RV Name and Conditionals - the format is: RVName|Condtionals: %A" splitNameAndConditionals

        let extractAndGetParameters (distributionNameAndParameters : string) : string * string list = 
            let splitDistributionAndParameters = splitToList distributionNameAndParameters [| "("; ")" |]
            (splitDistributionAndParameters.[0], splitToList splitDistributionAndParameters.[1] ",")
            
        match distributionParametersObserved with 

        // Case: Without Observations. Example: θ ~ Gamma(a,b)
        | distributionNameAndParameters when distributionNameAndParameters.Length = 1 ->
            let extractedDistributionAndParameters = extractAndGetParameters distributionNameAndParameters.[0]
            { Name             = name; 
              Conditionals     = conditionals; 
              Distribution     = (fst extractedDistributionAndParameters).ToLower();
              Observed         = None; 
              Parameters       = snd extractedDistributionAndParameters; }

        // Case: With Observations. Example: Y|θ ~ Poisson(θ) : observed
        | distributionNameAndParameters :: ":" :: observed ->
            let extractedDistributionAndParameters = extractAndGetParameters distributionNameAndParameters
            { Name             = name;
              Conditionals     = conditionals; 
              Distribution     = (fst extractedDistributionAndParameters).ToLower();
              Observed         = Some observed.Head; // Only 1 observed list permitted.
              Parameters       = snd extractedDistributionAndParameters; } 

        // Case: Error.
        | _ -> failwithf "Pattern not found for the model while parsing the distribution, parameters and optionally, the observed variables: %A" distributionParametersObserved 

    | _ -> failwithf "Pattern not found for the following line in the model - please check the syntax: %A" lineInModel

let parseModel (model : string) : ParsedBayesianModel = 
    model.Split('\n') 
    |> Array.map(parseLineOfModel)
    |> Array.toList

// Helper to print out the Parsed Bayesian Model
let printParsedModel (model : string) : unit = 
    let parsedModel = parseModel model
    printfn "Model: %A is represented as %A" model parsedModel

// Print out our simple 1-Parameter Model.
let model1 = @"θ ~ Gamma(a,b)
              Y|θ ~ Poisson(θ) : observed"
printParsedModel(model1)

// This model doesn't make sense but adding to test multiple conditionals.
let model2  = @"θ ~ Beta(unit,unit)
               γ ~ Gamma(a,b)
               Y|θ,γ ~ Binomial(n,θ) : observed"
printParsedModel(model2)

Model: "θ ~ Gamma(a,b)
              Y|θ ~ Poisson(θ) : observed" is represented as [{ Name = "θ"
   Conditionals = []
   Distribution = "gamma"
   Parameters = ["a"; "b"]
   Observed = None }; { Name = "Y"
                        Conditionals = ["θ"]
                        Distribution = "poisson"
                        Parameters = ["θ"]
                        Observed = Some "observed" }]
Model: "θ ~ Beta(unit,unit)
               γ ~ Gamma(a,b)
               Y|θ,γ ~ Binomial(n,θ) : observed" is represented as [{ Name = "θ"
   Conditionals = []
   Distribution = "beta"
   Parameters = ["unit"; "unit"]
   Observed = None }; { Name = "γ"
                        Conditionals = []
                        Distribution = "gamma"
                        Parameters = ["a"; "b"]
                        Observed = None }; { Name = "Y"
                                             Conditionals = ["θ"; "γ"]
                                             Distribution = "binomial"
                                             Parameters = ["n"; "θ"]
                                             Observed = Some "observed" }]

open System
open System.Collections.Generic

open MathNet.Numerics.Distributions

open Newtonsoft.Json

type Observed = float list

type ParameterList = 
    { Observed : float list; Parameters : Dictionary<string, float> } 

let deserializeParameters (paramsAsString : string) : ParameterList = 
    JsonConvert.DeserializeObject<ParameterList>(paramsAsString)

// Parameter List 1
let parameters1 = "{Parameters : {μ0 : 0, σ0 : 1, μ : 5, σ : 2, λ : 4}, observed : [4.2,0.235,2.11]}"
let deserializedParameters1 = deserializeParameters parameters1
printfn "Deserialized Parameters 1: %A" (deserializedParameters1)

// Parameter List 2
let parameters2 = "{Parameters: {λ : 2}}"
let deserializedParameters2 = deserializeParameters parameters2
// Applying the Deserialized Parameters to Sample from a Distribution
let exp = Exponential deserializedParameters2.Parameters.["λ"] 
printfn "Sampling from the Exponential Distribution with the λ = %A parameter: %A" exp (exp.Sample())

Deserialized Parameters 1: { Observed = [4.2; 0.235; 2.11]
  Parameters = seq [[μ0, 0]; [σ0, 1]; [μ, 5]; [σ, 2]; ...] }
Sampling from the Exponential Distribution with the λ = Exponential(λ = 2) parameter: 1.673949327

open MathNet.Numerics.Distributions

type DistributionType = 
    | Continuous
    | Discrete

type Parameter     = float
type DiscreteInput = int
type Input         = float

type DensityType = 
    | OneParameter         of (Parameter * Input -> float) 
    | OneParameterDiscrete of (Parameter * DiscreteInput -> float)
    | TwoParameter         of (Parameter * Parameter * Input -> float)

type DistributionInfo = { RVName             : string
                          DistributionType   : DistributionType 
                          DistributionName   : string 
                          Parameters         : float list
                          Density            : DensityType } with
                          
    static member Create (item : ParsedRandomVariable) 
                         (parameterList : ParameterList) : DistributionInfo = 
    // I know this is ugly but this functionality assumes the user enters the 
    // parameters in the order that's expected by the MathNet Numerics Library. 
    // Grab the parameters associated with this Random Variable.
    let rvParameters = 
        item.Parameters
        |> List.filter(parameterList.Parameters.ContainsKey) 
        |> List.map(fun item -> parameterList.Parameters.[item])

    // Extract Distribution Associated with the Parsed Random Variable.
    match item.Distribution with

    // 1 Parameter Distributions 
    | "exponential" -> 
        { RVName           = item.Name
          DistributionName = item.Distribution
          Parameters       = rvParameters
          DistributionType = DistributionType.Continuous
          Density          = OneParameter Exponential.PDF }
    | "poisson" ->
        { RVName           = item.Name
          DistributionName = item.Distribution
          Parameters       = rvParameters
          DistributionType = DistributionType.Discrete
          Density          = OneParameterDiscrete Poisson.PMF } 
    // 2 Parameter Distributions 
    | "normal" ->
        { RVName           = item.Name
          DistributionName = item.Distribution
          Parameters       = rvParameters
          DistributionType = DistributionType.Continuous
          Density          = TwoParameter Normal.PDF }
    | "gamma" ->
        { RVName           = item.Name
          DistributionName = item.Distribution
          Parameters       = rvParameters
          DistributionType = DistributionType.Continuous
          Density          = TwoParameter Gamma.PDF }
    | "beta" ->
        { RVName           = item.Name
          DistributionName = item.Distribution
          Parameters       = rvParameters
          DistributionType = DistributionType.Continuous
          Density          = TwoParameter Beta.PDF }
    | "continuousuniform" ->
        { RVName           = item.Name
          DistributionName = item.Distribution
          Parameters       = rvParameters
          DistributionType = DistributionType.Continuous
          Density          = TwoParameter ContinuousUniform.PDF }
    // Failure Case
    | _ -> failwithf "Distribution not registered: %A" item.Distribution                      
                          
    member this.ComputeOneParameterPDF (parameter : float) (input : float) : float =
        match this.Density with
        | OneParameter pdf -> pdf(parameter,input)
        | _ -> failwithf "Incorrect usage of function with a non One Parameter Density Type. Density given: %A" this.Density 
    member this.ComputeOneParameterDiscretePMF (parameter : float) (input : int) : float =
        match this.Density with
        | OneParameterDiscrete pmf -> pmf(parameter,input)
        | _ -> failwithf "Incorrect usage of function with a non One Parameter Discrete Density Type. Density given: %A" this.Density 
    member this.ComputeTwoParameterPDF (parameter1 : float) (parameter2 : float) (input : float) : float =
        match this.Density with
        | TwoParameter pdf -> pdf(parameter1, parameter2, input)
        | _ -> failwithf "Incorrect usage of function with a non Two Parameter Density Type. Density given: %A" this.Density

let getDistributionInfoForModel(model : string) (parameterList : string) : DistributionInfo list = 
    let parsedModel   = parseModel model
    let parameterList = deserializeParameters parameterList 
    parsedModel
    |> List.map(fun x -> DistributionInfo.Create x parameterList)

let getDensityOrProbabilityForModel (model : string) 
                                    (parameterList : string) 
                                    (data : float seq) : IDictionary<string, float seq> = 
    getDistributionInfoForModel model parameterList
    |> List.map(fun (e : DistributionInfo) -> 
        match e.Density with
        | OneParameter p ->
            let param   = List.exactlyOne e.Parameters 
            let results = data |> Seq.map(fun d -> e.ComputeOneParameterPDF param d)
            e.RVName, results
        | OneParameterDiscrete p ->
            let param   = List.exactlyOne e.Parameters 
            let results = data |> Seq.map(fun d -> e.ComputeOneParameterDiscretePMF param (int d))
            e.RVName, results
        | TwoParameter p ->
            let p2 : float list = e.Parameters |> List.take 2
            let results = data |> Seq.map(fun d -> e.ComputeTwoParameterPDF p2.[0] p2.[1] d)
            e.RVName, results)
    |> dict

// Exponential. 
let exponentialModel     = "x ~ Exponential(lambda)"
let exponentialParamList = "{Parameters: {lambda : 2, a : 2., b : 2.3}, Observed : []}"
let exponentialDummyData = ContinuousUniform.Samples(0., 200.) |> Seq.take 2000
let exponentialPdfs      = getDensityOrProbabilityForModel exponentialModel exponentialParamList exponentialDummyData
printfn "Exponential: %A" (exponentialPdfs.Values)

// Normal.
let normalModel     = "x ~ Normal(mu,sigma)"
let normalParamList = "{Parameters: {mu: 0., sigma : 1.}, Observed : []}"
let normalDummyData = Normal.Samples(0.0, 1.0) |> Seq.take 2000
let normalPdfs      = getDensityOrProbabilityForModel normalModel normalParamList normalDummyData 
printfn "Normal: %A" (normalPdfs.Values)

// Poisson.
let poissonModel     = "x ~ Poisson(theta)"
let poissonParamList = "{Parameters: {theta: 44}, Observed : []}"
let poissonDummyData = ContinuousUniform.Samples(0., 5.) |> Seq.take 2000 
let poissonPdfs      = getDensityOrProbabilityForModel poissonModel poissonParamList poissonDummyData 
printfn "Poisson: %A" (poissonPdfs.Values)

Exponential: seq
  [seq [5.021224175e-54; 6.790246353e-50; 2.044679134e-68; 2.02863983e-172; ...]]
Normal: seq [seq [0.3813916696; 0.1382035578; 0.1331370548; 0.2884060033; ...]]
Poisson: seq
  [seq [1.104713281e-15; 1.21518461e-14; 7.532136009e-17; 1.21518461e-14; ...]]

type BayesianNodeTypeInfo =
    | Observed of float list 
    | NonObserved

type BayesianNode = 
    { Name                 : string
      NodeType             : BayesianNodeTypeInfo 
      DistributionInfo     : DistributionInfo
      ParsedRandomVariable : ParsedRandomVariable } with

    static member Create(parsedRandomVariable : ParsedRandomVariable)
                        (parameterList : ParameterList) =

        let nodeType : BayesianNodeTypeInfo =
            match parsedRandomVariable.Observed with
            | Some _ -> BayesianNodeTypeInfo.Observed parameterList.Observed
            | None   -> BayesianNodeTypeInfo.NonObserved
        { Name                 = parsedRandomVariable.Name;
          NodeType             = nodeType;
          DistributionInfo     = DistributionInfo.Create parsedRandomVariable parameterList;
          ParsedRandomVariable = parsedRandomVariable; }

    member this.GetDependents (parsedBayesianModel : ParsedBayesianModel) : ParsedRandomVariable list =
        parsedBayesianModel 
        |> List.filter(fun x -> x.Conditionals |> List.contains(this.Name))

let lineOfModel = @"x ~ Exponential(lambda) : observed"
let paramList   = "{Parameters: {lambda : 2, a : 2., b : 2.3 }, observed : [1,2,3,55]}"
BayesianNode.Create (parseLineOfModel lineOfModel) (deserializeParameters paramList)

// Only to be used for a model with 2 nodes 
// i.e. one for the prior and one for the likelihood.
type SimpleBayesianNetworkModel = 
  { Name       : string
    Nodes      : IDictionary<string, BayesianNode>
    Prior      : BayesianNode
    Likelihood : BayesianNode } with

    member this.GetPriorProbability (input : float) : float = 
      let distributionInfo = this.Prior.DistributionInfo 
      match distributionInfo.Density with
      | OneParameter p ->
          let param = List.exactlyOne distributionInfo.Parameters 
          distributionInfo.ComputeOneParameterPDF param input 
      | OneParameterDiscrete p ->
          let param   = List.exactlyOne distributionInfo.Parameters 
          distributionInfo.ComputeOneParameterDiscretePMF param (int input)
      | TwoParameter p ->
          let p2 : float list = distributionInfo.Parameters |> List.take 2
          distributionInfo.ComputeTwoParameterPDF p2.[0] p2.[1] input

    member this.GetLikelihoodProbability (prior : float) : float = 
        let distributionInfo = this.Likelihood.DistributionInfo
        let observed         = match this.Likelihood.NodeType with
                               | Observed l -> l
                               | _ -> failwithf "Incorrectly constructed Simple Network Model. %A" this 
        let density = 
          match distributionInfo.Density with
          | OneParameter p -> 
              observed |> List.map(fun d -> distributionInfo.ComputeOneParameterPDF prior d)
          | OneParameterDiscrete p ->
              observed |> List.map(fun d -> distributionInfo.ComputeOneParameterDiscretePMF prior (int (Math.Ceiling d)))
          | TwoParameter p ->
              let p : float = distributionInfo.Parameters |> List.exactlyOne 
              observed |> List.map(fun d -> distributionInfo.ComputeTwoParameterPDF prior p d)
        density
        |> List.fold (*) 1.0

    member this.GetPosteriorWithoutScalingFactor (input: float) : float = 
      let priorPdf      = this.GetPriorProbability input 
      let likelihoodPdf = this.GetLikelihoodProbability priorPdf
      priorPdf * likelihoodPdf

    static member Construct (name : string)
                            (model : ParsedBayesianModel)
                            (parameterList : ParameterList) : SimpleBayesianNetworkModel = 

        // Construct all the modes of the model.
        let allNodes : (string * BayesianNode) list =
            model
            |> List.map(fun m -> m.Name, BayesianNode.Create m parameterList)

        let prior : BayesianNode =  
            allNodes
            |> List.filter(fun (_,m) -> m.NodeType = BayesianNodeTypeInfo.NonObserved) 
            |> List.map(fun (_,m) -> m)
            |> List.exactlyOne

        let likelihood : BayesianNode =  
            allNodes
            |> List.filter(fun (_,m) -> m.NodeType <> BayesianNodeTypeInfo.NonObserved)
            |> List.map(fun (_,m) -> m)
            |> List.exactlyOne

        { Name       = name; 
          Nodes      = dict allNodes; 
          Prior      = prior;
          Likelihood = likelihood; }

let model = @"x ~ Normal(μ,τ) 
              y|x ~ Normal(x,σ) : observed"
let parsedModel = parseModel model
let paramList = "{Parameters: {μ : 5, τ : 3.1622, σ : 1}, observed : [9.37,10.18,9.16,11.60,10.33]}"
let normalModel = 
    SimpleBayesianNetworkModel.Construct "Normal-Normal" parsedModel (deserializeParameters paramList)
normalModel

let model = @"x ~ Normal(μ,τ) 
              y|x ~ Normal(x,σ) : observed"
let parsedModel = parseModel model
let paramList = "{Parameters: {μ : 5, τ : 3.1622, σ : 1}, observed : [9.37,10.18,9.16,11.60,10.33]}"
let normalModel = 
    SimpleBayesianNetworkModel.Construct "Normal-Normal" parsedModel (deserializeParameters paramList)
    
// Prior Check
let priorFromModel = normalModel.GetPriorProbability 0.0
let computedPrior  = Normal.PDF(5., 3.1622, 0.0)
printfn "Normal Model Prior: %A | Prior from formula: %A" priorFromModel computedPrior

// Likelihood Check
let likelihoodFromModel = normalModel.GetLikelihoodProbability 0.0
let observed : float list = 
    match normalModel.Likelihood.NodeType with
    | Observed x -> x
    | _ -> failwithf "Incorrectly constructed Likelihood"
    
let likelihoodListFromFormula : float list = 
    observed |> List.map(fun (o:float) -> Normal.PDF(computedPrior, 1., o))
let likelihoodFromFormula = likelihoodListFromFormula |> List.fold (*) 1.

printfn "Normal Model Likelihood: %A | Likelihood from formula: %A" 
    (normalModel.GetLikelihoodProbability 0.03614314702) 
    (likelihoodFromFormula)

Normal Model Prior: 0.03614314702 | Prior from formula: 0.03614314702
Normal Model Likelihood: 4.158402902e-114 | Likelihood from formula: 4.158402901e-114

type ConvergenceCriteria = 
    | IterativeConvergence of int // Number of Iterations

type ProposalDistribution =
    | Normal of float  // Normal( current, delta )
    | ContinuousUniform of float // ContinuousUniform( current - delta, current + delta )
    | PositiveContinuousUniform of float // ( x ~ Uniform( current - delta, current + delta ) if x <= 0 then 0.1 else x )

type MCMCInferenceStep = 
    | SymmetricMetropolisHastings of ProposalDistribution * float

type MCMCChain =
    { Id             : int
      AcceptanceRate : float // # of Acceptances / Total # of Chain Items
      StepValues     : float seq }

type MCMCRequest = 
    { StepFunction         : MCMCInferenceStep
      ConvergenceCriteria  : ConvergenceCriteria
      BurnInIterationsPct  : float 
      Chains               : int 
      PrintDebug           : bool }

type MCMCResult =
    { Chains         : MCMCChain seq 
      MCMCRequest    : MCMCRequest }

open System
open MathNet.Numerics.Distributions
open XPlot.Plotly

type AcceptanceRejection =
    | Acceptance of float // Acceptance of the Proposal
    | Rejection  of float // Rejection of the Proposal

let doSymmetricMetropolisHastings (request : MCMCRequest) 
                                  (iterations : int)
                                  (proposalDistribution : ProposalDistribution)
                                  (initialValue : float)
                                  (simpleBayesianModel : SimpleBayesianNetworkModel) : MCMCResult =

    let getChain (id : int) (request : MCMCRequest) 
                 (iterations : int) (simpleBayesianModel : SimpleBayesianNetworkModel) : MCMCChain =

        let burnin         = int (Math.Ceiling(request.BurnInIterationsPct / 100. * float iterations))
        let zeroOneUniform = ContinuousUniform()

        let mutable current = simpleBayesianModel.GetPosteriorWithoutScalingFactor initialValue

        let matchedProposalDistribution (input : float) :  float =
            match proposalDistribution with 
            | Normal delta            -> Normal(input, delta).Sample()
            | ContinuousUniform delta -> ContinuousUniform(input - delta, input + delta).Sample()
            | PositiveContinuousUniform delta -> 
                let u = ContinuousUniform(input - delta, input + delta).Sample()
                if u <= 0. then input + 0.1 else u

        let step (iteration : int) : AcceptanceRejection =
            let proposed           = matchedProposalDistribution current
            let currentProbability = simpleBayesianModel.GetPosteriorWithoutScalingFactor current
            let proposedProbabilty = simpleBayesianModel.GetPosteriorWithoutScalingFactor proposed
            let acceptanceRatio    = Math.Min(currentProbability / proposedProbabilty, 1.)
            let uniformDraw        = zeroOneUniform.Sample()

            if request.PrintDebug then 
                printfn "Chain: %A Iteration: %A - Current Probability: %A | Proposed Probability: %A | AcceptanceRatio: %A | Uniform Draw: %A | Current: %A | Proposed: %A" 
                    id iteration currentProbability proposedProbabilty acceptanceRatio uniformDraw current proposed

            if uniformDraw < acceptanceRatio then (current <- proposed; Acceptance proposed) 
            else Rejection current
    
        let stepResults : AcceptanceRejection seq =
            seq {1..iterations}
            |> Seq.map step
            |> Seq.skip burnin

        let getAcceptanceRateAndStepValues : float * float seq =

            // Compute the Acceptance Rate.
            let acceptanceRate : float = 
                let totalBurninWithoutBurnin : float = float (Seq.length stepResults)
                let totalNumberOfAcceptances : float = 
                    stepResults 
                    |> Seq.filter(fun x -> 
                        match x with
                        | Acceptance x -> true
                        | _ -> false)
                    |> Seq.length
                    |> float
                totalNumberOfAcceptances / totalBurninWithoutBurnin

            // Grab the Step Values that'll approximate the posterior.
            let stepValues =
                stepResults |> Seq.map(fun s ->
                    match s with
                    | Acceptance v -> v
                    | Rejection v -> v)

            acceptanceRate, stepValues

        let acceptanceRate, stepValues = getAcceptanceRateAndStepValues 

        { Id             = id
          AcceptanceRate = acceptanceRate
          StepValues     = stepValues} 

    let chains : MCMCChain seq =
        seq {1..request.Chains}
        |> Seq.map(fun id ->  getChain id request iterations simpleBayesianModel)

    { Chains      = chains
      MCMCRequest =  request }

let runMCMC (request : MCMCRequest) 
            (model : SimpleBayesianNetworkModel) : MCMCResult =
    match request.StepFunction with
    | SymmetricMetropolisHastings (proposalDistribution,initialValue) ->
        match request.ConvergenceCriteria with 
        | IterativeConvergence iterations -> 
            doSymmetricMetropolisHastings request iterations proposalDistribution initialValue model
        | _ -> failwith "You need to pass in the number of iterations for the Metropolis-Hastings algorithm"
    | _ -> failwithf "Step Function Not Registered: %A" request.StepFunction

let model = @"x ~ Exponential(a) 
              y|x ~ Exponential(x) : observed"
let parsedModel = parseModel model
let paramList = "{Parameters: {a : 2., b : 2.}, observed : [1,2,3,4,4,2,5,6,7,3,2,3,4,5,6,1,2,3,4,4,4,4]}"
let simpleModel = SimpleBayesianNetworkModel.Construct "Exponential Model" parsedModel (deserializeParameters paramList)

let request : MCMCRequest = 
    { StepFunction        = SymmetricMetropolisHastings (ProposalDistribution.PositiveContinuousUniform 0.002, 0.001)
      ConvergenceCriteria = IterativeConvergence 10000
      BurnInIterationsPct = 0.2
      PrintDebug          = false 
      Chains              = 4 }

let mcmc       = runMCMC request simpleModel
let firstChain = Seq.head mcmc.Chains

Histogram(x = firstChain.StepValues)
|> Chart.Plot
|> Chart.WithTitle "Exponential Prior and Likelihood"

#!fsharp
let model = @"x ~ Normal(μ,τ) 
              y|x ~ Exponential(x) : observed"
let parsedModel = parseModel model
let paramList = "{Parameters: {μ : 5, τ : 3.1622, σ : 1}, observed : [9.37,10.18,9.16,11.60,10.33]}"
let simpleModel = 
    SimpleBayesianNetworkModel.Construct "Normal-Exponential" parsedModel (deserializeParameters paramList)


let request : MCMCRequest = 
    { StepFunction        = SymmetricMetropolisHastings (ProposalDistribution.PositiveContinuousUniform 0.003, 0.0001)
      ConvergenceCriteria = IterativeConvergence 50000
      BurnInIterationsPct = 0.2
      PrintDebug          = false 
      Chains              = 4 }

let mcmc       = runMCMC request simpleModel
let firstChain = Seq.head mcmc.Chains

Histogram(x = firstChain.StepValues)
|> Chart.Plot
|> Chart.WithTitle "Normal Prior with Exponential Likelihood"
|> Chart.WithWidth 700
|> Chart.WithHeight 500