diff --git a/previews/PR150/BUGS_notes/index.html b/previews/PR150/BUGS_notes/index.html index f5572f392..4a9852a47 100644 --- a/previews/PR150/BUGS_notes/index.html +++ b/previews/PR150/BUGS_notes/index.html @@ -1,2 +1,2 @@ -Notes on BUGS Implementations · JuliaBUGS.jl
Edit on GitHub

Miscellaneous Notes on BUGS

Here are some exert from BUGS Developer Manual and notes on the original BUGS implementations.

Lexing

  • The BUGS language has the convention that if a name is followed immediately by a round bracket, that is by a "(", then the names is a reserved name in the BUGS language and does not represent a variable in the model.
  • By scanning the stream of tokens that constitute a BUGS language model the names of all the variables in the model can be found.

Table of Names

  • The BUGS language compiler expands all the for loops in the model and records the value of the indices of each use of a tensor on the left hand side of each relation.
  • The range of each index, for a tensor, is set at the maximum value observed value of the index and added to the name table. There is one exception to this procedure for finding index bounds: names that are data, that is in the data source, have the ranges of their indices fixed in the data source.
  • Each scalar and each component of a tensor used on the right hand side of a relation must occur either on the left hand side of a relation and or in a data source.

Data Transformations

If the compiler can prove that a logical assignment can be evaluated to a constant then the assignment is called a data transformation. This occurs if an assignment's right hand side does not depend on any variable quantities. The BUGS language has a general rule that there must only be one assignment statement for each scalar or component of a tensor. This rule is slightly relaxed for data transformations. The language allows a logical assignment and a stochastic assignment to the same scalar or tensor component if and only if the logical assignment is a data transformation.

Generated Quantities

Only need to be evaluated after the inference algorithm has finished its task.

  • Generally, these are leaf nodes that logical variables
  • In the case of stochastic variables that are leaf nodes, do “forward sampling”, also part of the generated Quantities

Computation

  • All the nodes in the graphical model representing logical relations are placed into an

array and sorted by their nesting level with the first array entries only depending on quantities defined by stochastic relations. Traversing this array and evaluating nodes gives up to date values to all logical relations.

Types

The BUGS compiler uses the properties of the distribution on the right-hand side of a stochastic assignment statement to make deductions about the variable on the left-hand side. For example, r ~ dbin(p, n) implies that r is integer-valued, while x ~ dnorm(mu, tau) implies that x is real-valued.

Some distributions are real-valued but have support on a restricted range of the reals. For example, p ~ dbeta(a, b) implies that p is real-valued with support on the unit interval, while x ~ dgamma(r, lambda) implies that x is real-valued but with support on the positive real line.

There are two multivariate distributions in the BUGS language, the Dirichlet and the Wishart, that have support on a complex subspace of the reals. The Dirichlet has support on the unit simplex, while the Wishart has support on symmetric positive definite matrices.

The BUGS compiler tries to infer if logical relations return an integer value by looking at whether their parents are integer-valued and the operators that combine the values of their parents into the return value. For example, in the cure model example above, the logical relation state1[i] <- state[i] + 1 is integer-valued because state[i] is a Bernoulli variable and therefore integer, the literal 1 is integer, and the sum of two integers is an integer.

When the BUGS system reads in data from a data source, it can tag whether the number read is an integer or a real and propagate this information to logical relations. Again, using the cure model as an example, the statement t[i] <- x[i] + y[i] is integer-valued because both x and y are data and are given as integers in the data source.

One special type of data is constants: that is just numbers with no associated distribution. Constants have many uses in BUGS language models, but one of the most important is as covariates. A model can contain a large number of constants that are used as covariates. Because of the possible large numbers of these covariate-type constants, they are given special treatment by the BUGS compiler. If a name read in from a data source is only used on the right-hand side of logical relations, no nodes in the graphical model are created to hold its values; they are directly incorporated in the objects that represent the right-hand sides of the logical relations.

For example, the large Methadone model contains the regression:

mu.indexed[i] <- beta[1] * x1[i] + beta[2] * x2[i] + beta[3] * x3[i] + beta[4] * x4[i] + region.effect[region.indexed[i]] + source.effect[region.indexed[i]] * source.indexed[i] + person.effect[person.indexed[i]]

where i ranges from 1 to 240776. Not having to create a node in the graphical model to represent x1, x2, x3, x4, region.indexed, source.index, and person.indexed saves a large amount of space.

In the BUGS language, the type information is fine-grained: each component of a tensor can have different type information. This is quite distinct from the situation in STAN and can make it much easier to specify a statistical model. One common case is where some components of a tensor have been observed while other components need to be estimated. The STAN documentation suggests workarounds for these situations, but these are somewhat complex.

  • The type propagation is interesting and maybe useful. But we don’t necessarily need to implement a type system. A dirty way to get type information is simply do a dry run with some tricks.

Work flow

The statistical model and data are presented to the BUGS system in a series of stages. In the first stage the model text is parsed into a tree and the name table constructed. The data is then loaded and checked against the model. The data can be split over a number of source. Once all the data has been loaded the model is compiled. Compiling builds the graphical model and does a large number of checks on the consistency of the model. Finally initial values can be given or generated for the model.

The compiler creates a node in the graphical model for each scalar name and each component of a tensor name in the BUGS language model. The compiler checks that only one node is created for each scalar name or component of a tensor name.

Reading in a data source causes the compiler to create special nodes called constant nodes to hold the values of the data.

The compiler processes logical relations before stochastic relations. Any logical relations that only have constant nodes on their right hand side become new constant nodes with the appropriate fixed value. Even if a logical relation can not be reduced to a constant some parts of the relation might be reduced to constants.

Any constant nodes that have an associated stochastic relation become data nodes in the graphical model.

Logical relations in the BUGS Language

The OpenBUGS software compiles a description of a statistical model in the BUGS language into a graph of objects. Each relation in the statistical model gives rise to a node in the graph of objects. Each distinct type of relation in the statistical model is represented by a node of a distinct class. For stochastic relations there is a fixed set of distributions that can be used in the modelling. For logical relations the situation is more complex. The software can use arbitrary logical expressions build out of a fixed set of basic operators and functions. For each distinct logical expression a new software source code module is written to implement a class to represent that logical expression in the graph of objects. The software module is then compiled using the Components Pascal compiler and the executable code merged into the running OpenBUGS software using the run time loading linker.

The BUGS language description of a statistical model is parsed into a list of trees. The sub-trees that represent logical relations in the statistical model are first converted into a stack based representation and then into Component Pascal source code. The source code is generated in module BugsCPWrite and the source code is then compiled in module BugsCPCompiler. Usually the generated source code is not displayed. Checking the Verbose option in the Info menu will cause each each source code module generated by the OpenBUGS software to be displayed in a separate window.

One advantage of a stack based representation of an expression is that it is straight forward to use it to derive source code that calculates the derivative of the expression with respect to its arguments. This part of the source code generation is carried out in module BugsCPWrite in procedure WriteEvaluateDiffMethod. Each operator in the stack representation of the logical expression causes a snippet of Component Pascal code to be written. These code snippets are generally very simple with those of binary operators slightly more complex than those of unitary operators. Each binary operators can emit three different code snippets: the general case and two special snippets depending on whether the left or right operands are numerical constants. The only complex code snippet is when an operand that is a logical relation in the statistical model is pushed onto the stack – the > case of nested logical relations. In this case the nested logical relation will have its own code to calculate derivatives and these values can be passed up the nesting level.

The OpenBUGS software now uses a backward mode scheme to calculate the value of logical nodes in the statistical model. All the logical nodes in the statistical model are held in a global array and sorted according to their nesting level with unnested nodes at the start of the array. To evaluate all the logical nodes in the statistical model this array is then traversed and each logical node evaluated and the value stored in the node. The same scheme is used to calculate derivatives.

The graphs derived from the BUGS language representation of statistical models are generally sparse. The OpenBUGS software uses conditional independence arguments to exploit sparsity in the stochastic parts of the model. There is also a sparsity structure in logical relations.Each logical relation will often depend on just a few stochastic parents and derivatives with respect to other stochastic nodes in the model will be structurally zero. Each logical node has an associated array of stochastic parents for which the derivatives are non zero. Moving up the level of nesting the number of parents can grow. Dealing with this issue leads to the complexity in the code snippet for the operator that pushes a logical node onto the stack. These issues can be seen in the non-linear random effects model called Orange trees in volume II of the OpenBUGS examples. In this model eta[i,] is a function of phi[i,1], phi[i,2] and phi[i,3] where the phi are also logical functions of the stochastic theta[i,].

One refinement of the backward mode scheme used to calculate the value of logical nodes is to consider separately any logical nodes in the statistical model which are only used for prediction and do not affect the calculation of the joint probability distribution. These nodes need only be evaluated once per iteration of the inference algorithm. Examples of such nodes are sigma[k] and sigma.C in the Orange trees example. There is no need to evaluate the derivatives of these prediction nodes.

The workings of the backward mode scheme are easy to visualize when the inference algorithm updates all the stochastic nodes in the statistical model in one block. Local versions of the backward mode scheme can be used when the inference algorithm works on single nodes or when a small blocks of nodes are updated. Each stochastic node is given its own vector of logical nodes that depend on it either directly or via other logical nodes and this vector is sorted by nesting level. Each updater that works on small blocks of nodes contains a vector of logical nodes which is the union of the vectors of dependent logical nodes for each of its components.

The idea of the backward mode scheme for evaluating logical nodes can be used with caching in Metropolis Hastings sampling. First the vector of logical nodes depending on the relevant stochastic node(s) is evaluated and their values cached. The log of the conditional distribution is then calculated. Next a new value of the stochastic node is proposed. The vector of logical nodes is re-evaluated and the log of the > conditional distribution calculated. If the proposed value is rejected then the cache is used to set the vector of logical nodes back to its old values.

The OpenBUGS software also calculates what class of function each logical node is in terms of its stochastic parents. If the software can prove for example that a logical node is a linear function of its parents more efficient sampling algorithms can be used. If a linear relation can be proved then the calculation of derivatives can also be optimized in some cases because they will be constant and so only need to be calculated once. Generalized linear models are implemented in a way that allows fast calculation of derivatives. The structure of the algorithm to classify the functional form of logical nodes is very similar to that for derivatives and uses a backward mode scheme

BUGS separates management of logical and stochastic variables, essentially two graphs. Logical variables are stored in an array and values are updated with values in earlier positions of the array.

+Notes on BUGS Implementations · JuliaBUGS.jl
Edit on GitHub

Miscellaneous Notes on BUGS

Here are some exert from BUGS Developer Manual and notes on the original BUGS implementations.

Lexing

  • The BUGS language has the convention that if a name is followed immediately by a round bracket, that is by a "(", then the names is a reserved name in the BUGS language and does not represent a variable in the model.
  • By scanning the stream of tokens that constitute a BUGS language model the names of all the variables in the model can be found.

Table of Names

  • The BUGS language compiler expands all the for loops in the model and records the value of the indices of each use of a tensor on the left hand side of each relation.
  • The range of each index, for a tensor, is set at the maximum value observed value of the index and added to the name table. There is one exception to this procedure for finding index bounds: names that are data, that is in the data source, have the ranges of their indices fixed in the data source.
  • Each scalar and each component of a tensor used on the right hand side of a relation must occur either on the left hand side of a relation and or in a data source.

Data Transformations

If the compiler can prove that a logical assignment can be evaluated to a constant then the assignment is called a data transformation. This occurs if an assignment's right hand side does not depend on any variable quantities. The BUGS language has a general rule that there must only be one assignment statement for each scalar or component of a tensor. This rule is slightly relaxed for data transformations. The language allows a logical assignment and a stochastic assignment to the same scalar or tensor component if and only if the logical assignment is a data transformation.

Generated Quantities

Only need to be evaluated after the inference algorithm has finished its task.

  • Generally, these are leaf nodes that logical variables
  • In the case of stochastic variables that are leaf nodes, do “forward sampling”, also part of the generated Quantities

Computation

  • All the nodes in the graphical model representing logical relations are placed into an

array and sorted by their nesting level with the first array entries only depending on quantities defined by stochastic relations. Traversing this array and evaluating nodes gives up to date values to all logical relations.

Types

The BUGS compiler uses the properties of the distribution on the right-hand side of a stochastic assignment statement to make deductions about the variable on the left-hand side. For example, r ~ dbin(p, n) implies that r is integer-valued, while x ~ dnorm(mu, tau) implies that x is real-valued.

Some distributions are real-valued but have support on a restricted range of the reals. For example, p ~ dbeta(a, b) implies that p is real-valued with support on the unit interval, while x ~ dgamma(r, lambda) implies that x is real-valued but with support on the positive real line.

There are two multivariate distributions in the BUGS language, the Dirichlet and the Wishart, that have support on a complex subspace of the reals. The Dirichlet has support on the unit simplex, while the Wishart has support on symmetric positive definite matrices.

The BUGS compiler tries to infer if logical relations return an integer value by looking at whether their parents are integer-valued and the operators that combine the values of their parents into the return value. For example, in the cure model example above, the logical relation state1[i] <- state[i] + 1 is integer-valued because state[i] is a Bernoulli variable and therefore integer, the literal 1 is integer, and the sum of two integers is an integer.

When the BUGS system reads in data from a data source, it can tag whether the number read is an integer or a real and propagate this information to logical relations. Again, using the cure model as an example, the statement t[i] <- x[i] + y[i] is integer-valued because both x and y are data and are given as integers in the data source.

One special type of data is constants: that is just numbers with no associated distribution. Constants have many uses in BUGS language models, but one of the most important is as covariates. A model can contain a large number of constants that are used as covariates. Because of the possible large numbers of these covariate-type constants, they are given special treatment by the BUGS compiler. If a name read in from a data source is only used on the right-hand side of logical relations, no nodes in the graphical model are created to hold its values; they are directly incorporated in the objects that represent the right-hand sides of the logical relations.

For example, the large Methadone model contains the regression:

mu.indexed[i] <- beta[1] * x1[i] + beta[2] * x2[i] + beta[3] * x3[i] + beta[4] * x4[i] + region.effect[region.indexed[i]] + source.effect[region.indexed[i]] * source.indexed[i] + person.effect[person.indexed[i]]

where i ranges from 1 to 240776. Not having to create a node in the graphical model to represent x1, x2, x3, x4, region.indexed, source.index, and person.indexed saves a large amount of space.

In the BUGS language, the type information is fine-grained: each component of a tensor can have different type information. This is quite distinct from the situation in STAN and can make it much easier to specify a statistical model. One common case is where some components of a tensor have been observed while other components need to be estimated. The STAN documentation suggests workarounds for these situations, but these are somewhat complex.

  • The type propagation is interesting and maybe useful. But we don’t necessarily need to implement a type system. A dirty way to get type information is simply do a dry run with some tricks.

Work flow

The statistical model and data are presented to the BUGS system in a series of stages. In the first stage the model text is parsed into a tree and the name table constructed. The data is then loaded and checked against the model. The data can be split over a number of source. Once all the data has been loaded the model is compiled. Compiling builds the graphical model and does a large number of checks on the consistency of the model. Finally initial values can be given or generated for the model.

The compiler creates a node in the graphical model for each scalar name and each component of a tensor name in the BUGS language model. The compiler checks that only one node is created for each scalar name or component of a tensor name.

Reading in a data source causes the compiler to create special nodes called constant nodes to hold the values of the data.

The compiler processes logical relations before stochastic relations. Any logical relations that only have constant nodes on their right hand side become new constant nodes with the appropriate fixed value. Even if a logical relation can not be reduced to a constant some parts of the relation might be reduced to constants.

Any constant nodes that have an associated stochastic relation become data nodes in the graphical model.

Logical relations in the BUGS Language

The OpenBUGS software compiles a description of a statistical model in the BUGS language into a graph of objects. Each relation in the statistical model gives rise to a node in the graph of objects. Each distinct type of relation in the statistical model is represented by a node of a distinct class. For stochastic relations there is a fixed set of distributions that can be used in the modelling. For logical relations the situation is more complex. The software can use arbitrary logical expressions build out of a fixed set of basic operators and functions. For each distinct logical expression a new software source code module is written to implement a class to represent that logical expression in the graph of objects. The software module is then compiled using the Components Pascal compiler and the executable code merged into the running OpenBUGS software using the run time loading linker.

The BUGS language description of a statistical model is parsed into a list of trees. The sub-trees that represent logical relations in the statistical model are first converted into a stack based representation and then into Component Pascal source code. The source code is generated in module BugsCPWrite and the source code is then compiled in module BugsCPCompiler. Usually the generated source code is not displayed. Checking the Verbose option in the Info menu will cause each each source code module generated by the OpenBUGS software to be displayed in a separate window.

One advantage of a stack based representation of an expression is that it is straight forward to use it to derive source code that calculates the derivative of the expression with respect to its arguments. This part of the source code generation is carried out in module BugsCPWrite in procedure WriteEvaluateDiffMethod. Each operator in the stack representation of the logical expression causes a snippet of Component Pascal code to be written. These code snippets are generally very simple with those of binary operators slightly more complex than those of unitary operators. Each binary operators can emit three different code snippets: the general case and two special snippets depending on whether the left or right operands are numerical constants. The only complex code snippet is when an operand that is a logical relation in the statistical model is pushed onto the stack – the > case of nested logical relations. In this case the nested logical relation will have its own code to calculate derivatives and these values can be passed up the nesting level.

The OpenBUGS software now uses a backward mode scheme to calculate the value of logical nodes in the statistical model. All the logical nodes in the statistical model are held in a global array and sorted according to their nesting level with unnested nodes at the start of the array. To evaluate all the logical nodes in the statistical model this array is then traversed and each logical node evaluated and the value stored in the node. The same scheme is used to calculate derivatives.

The graphs derived from the BUGS language representation of statistical models are generally sparse. The OpenBUGS software uses conditional independence arguments to exploit sparsity in the stochastic parts of the model. There is also a sparsity structure in logical relations.Each logical relation will often depend on just a few stochastic parents and derivatives with respect to other stochastic nodes in the model will be structurally zero. Each logical node has an associated array of stochastic parents for which the derivatives are non zero. Moving up the level of nesting the number of parents can grow. Dealing with this issue leads to the complexity in the code snippet for the operator that pushes a logical node onto the stack. These issues can be seen in the non-linear random effects model called Orange trees in volume II of the OpenBUGS examples. In this model eta[i,] is a function of phi[i,1], phi[i,2] and phi[i,3] where the phi are also logical functions of the stochastic theta[i,].

One refinement of the backward mode scheme used to calculate the value of logical nodes is to consider separately any logical nodes in the statistical model which are only used for prediction and do not affect the calculation of the joint probability distribution. These nodes need only be evaluated once per iteration of the inference algorithm. Examples of such nodes are sigma[k] and sigma.C in the Orange trees example. There is no need to evaluate the derivatives of these prediction nodes.

The workings of the backward mode scheme are easy to visualize when the inference algorithm updates all the stochastic nodes in the statistical model in one block. Local versions of the backward mode scheme can be used when the inference algorithm works on single nodes or when a small blocks of nodes are updated. Each stochastic node is given its own vector of logical nodes that depend on it either directly or via other logical nodes and this vector is sorted by nesting level. Each updater that works on small blocks of nodes contains a vector of logical nodes which is the union of the vectors of dependent logical nodes for each of its components.

The idea of the backward mode scheme for evaluating logical nodes can be used with caching in Metropolis Hastings sampling. First the vector of logical nodes depending on the relevant stochastic node(s) is evaluated and their values cached. The log of the conditional distribution is then calculated. Next a new value of the stochastic node is proposed. The vector of logical nodes is re-evaluated and the log of the > conditional distribution calculated. If the proposed value is rejected then the cache is used to set the vector of logical nodes back to its old values.

The OpenBUGS software also calculates what class of function each logical node is in terms of its stochastic parents. If the software can prove for example that a logical node is a linear function of its parents more efficient sampling algorithms can be used. If a linear relation can be proved then the calculation of derivatives can also be optimized in some cases because they will be constant and so only need to be calculated once. Generalized linear models are implemented in a way that allows fast calculation of derivatives. The structure of the algorithm to classify the functional form of logical nodes is very similar to that for derivatives and uses a backward mode scheme

BUGS separates management of logical and stochastic variables, essentially two graphs. Logical variables are stored in an array and values are updated with values in earlier positions of the array.

diff --git a/previews/PR150/R_interface/index.html b/previews/PR150/R_interface/index.html index 16d4da983..f5f5a40ac 100644 --- a/previews/PR150/R_interface/index.html +++ b/previews/PR150/R_interface/index.html @@ -123,4 +123,4 @@ 169.0 216.0 261.0 295.0 333.0 157.0 205.0 248.0 289.0 316.0 137.0 180.0 219.0 258.0 291.0 - 153.0 200.0 244.0 286.0 324.0

Please always verify the data before using.

+ 153.0 200.0 244.0 286.0 324.0

Please always verify the data before using.

diff --git a/previews/PR150/api/index.html b/previews/PR150/api/index.html index 9066d1615..9d2aff7bf 100644 --- a/previews/PR150/api/index.html +++ b/previews/PR150/api/index.html @@ -1,2 +1,2 @@ -General · JuliaBUGS.jl
Edit on GitHub

API

JuliaBUGS.Parser.@bugs — Macro
@bugs(prog::String, replace_period=true, no_enclosure=false)

Produce similar output as @bugs, but takes a string as input. This is useful for parsing original BUGS programs.

Arguments

  • prog::String: The BUGS program code as a string.
  • replace_period::Bool: If true, periods in the BUGS code will be replaced (default true).
  • no_enclosure::Bool: If true, the parser will not expect the program to be wrapped between model{ } (default false).
source
JuliaBUGS.compile — Function
compile(model_def[, data, initializations])

Compile a BUGS model into a log density problem.

Arguments

  • model_def::Expr: The BUGS model definition.
  • data::NamedTuple or AbstractDict: The data to be used in the model. If none is passed, the data will be assumed to be empty.
  • initializations::NamedTuple or AbstractDict: The initial values for the model parameters. If none is passed, the parameters will be assumed to be initialized to zero.
  • is_transformed::Bool=true: If true, the model parameters during inference will be transformed to the unconstrained space.

Returns

  • A BUGSModel object representing the compiled model.
source
JuliaBUGS.BUGSModel — Type
BUGSModel

The BUGSModel object is used for inference and represents the output of compilation. It fully implements the LogDensityProblems.jl interface.

Fields

  • transformed::Bool: Indicates whether the model parameters are in the transformed space.
  • untransformed_param_length::Int: The length of the parameters vector in the original space.
  • transformed_param_length::Int: The length of the parameters vector in the transformed space.
  • untransformed_var_lengths::Dict{VarName,Int}: A dictionary mapping the names of the variables to their lengths in the original space.
  • transformed_var_lengths::Dict{VarName,Int}: A dictionary mapping the names of the variables to their lengths in the transformed space.
  • varinfo::SimpleVarInfo: An instance of DynamicPPL.SimpleVarInfo, which is a dictionary-like data structure that maps both data and values of variables in the model to the corresponding values.
  • parameters::Vector{VarName}: A vector containing the names of the parameters in the model, defined as stochastic variables that are not observed. This vector should be consistent with sorted_nodes.
  • sorted_nodes::Vector{VarName}: A vector containing the names of all the variables in the model, sorted in topological order. In the case of a conditioned model, sorted_nodes include all the variables in parameters and the variables in the Markov blanket of parameters.
  • g::BUGSGraph: An instance of BUGSGraph, representing the dependency graph of the model.
  • base_model::Union{BUGSModel,Nothing}: If not Nothing, the model is a conditioned model; otherwise, it's the model returned by compile.
source
JuliaBUGS.ConcreteNodeInfo — Type
ConcreteNodeInfo

Defines the information stored in each node of the BUGS graph, encapsulating the essential characteristics and functions associated with a node within the BUGS model's dependency graph.

Fields

  • node_type::VariableTypes: Specifies whether the node is a stochastic or logical variable.
  • node_function_expr::Expr: The node function expression.
  • node_args::Vector{VarName}: A vector containing the names of the variables that are arguments to the node function.
source
+General · JuliaBUGS.jl
Edit on GitHub

API

JuliaBUGS.Parser.@bugs — Macro
@bugs(prog::String, replace_period=true, no_enclosure=false)

Produce similar output as @bugs, but takes a string as input. This is useful for parsing original BUGS programs.

Arguments

  • prog::String: The BUGS program code as a string.
  • replace_period::Bool: If true, periods in the BUGS code will be replaced (default true).
  • no_enclosure::Bool: If true, the parser will not expect the program to be wrapped between model{ } (default false).
source
JuliaBUGS.compile — Function
compile(model_def[, data, initializations])

Compile a BUGS model into a log density problem.

Arguments

  • model_def::Expr: The BUGS model definition.
  • data::NamedTuple or AbstractDict: The data to be used in the model. If none is passed, the data will be assumed to be empty.
  • initializations::NamedTuple or AbstractDict: The initial values for the model parameters. If none is passed, the parameters will be assumed to be initialized to zero.
  • is_transformed::Bool=true: If true, the model parameters during inference will be transformed to the unconstrained space.

Returns

  • A BUGSModel object representing the compiled model.
source
JuliaBUGS.BUGSModel — Type
BUGSModel

The BUGSModel object is used for inference and represents the output of compilation. It fully implements the LogDensityProblems.jl interface.

Fields

  • transformed::Bool: Indicates whether the model parameters are in the transformed space.
  • untransformed_param_length::Int: The length of the parameters vector in the original space.
  • transformed_param_length::Int: The length of the parameters vector in the transformed space.
  • untransformed_var_lengths::Dict{VarName,Int}: A dictionary mapping the names of the variables to their lengths in the original space.
  • transformed_var_lengths::Dict{VarName,Int}: A dictionary mapping the names of the variables to their lengths in the transformed space.
  • varinfo::SimpleVarInfo: An instance of DynamicPPL.SimpleVarInfo, which is a dictionary-like data structure that maps both data and values of variables in the model to the corresponding values.
  • parameters::Vector{VarName}: A vector containing the names of the parameters in the model, defined as stochastic variables that are not observed. This vector should be consistent with sorted_nodes.
  • sorted_nodes::Vector{VarName}: A vector containing the names of all the variables in the model, sorted in topological order. In the case of a conditioned model, sorted_nodes include all the variables in parameters and the variables in the Markov blanket of parameters.
  • g::BUGSGraph: An instance of BUGSGraph, representing the dependency graph of the model.
  • base_model::Union{BUGSModel,Nothing}: If not Nothing, the model is a conditioned model; otherwise, it's the model returned by compile.
source
JuliaBUGS.ConcreteNodeInfo — Type
ConcreteNodeInfo

Defines the information stored in each node of the BUGS graph, encapsulating the essential characteristics and functions associated with a node within the BUGS model's dependency graph.

Fields

  • node_type::VariableTypes: Specifies whether the node is a stochastic or logical variable.
  • node_function_expr::Expr: The node function expression.
  • node_args::Vector{VarName}: A vector containing the names of the variables that are arguments to the node function.
source
diff --git a/previews/PR150/bugs_examples/Leuk/index.html b/previews/PR150/bugs_examples/Leuk/index.html index c260ddffc..69acd9fa5 100644 --- a/previews/PR150/bugs_examples/Leuk/index.html +++ b/previews/PR150/bugs_examples/Leuk/index.html @@ -31,4 +31,4 @@ dL0.star[j] <- r * (t[j + 1] - t[j]) } beta ~ dnorm(0.0,0.000001) -}

Reference Results

VariableMeanMedianStandard DeviationMonte Carlo Error2.5% Value97.5% ValueStartSampleESS
S.placebo[1]0.92640.93740.049893.349E-40.80290.990910012000022184
S.placebo[17]0.044310.033440.039092.698E-40.0024780.148710012000020992
S.treat[1]0.98260.98630.014131.074E-40.94570.998210012000017315
S.treat[17]0.47670.47630.11980.0010090.24740.708610012000014104
beta1.5391.5240.42110.00340.74752.38810012000015340
+}

Reference Results

VariableMeanMedianStandard DeviationMonte Carlo Error2.5% Value97.5% ValueStartSampleESS
S.placebo[1]0.92640.93740.049893.349E-40.80290.990910012000022184
S.placebo[17]0.044310.033440.039092.698E-40.0024780.148710012000020992
S.treat[1]0.98260.98630.014131.074E-40.94570.998210012000017315
S.treat[17]0.47670.47630.11980.0010090.24740.708610012000014104
beta1.5391.5240.42110.00340.74752.38810012000015340
diff --git a/previews/PR150/bugs_lang/index.html b/previews/PR150/bugs_lang/index.html index 0a9c6d36b..c0a195835 100644 --- a/previews/PR150/bugs_lang/index.html +++ b/previews/PR150/bugs_lang/index.html @@ -25,4 +25,4 @@ x[a + 1] ~ dnorm(0, 1) } -list(y=c(1, 2, 3))

The reason for this is that the transformed variable and indices are evaluated simultaneously, and the transformed variable is not available at the time of indexing. This is a limitation of the current implementation.

+list(y=c(1, 2, 3))

The reason for this is that the transformed variable and indices are evaluated simultaneously, and the transformed variable is not available at the time of indexing. This is a limitation of the current implementation.

diff --git a/previews/PR150/differences/index.html b/previews/PR150/differences/index.html index 79754813e..25641bba7 100644 --- a/previews/PR150/differences/index.html +++ b/previews/PR150/differences/index.html @@ -1,2 +1,2 @@ -Differences Between BUGS and JuliaBUGS · JuliaBUGS.jl
Edit on GitHub

Differences Between BUGS and JuliaBUGS

Implicit Indexing

In BUGS, x[, ] is used for implicit indexing, which selects all elements from both the first and second dimensions. In JuliaBUGS, users must explicitly use Colon (:) like x[:, :] when using the @bugs macro. The @bugs macro will insert a Colon when given x[], however, the Julia parser will throw an error if given x[, ]. The original BUGS parser will automatically insert a Colon (:) when it encounters x[, ].

+Differences Between BUGS and JuliaBUGS · JuliaBUGS.jl
Edit on GitHub

Differences Between BUGS and JuliaBUGS

Implicit Indexing

In BUGS, x[, ] is used for implicit indexing, which selects all elements from both the first and second dimensions. In JuliaBUGS, users must explicitly use Colon (:) like x[:, :] when using the @bugs macro. The @bugs macro will insert a Colon when given x[], however, the Julia parser will throw an error if given x[, ]. The original BUGS parser will automatically insert a Colon (:) when it encounters x[, ].

diff --git a/previews/PR150/distributions/index.html b/previews/PR150/distributions/index.html index efab93447..834044dd3 100644 --- a/previews/PR150/distributions/index.html +++ b/previews/PR150/distributions/index.html @@ -1,3 +1,3 @@ -Distributions · JuliaBUGS.jl
Edit on GitHub
JuliaBUGS.BUGSPrimitives.dnorm — Function
dnorm(μ, τ)

Returns an instance of Normal with mean $μ$ and standard deviation $\frac{1}{√τ}$.

\[p(x|μ,τ) = \sqrt{\frac{τ}{2π}} e^{-τ \frac{(x-μ)^2}{2}}\]

source
JuliaBUGS.BUGSPrimitives.dlogis — Function
dlogis(μ, τ)

Return an instance of Logistic with location parameter $μ$ and scale parameter $\frac{1}{√τ}$.

\[p(x|μ,τ) = \frac{\sqrt{τ} e^{-\sqrt{τ}(x-μ)}}{(1+e^{-\sqrt{τ}(x-μ)})^2}\]

source
JuliaBUGS.BUGSPrimitives.dt — Function
dt(μ, τ, ν)

If $μ = 0$ and $σ = 1$, the function returns an instance of TDist with $ν$ degrees of freedom, location $μ$, and scale $σ = \frac{1}{\sqrt{τ}}$. Otherwise, it returns an instance of TDistShiftedScaled.

\[p(x|ν,μ,σ) = \frac{Γ((ν+1)/2)}{Γ(ν/2) \sqrt{νπσ}} -\left(1+\frac{1}{ν}\left(\frac{x-μ}{σ}\right)^2\right)^{-\frac{ν+1}{2}}\]

source
JuliaBUGS.BUGSPrimitives.TDistShiftedScaled — Type
TDistShiftedScaled(ν, μ, σ)

Student's t-distribution with $ν$ degrees of freedom, location $μ$, and scale $σ$.

This struct allows for a shift (determined by $μ$) and a scale (determined by $σ$) of the standard Student's t-distribution provided by the Distributions.jl package.

Only pdf and logpdf are implemented for this distribution.

See Also

TDist

source
JuliaBUGS.BUGSPrimitives.dflat — Function
dflat()

Returns an instance of Flat or TruncatedFlat if truncated.

Flat represents a flat (uniform) prior over the real line, which is an improper distribution. And TruncatedFlat represents a truncated version of the Flat distribution.

Only pdf, logpdf, minimum, and maximum are implemented for these Distributions.

When use in a model, the parameters always need to be initialized.

source
JuliaBUGS.BUGSPrimitives.dweib — Function
dweib(a, b)

Returns an instance of Weibull distribution object with shape parameter $a$ and scale parameter $\frac{1}{b}$.

The Weibull distribution is a common model for event times. The hazard or instantaneous risk of the event is $abx^{a-1}$. For $a < 1$ the hazard decreases with $x$; for $a > 1$ it increases. $a = 1$ results in the exponential distribution with constant hazard.

\[p(x|a,b) = abx^{a-1}e^{-b x^a}\]

source
JuliaBUGS.BUGSPrimitives.dgev — Function
dgev(μ, σ, η)

Returns an instance of GeneralizedExtremeValue with location $μ$, scale $σ$, and shape $η$.

\[p(x|μ,σ,η) = \frac{1}{σ} \left(1 + η \frac{x - μ}{σ}\right)^{-\frac{1}{η} - 1} e^{-\left(1 + η \frac{x - μ}{σ}\right)^{-\frac{1}{η}}}\]

where $\frac{η(x - μ)}{σ} > -1$.

source
JuliaBUGS.BUGSPrimitives.df — Function
df(n, m, μ=0, τ=1)

Returns an instance of F-distribution object with $n$ and $m$ degrees of freedom, location $μ$, and scale $τ$. This function is only valid when $μ = 0$ and $τ = 1$,

\[p(x|n, m, μ, τ) = \frac{\Gamma\left(\frac{n+m}{2}\right)}{\Gamma\left(\frac{n}{2}\right) \Gamma\left(\frac{m}{2}\right)} \left(\frac{n}{m}\right)^{\frac{n}{2}} \sqrt{τ} \left(\sqrt{τ}(x - μ)\right)^{\frac{n}{2}-1} \left(1 + \frac{n \sqrt{τ}(x-μ)}{m}\right)^{-\frac{n+m}{2}}\]

where $\frac{n \sqrt{τ} (x - μ)}{m} > -1$.

source
JuliaBUGS.BUGSPrimitives.dmnorm — Function
dmnorm(μ::AbstractVector, T::AbstractMatrix)

Returns an instance of Multivariate Normal with mean vector μ and covariance matrix $T^{-1}$.

\[p(x|μ,T) = (2π)^{-k/2} |T|^{1/2} e^{-1/2 (x-μ)' T (x-μ)}\]

where $k$ is the dimension of x.

source
JuliaBUGS.BUGSPrimitives.dmt — Function
dmt(μ::AbstractVector, T::AbstractMatrix, k)

Returns an instance of Multivariate T with mean vector $μ$, scale matrix $T^{-1}$, and $k$ degrees of freedom.

\[p(x|k,μ,Σ) = \frac{\Gamma((k+d)/2)}{\Gamma(k/2) (k\pi)^{p/2} |Σ|^{1/2}} \left(1 + \frac{1}{k} (x-μ)^T Σ^{-1} (x-μ)\right)^{-\frac{k+p}{2}}\]

where $p$ is the dimension of $x$.

source
JuliaBUGS.BUGSPrimitives.dwish — Function
dwish(R::AbstractMatrix, k)

Returns an instance of Wishart with $k$ degrees of freedom and the scale matrix $T^{-1}$.

\[p(X|R,k) = |R|^{k/2} |X|^{(k-p-1)/2} e^{-(1/2) tr(RX)} / (2^{kp/2} Γ_p(k/2))\]

where $p$ is the dimension of $X$, and it should be less than or equal to $k$.

source
JuliaBUGS.BUGSPrimitives.ddirich — Function
ddirich(θ::AbstractVector)

Return an instance of Dirichlet with parameters $θ_i$.

\[p(x|θ) = \frac{Γ(\sum θ)}{∏ Γ(θ)} ∏ x_i^{θ_i - 1}\]

where $\theta_i > 0, x_i \in [0, 1], \sum_i x_i = 1$

source
JuliaBUGS.BUGSPrimitives.dbin — Function
dbin(p, n)

Returns an instance of Binomial with number of trials n and success probability p.

\[p(x|n,p) = \binom{n}{x} p^x (1 - p)^{n-x}\]

end

where $\theta \in [0, 1], n \in \mathbb{Z}^+,$ and $x = 0, \ldots, n$.

source
JuliaBUGS.BUGSPrimitives.dhyper — Function
dhyper(n₁, n₂, m₁, ψ=1)

Returns an instance of Hypergeometric. This distribution is used when sampling without replacement from a population consisting of $n₁$ successes and $n₂$ failures, with $m₁$ being the number of trials or the sample size. The function currently only allows for $ψ = 1$.

\[p(x | n₁, n₂, m₁, \psi) = \frac{\binom{n₁}{x} \binom{n₂}{m₁ - x} \psi^x}{\sum_{i=u_0}^{u_1} \binom{n1}{i} \binom{n2}{m₁ - i} \psi^i}\]

where $u_0 = \max(0, m₁-n₂), u_1 = \min(n₁,m₁),$ and $u_0 \leq x \leq u_1$

source
+Distributions · JuliaBUGS.jl
Edit on GitHub
JuliaBUGS.BUGSPrimitives.dnorm — Function
dnorm(μ, τ)

Returns an instance of Normal with mean $μ$ and standard deviation $\frac{1}{√τ}$.

\[p(x|μ,τ) = \sqrt{\frac{τ}{2π}} e^{-τ \frac{(x-μ)^2}{2}}\]

source
JuliaBUGS.BUGSPrimitives.dlogis — Function
dlogis(μ, τ)

Return an instance of Logistic with location parameter $μ$ and scale parameter $\frac{1}{√τ}$.

\[p(x|μ,τ) = \frac{\sqrt{τ} e^{-\sqrt{τ}(x-μ)}}{(1+e^{-\sqrt{τ}(x-μ)})^2}\]

source
JuliaBUGS.BUGSPrimitives.dt — Function
dt(μ, τ, ν)

If $μ = 0$ and $σ = 1$, the function returns an instance of TDist with $ν$ degrees of freedom, location $μ$, and scale $σ = \frac{1}{\sqrt{τ}}$. Otherwise, it returns an instance of TDistShiftedScaled.

\[p(x|ν,μ,σ) = \frac{Γ((ν+1)/2)}{Γ(ν/2) \sqrt{νπσ}} +\left(1+\frac{1}{ν}\left(\frac{x-μ}{σ}\right)^2\right)^{-\frac{ν+1}{2}}\]

source
JuliaBUGS.BUGSPrimitives.TDistShiftedScaled — Type
TDistShiftedScaled(ν, μ, σ)

Student's t-distribution with $ν$ degrees of freedom, location $μ$, and scale $σ$.

This struct allows for a shift (determined by $μ$) and a scale (determined by $σ$) of the standard Student's t-distribution provided by the Distributions.jl package.

Only pdf and logpdf are implemented for this distribution.

See Also

TDist

source
JuliaBUGS.BUGSPrimitives.dflat — Function
dflat()

Returns an instance of Flat or TruncatedFlat if truncated.

Flat represents a flat (uniform) prior over the real line, which is an improper distribution. And TruncatedFlat represents a truncated version of the Flat distribution.

Only pdf, logpdf, minimum, and maximum are implemented for these Distributions.

When use in a model, the parameters always need to be initialized.

source
JuliaBUGS.BUGSPrimitives.dweib — Function
dweib(a, b)

Returns an instance of Weibull distribution object with shape parameter $a$ and scale parameter $\frac{1}{b}$.

The Weibull distribution is a common model for event times. The hazard or instantaneous risk of the event is $abx^{a-1}$. For $a < 1$ the hazard decreases with $x$; for $a > 1$ it increases. $a = 1$ results in the exponential distribution with constant hazard.

\[p(x|a,b) = abx^{a-1}e^{-b x^a}\]

source
JuliaBUGS.BUGSPrimitives.dgev — Function
dgev(μ, σ, η)

Returns an instance of GeneralizedExtremeValue with location $μ$, scale $σ$, and shape $η$.

\[p(x|μ,σ,η) = \frac{1}{σ} \left(1 + η \frac{x - μ}{σ}\right)^{-\frac{1}{η} - 1} e^{-\left(1 + η \frac{x - μ}{σ}\right)^{-\frac{1}{η}}}\]

where $\frac{η(x - μ)}{σ} > -1$.

source
JuliaBUGS.BUGSPrimitives.df — Function
df(n, m, μ=0, τ=1)

Returns an instance of F-distribution object with $n$ and $m$ degrees of freedom, location $μ$, and scale $τ$. This function is only valid when $μ = 0$ and $τ = 1$,

\[p(x|n, m, μ, τ) = \frac{\Gamma\left(\frac{n+m}{2}\right)}{\Gamma\left(\frac{n}{2}\right) \Gamma\left(\frac{m}{2}\right)} \left(\frac{n}{m}\right)^{\frac{n}{2}} \sqrt{τ} \left(\sqrt{τ}(x - μ)\right)^{\frac{n}{2}-1} \left(1 + \frac{n \sqrt{τ}(x-μ)}{m}\right)^{-\frac{n+m}{2}}\]

where $\frac{n \sqrt{τ} (x - μ)}{m} > -1$.

source
JuliaBUGS.BUGSPrimitives.dmnorm — Function
dmnorm(μ::AbstractVector, T::AbstractMatrix)

Returns an instance of Multivariate Normal with mean vector μ and covariance matrix $T^{-1}$.

\[p(x|μ,T) = (2π)^{-k/2} |T|^{1/2} e^{-1/2 (x-μ)' T (x-μ)}\]

where $k$ is the dimension of x.

source
JuliaBUGS.BUGSPrimitives.dmt — Function
dmt(μ::AbstractVector, T::AbstractMatrix, k)

Returns an instance of Multivariate T with mean vector $μ$, scale matrix $T^{-1}$, and $k$ degrees of freedom.

\[p(x|k,μ,Σ) = \frac{\Gamma((k+d)/2)}{\Gamma(k/2) (k\pi)^{p/2} |Σ|^{1/2}} \left(1 + \frac{1}{k} (x-μ)^T Σ^{-1} (x-μ)\right)^{-\frac{k+p}{2}}\]

where $p$ is the dimension of $x$.

source
JuliaBUGS.BUGSPrimitives.dwish — Function
dwish(R::AbstractMatrix, k)

Returns an instance of Wishart with $k$ degrees of freedom and the scale matrix $T^{-1}$.

\[p(X|R,k) = |R|^{k/2} |X|^{(k-p-1)/2} e^{-(1/2) tr(RX)} / (2^{kp/2} Γ_p(k/2))\]

where $p$ is the dimension of $X$, and it should be less than or equal to $k$.

source
JuliaBUGS.BUGSPrimitives.ddirich — Function
ddirich(θ::AbstractVector)

Return an instance of Dirichlet with parameters $θ_i$.

\[p(x|θ) = \frac{Γ(\sum θ)}{∏ Γ(θ)} ∏ x_i^{θ_i - 1}\]

where $\theta_i > 0, x_i \in [0, 1], \sum_i x_i = 1$

source
JuliaBUGS.BUGSPrimitives.dbin — Function
dbin(p, n)

Returns an instance of Binomial with number of trials n and success probability p.

\[p(x|n,p) = \binom{n}{x} p^x (1 - p)^{n-x}\]

end

where $\theta \in [0, 1], n \in \mathbb{Z}^+,$ and $x = 0, \ldots, n$.

source
JuliaBUGS.BUGSPrimitives.dhyper — Function
dhyper(n₁, n₂, m₁, ψ=1)

Returns an instance of Hypergeometric. This distribution is used when sampling without replacement from a population consisting of $n₁$ successes and $n₂$ failures, with $m₁$ being the number of trials or the sample size. The function currently only allows for $ψ = 1$.

\[p(x | n₁, n₂, m₁, \psi) = \frac{\binom{n₁}{x} \binom{n₂}{m₁ - x} \psi^x}{\sum_{i=u_0}^{u_1} \binom{n1}{i} \binom{n2}{m₁ - i} \psi^i}\]

where $u_0 = \max(0, m₁-n₂), u_1 = \min(n₁,m₁),$ and $u_0 \leq x \leq u_1$

source
diff --git a/previews/PR150/example/index.html b/previews/PR150/example/index.html index 037e4901d..5fb7a77df 100644 --- a/previews/PR150/example/index.html +++ b/previews/PR150/example/index.html @@ -134,4 +134,4 @@ init_params = [initial_θ for _ = 1:n_chains], # each chain has its own initial parameters discard_initial = n_adapts, progress = false, # Base.TTY creating problems in distributed setting -)

In this case, we pass two additional arguments to AbstractMCMC.sample:

Note that the init_params argument is now a vector of initial parameters for each chain. Sometimes the progress logger can cause problems in distributed setting, so we can disable it by setting progress = false.

More Examples

We have transcribed all the examples from the first volume of the BUGS Examples (original and transcribed). All programs and data are included, and can be compiled using the steps described in the tutorial above.

+)

In this case, we pass two additional arguments to AbstractMCMC.sample:

Note that the init_params argument is now a vector of initial parameters for each chain. Sometimes the progress logger can cause problems in distributed setting, so we can disable it by setting progress = false.

More Examples

We have transcribed all the examples from the first volume of the BUGS Examples (original and transcribed). All programs and data are included, and can be compiled using the steps described in the tutorial above.

diff --git a/previews/PR150/functions/index.html b/previews/PR150/functions/index.html index a380d1b3c..90a6c20d9 100644 --- a/previews/PR150/functions/index.html +++ b/previews/PR150/functions/index.html @@ -1,16 +1,5 @@ -Functions · JuliaBUGS.jl
Edit on GitHub

Most of the functions from BUGS have been implemented.

JuliaBUGS directly utilizes functions from the Julia Standard Library when they share the same names and functionalities. For functions not available in the Julia Standard Library and other popular libraries, we have developed equivalents within JuliaBUGS.BUGSPrimitives.

Function defined in Julia Standard Library

Warning

JuliaBUGS does not allow keyword arguments syntax

Base.abs — Function
abs(x)

The absolute value of x.

When abs is applied to signed integers, overflow may occur, resulting in the return of a negative value. This overflow occurs only when abs is applied to the minimum representable value of a signed integer. That is, when x == typemin(typeof(x)), abs(x) == x < 0, not -x as might be expected.

See also: abs2, unsigned, sign.

Examples

julia> abs(-3)
-3
-
-julia> abs(1 + im)
-1.4142135623730951
-
-julia> abs.(Int8[-128 -127 -126 0 126 127])  # overflow at typemin(Int8)
-1×6 Matrix{Int8}:
- -128  127  126  0  126  127
-
-julia> maximum(abs, [1, -2, 3, -4])
-4
source
Base.exp — Method
exp(x)

Compute the natural base exponential of x, in other words $ℯ^x$.

See also exp2, exp10 and cis.

Examples

julia> exp(1.0)
+Functions · JuliaBUGS.jl
Edit on GitHub
Most of the functions from BUGS have been implemented. 
JuliaBUGS directly utilizes functions from the Julia Standard Library when they share the same names and functionalities. For functions not available in the Julia Standard Library and other popular libraries, we have developed equivalents within JuliaBUGS.BUGSPrimitives.
Function defined in Julia Standard Library
Missing docstring.
Missing docstring for abs(x::Real). Check Documenter's build log for details.
Base.exp — Method
exp(x)
Compute the natural base exponential of x, in other words $ℯ^x$.
See also exp2, exp10 and cis.
Examples
julia> exp(1.0)
 2.718281828459045
 
 julia> exp(im * pi) ≈ cis(pi)
@@ -82,51 +71,8 @@
 julia> sort(A, dims = 2)
 2×2 Matrix{Int64}:
  3  4
- 1  2
source
Base.sin — Method
sin(x)
Compute sine of x, where x is in radians.
See also sind, sinpi, sincos, cis, asin.
Examples
julia> round.(sin.(range(0, 2pi, length=9)'), digits=3)
-1×9 Matrix{Float64}:
- 0.0  0.707  1.0  0.707  0.0  -0.707  -1.0  -0.707  -0.0
-
-julia> sind(45)
-0.7071067811865476
-
-julia> sinpi(1/4)
-0.7071067811865475
-
-julia> round.(sincos(pi/6), digits=3)
-(0.5, 0.866)
-
-julia> round(cis(pi/6), digits=3)
-0.866 + 0.5im
-
-julia> round(exp(im*pi/6), digits=3)
-0.866 + 0.5im
source
Base.tan — Method
tan(x)
Compute tangent of x, where x is in radians.
source
Base.asin — Method
asin(x)
Compute the inverse sine of x, where the output is in radians.
See also asind for output in degrees.
Examples
julia> asin.((0, 1/2, 1))
-(0.0, 0.5235987755982989, 1.5707963267948966)
-
-julia> asind.((0, 1/2, 1))
-(0.0, 30.000000000000004, 90.0)
source
Base.acos — Method
acos(x)
Compute the inverse cosine of x, where the output is in radians
source
Base.atan — Method
atan(y)
-atan(y, x)
Compute the inverse tangent of y or y/x, respectively.
For one argument, this is the angle in radians between the positive x-axis and the point (1, y), returning a value in the interval $[-\pi/2, \pi/2]$.
For two arguments, this is the angle in radians between the positive x-axis and the point (x, y), returning a value in the interval $[-\pi, \pi]$. This corresponds to a standard atan2 function. Note that by convention atan(0.0,x) is defined as $\pi$ and atan(-0.0,x) is defined as $-\pi$ when x < 0.
See also atand for degrees.
Examples
julia> rad2deg(atan(-1/√3))
--30.000000000000004
-
-julia> rad2deg(atan(-1, √3))
--30.000000000000004
-
-julia> rad2deg(atan(1, -√3))
-150.0
source
Statistics.mean — Method
mean(A::AbstractArray; dims)
Compute the mean of an array over the given dimensions.
Julia 1.1
mean for empty arrays requires at least Julia 1.1.
Examples
julia> using Statistics
-
-julia> A = [1 2; 3 4]
-2×2 Matrix{Int64}:
- 1  2
- 3  4
-
-julia> mean(A, dims=1)
-1×2 Matrix{Float64}:
- 2.0  3.0
-
-julia> mean(A, dims=2)
-2×1 Matrix{Float64}:
- 1.5
- 3.5
Function defined in LogExpFunctions
Missing docstring.
Missing docstring for !!! warning. Check Documenter's build log for details.
Function defined in LogExpFunctions
LogExpFunctions.logit — Function
logit(x)
 
The logit or log-odds transformation, defined as
\[\operatorname{logit}(x) = \log\left(\frac{x}{1-x}\right)\]
for $0 < x < 1$.
Its inverse is the logistic function.
LogExpFunctions.logistic — Function
logistic(x)
-
The logistic sigmoid function mapping a real number to a value in the interval $[0,1]$,
\[\sigma(x) = \frac{1}{e^{-x} + 1} = \frac{e^x}{1+e^x}.\]
Its inverse is the logit function.
Function defined in JuliaBUGS.BUGSPrimitives

+

The logistic sigmoid function mapping a real number to a value in the interval $[0,1]$,

\[\sigma(x) = \frac{1}{e^{-x} + 1} = \frac{e^x}{1+e^x}.\]

Its inverse is the logit function.

Function defined in JuliaBUGS.BUGSPrimitives

diff --git a/previews/PR150/graph_plotting/index.html b/previews/PR150/graph_plotting/index.html index 7871bb7d0..6d1f2a868 100644 --- a/previews/PR150/graph_plotting/index.html +++ b/previews/PR150/graph_plotting/index.html @@ -30,4 +30,4 @@ model = compile(model_def, NamedTuple(), inits)

TikzGraphs.jl.

using TikzGraphs
 TikzGraphs.plot(model)

TikzGraphs

GraphPlot.jl

using GraphPlot
 gplot(model)

GraphPlot

GraphMakie.jl

using GLMakie, GraphMakie
-graphplot(model)

GraphMakie

+graphplot(model)

GraphMakie

diff --git a/previews/PR150/index.html b/previews/PR150/index.html index 4ad256ac3..603fd65e7 100644 --- a/previews/PR150/index.html +++ b/previews/PR150/index.html @@ -1,2 +1,2 @@ -Introduction · JuliaBUGS.jl
Edit on GitHub

Introduction

JuliaBUGS is a graph-based probabilistic programming language and a component of the Turing ecosystem. The package aims to support modelling and inference for probabilistic programs written in the BUGS language.

This project is still in its early stage, with many key components needing to be completed.

Please refer to the example for usage information and a complete example.

What is BUGS?

The BUGS (Bayesian inference Using Gibbs Sampling) system is a probabilistic programming framework designed for specifying directed graphical models. Unlike certain other probabilistic programming languages (PPLs), such as Turing.jl or Pyro, the focus of BUGS is on specifying declarative relationships between nodes in a graph, which can be either logical or stochastic. This means that explicit declarations of variables, inputs, outputs, etc., are not required, and the order of statements is not critical.

The BUGS Approach and Benefits

Loops in BUGS are essentially a form of "plate notation," offering a concise way to express repetitive statements across many constant indices. Variables in BUGS are either the names of nodes within the program or constant parts of the "data" that must be combined with a model for instantiation.

A BUGS model provides a comprehensive representation of the relationships and dependencies among a set of variables within a Bayesian framework. Our goal is to support BUGS programs as much as possible while also incorporating Julia-specific syntax enhancements.

The key advantage of utilizing such a graph-based approach is the clarity it provides in understanding the dependencies and relationships within a complex system. These graphical models allow users to explicitly state the conditional dependencies between variables. This makes the model's structure and assumptions transparent, aiding both in the development and interpretation stages. Furthermore, using such a graphical approach makes it easier to apply advanced algorithms for model inference, as it enables more efficient computation by identifying and exploiting the structure of the model.

+Introduction · JuliaBUGS.jl
Edit on GitHub

Introduction

JuliaBUGS is a graph-based probabilistic programming language and a component of the Turing ecosystem. The package aims to support modelling and inference for probabilistic programs written in the BUGS language.

This project is still in its early stage, with many key components needing to be completed.

Please refer to the example for usage information and a complete example.

What is BUGS?

The BUGS (Bayesian inference Using Gibbs Sampling) system is a probabilistic programming framework designed for specifying directed graphical models. Unlike certain other probabilistic programming languages (PPLs), such as Turing.jl or Pyro, the focus of BUGS is on specifying declarative relationships between nodes in a graph, which can be either logical or stochastic. This means that explicit declarations of variables, inputs, outputs, etc., are not required, and the order of statements is not critical.

The BUGS Approach and Benefits

Loops in BUGS are essentially a form of "plate notation," offering a concise way to express repetitive statements across many constant indices. Variables in BUGS are either the names of nodes within the program or constant parts of the "data" that must be combined with a model for instantiation.

A BUGS model provides a comprehensive representation of the relationships and dependencies among a set of variables within a Bayesian framework. Our goal is to support BUGS programs as much as possible while also incorporating Julia-specific syntax enhancements.

The key advantage of utilizing such a graph-based approach is the clarity it provides in understanding the dependencies and relationships within a complex system. These graphical models allow users to explicitly state the conditional dependencies between variables. This makes the model's structure and assumptions transparent, aiding both in the development and interpretation stages. Furthermore, using such a graphical approach makes it easier to apply advanced algorithms for model inference, as it enables more efficient computation by identifying and exploiting the structure of the model.

diff --git a/previews/PR150/parser/index.html b/previews/PR150/parser/index.html index 91580b6d0..ae57b5c2c 100644 --- a/previews/PR150/parser/index.html +++ b/previews/PR150/parser/index.html @@ -1,2 +1,2 @@ -Parser · JuliaBUGS.jl
Edit on GitHub

BUGS Parser

The macro @bugs produces a Julia Expr object that represents the BUGS model definition.

If the input is a String, it's assumed to be a program in the original BUGS language. In this case, the macro will first convert the program to an equivalent Julia program, then use the Julia parser to parse the program into an Expr object.

Both model definitions written in Julia and those written in the original BUGS and subsequently parsed are now represented as a Julia Expr object. These objects go through syntax checking and post-processing to create the input for the compile function.

Below, we describe how the original BUGS program is translated to an equivalent Julia program and detail the post-processing done to the Expr object.

BUGS to Julia Translation

In this section, we refer to the translation program as the "parser" and the translating process as "parsing". Although the parser doesn't produce a syntax tree, it does follow the form of a recursive descent parser, building a Julia program in the form of a vector of tokens rather than a syntax tree.

This general implementation is heavily inspired by JuliaSyntax.jl, the official parser for Julia since version 1.10.

The BUGS parser implemented here takes a token stream with a recursive descent structure and checks the program's correctness. Here's how it works:

  1. Use tokenize to obtain the token vector.
  2. Inspect the tokens and build the Julia version of the program as a vector of tokens.
  3. Push the token to the Julia version of the program vector when appropriate.
  4. Detect errors and make necessary alterations to tokens, such as deletion, combination, or replacement.

During the recursive descent, BUGS syntax tokens will be translated into Julia syntax tokens. Some tokens will remain as they are, while others will be transformed, removed, or new tokens may be added.

The parser will throw an error if it encounters a program that does not adhere to strict BUGS syntax.

Some Notes on Error Recovery

The current error recovery is ad hoc and primarily rudimentary. If the program is correct, it will produce the correct result. If the program is syntactically or semantically incorrect, the token stream will not be pushed forward, resulting in failure.

The failure detection mechanism checks if two errors occur with the same "current token". If they do, the parser stops and reports the error. This ensures that the parser won't incorrectly parse a flawed program.

+Parser · JuliaBUGS.jl
Edit on GitHub

BUGS Parser

The macro @bugs produces a Julia Expr object that represents the BUGS model definition.

If the input is a String, it's assumed to be a program in the original BUGS language. In this case, the macro will first convert the program to an equivalent Julia program, then use the Julia parser to parse the program into an Expr object.

Both model definitions written in Julia and those written in the original BUGS and subsequently parsed are now represented as a Julia Expr object. These objects go through syntax checking and post-processing to create the input for the compile function.

Below, we describe how the original BUGS program is translated to an equivalent Julia program and detail the post-processing done to the Expr object.

BUGS to Julia Translation

In this section, we refer to the translation program as the "parser" and the translating process as "parsing". Although the parser doesn't produce a syntax tree, it does follow the form of a recursive descent parser, building a Julia program in the form of a vector of tokens rather than a syntax tree.

This general implementation is heavily inspired by JuliaSyntax.jl, the official parser for Julia since version 1.10.

The BUGS parser implemented here takes a token stream with a recursive descent structure and checks the program's correctness. Here's how it works:

  1. Use tokenize to obtain the token vector.
  2. Inspect the tokens and build the Julia version of the program as a vector of tokens.
  3. Push the token to the Julia version of the program vector when appropriate.
  4. Detect errors and make necessary alterations to tokens, such as deletion, combination, or replacement.

During the recursive descent, BUGS syntax tokens will be translated into Julia syntax tokens. Some tokens will remain as they are, while others will be transformed, removed, or new tokens may be added.

The parser will throw an error if it encounters a program that does not adhere to strict BUGS syntax.

Some Notes on Error Recovery

The current error recovery is ad hoc and primarily rudimentary. If the program is correct, it will produce the correct result. If the program is syntactically or semantically incorrect, the token stream will not be pushed forward, resulting in failure.

The failure detection mechanism checks if two errors occur with the same "current token". If they do, the parser stops and reports the error. This ensures that the parser won't incorrectly parse a flawed program.

diff --git a/previews/PR150/pitfalls/index.html b/previews/PR150/pitfalls/index.html index 9108dfbdf..00d1590c0 100644 --- a/previews/PR150/pitfalls/index.html +++ b/previews/PR150/pitfalls/index.html @@ -12,4 +12,4 @@ } p[1] <- 0.5 p[2] <- 0.5 -}

For a variable to be used as an observation in loop bounds or indexing, it must be part of the provided data, not a transformed variable.

This behavior is maintained in the current version of JuliaBUGS, although it was prohibited in the earlier SymbolicPPL.

Possible Check Implementation in JuliaBUGS

Implementing a check for this behavior in JuliaBUGS is feasible. A simplistic approach could be to invalidate (e.g., mark as missing) all observations after the first pass and verify if any are used in loop bounds or indexing. However, there is currently no plan to implement this check.

+}

For a variable to be used as an observation in loop bounds or indexing, it must be part of the provided data, not a transformed variable.

This behavior is maintained in the current version of JuliaBUGS, although it was prohibited in the earlier SymbolicPPL.

Possible Check Implementation in JuliaBUGS

Implementing a check for this behavior in JuliaBUGS is feasible. A simplistic approach could be to invalidate (e.g., mark as missing) all observations after the first pass and verify if any are used in loop bounds or indexing. However, there is currently no plan to implement this check.

diff --git a/previews/PR150/search/index.html b/previews/PR150/search/index.html index 0673b7524..ab03b4f4d 100644 --- a/previews/PR150/search/index.html +++ b/previews/PR150/search/index.html @@ -1,2 +1,2 @@ -Search · JuliaBUGS.jl

Loading search...

    +Search · JuliaBUGS.jl

    Loading search...

      diff --git a/previews/PR150/search_index.js b/previews/PR150/search_index.js index 5f0619d15..58134077c 100644 --- a/previews/PR150/search_index.js +++ b/previews/PR150/search_index.js @@ -1,3 +1,3 @@ var documenterSearchIndex = {"docs": -[{"location":"pitfalls/#Understanding-Pitfalls-in-Model-Definitions","page":"Pitfalls","title":"Understanding Pitfalls in Model Definitions","text":"","category":"section"},{"location":"pitfalls/#Consequence-of-Observations-on-Model-Parameters","page":"Pitfalls","title":"Consequence of Observations on Model Parameters","text":"","category":"section"},{"location":"pitfalls/","page":"Pitfalls","title":"Pitfalls","text":"When providing observations for the parameters of a model, the dependencies may become disrupted. Consider the following example written in Julia:","category":"page"},{"location":"pitfalls/","page":"Pitfalls","title":"Pitfalls","text":"model_def = @bugs begin\n a ~ Normal(0, 1)\n b ~ Normal(0, 1)\n c ~ Normal(a, b)\nend\n\ndata = (a=1.0, b=2.0)","category":"page"},{"location":"pitfalls/","page":"Pitfalls","title":"Pitfalls","text":"In this scenario, the generated graph will lack the edges a -> c and b -> c, leading the node function of c to become c ~ Normal(1.0, 2.0).","category":"page"},{"location":"pitfalls/#Ambiguity-Between-Data-and-Observed-Stochastic-Variables","page":"Pitfalls","title":"Ambiguity Between Data and Observed Stochastic Variables","text":"","category":"section"},{"location":"pitfalls/","page":"Pitfalls","title":"Pitfalls","text":"A subtle and possibly contentious feature of BUGS syntax is that the observation value of a stochastic variable is treated identically to any model parameters supplied in the data. Here's a legal example in BUGS:","category":"page"},{"location":"pitfalls/","page":"Pitfalls","title":"Pitfalls","text":"model {\n N ~ dcat(p[])\n for (i in 1:N) {\n y[i] ~ dnorm(mu, tau)\n }\n p[1] <- 0.5\n p[2] <- 0.5\n}","category":"page"},{"location":"pitfalls/","page":"Pitfalls","title":"Pitfalls","text":"For a variable to be used as an observation in loop bounds or indexing, it must be part of the provided data, not a transformed variable.","category":"page"},{"location":"pitfalls/","page":"Pitfalls","title":"Pitfalls","text":"This behavior is maintained in the current version of JuliaBUGS, although it was prohibited in the earlier SymbolicPPL.","category":"page"},{"location":"pitfalls/#Possible-Check-Implementation-in-JuliaBUGS","page":"Pitfalls","title":"Possible Check Implementation in JuliaBUGS","text":"","category":"section"},{"location":"pitfalls/","page":"Pitfalls","title":"Pitfalls","text":"Implementing a check for this behavior in JuliaBUGS is feasible. A simplistic approach could be to invalidate (e.g., mark as missing) all observations after the first pass and verify if any are used in loop bounds or indexing. However, there is currently no plan to implement this check.","category":"page"},{"location":"distributions/","page":"Distributions","title":"Distributions","text":"dnorm\ndlogis\ndt\nTDistShiftedScaled\nddexp\ndflat\nFlat\nTruncatedFlat\ndexp\ndchisqr\ndweib\ndlnorm\ndgamma\ndpar\ndgev\ndgpar\ndf\ndunif\ndbeta\ndmnorm\ndmt\ndwish\nddirich\ndbern\ndbin\ndcat\ndpois\ndgeom\ndnegbin\ndbetabin\ndhyper\ndmulti","category":"page"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dnorm","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dnorm","text":"dnorm(μ, τ)\n\nReturns an instance of Normal with mean μ and standard deviation frac1τ. \n\np(xμτ) = sqrtfracτ2π e^-τ frac(x-μ)^22\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dlogis","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dlogis","text":"dlogis(μ, τ)\n\nReturn an instance of Logistic with location parameter μ and scale parameter frac1τ.\n\np(xμτ) = fracsqrtτ e^-sqrtτ(x-μ)(1+e^-sqrtτ(x-μ))^2\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dt","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dt","text":"dt(μ, τ, ν)\n\nIf μ = 0 and σ = 1, the function returns an instance of TDist with ν degrees of freedom, location μ, and scale σ = frac1sqrtτ. Otherwise, it returns an instance of TDistShiftedScaled.\n\np(xνμσ) = fracΓ((ν+1)2)Γ(ν2) sqrtνπσ\nleft(1+frac1νleft(fracx-μσright)^2right)^-fracν+12\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.TDistShiftedScaled","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.TDistShiftedScaled","text":"TDistShiftedScaled(ν, μ, σ)\n\nStudent's t-distribution with ν degrees of freedom, location μ, and scale σ. \n\nThis struct allows for a shift (determined by μ) and a scale (determined by σ) of the standard Student's t-distribution provided by the Distributions.jl package. \n\nOnly pdf and logpdf are implemented for this distribution.\n\nSee Also\n\nTDist\n\n\n\n\n\n","category":"type"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.ddexp","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.ddexp","text":"ddexp(μ, τ)\n\nReturn an instance of Laplace (Double Exponential) with location μ and scale frac1sqrtτ.\n\np(xμτ) = fracsqrtτ2 e^-sqrtτ x-μ\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dflat","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dflat","text":"dflat()\n\nReturns an instance of Flat or TruncatedFlat if truncated.\n\nFlat represents a flat (uniform) prior over the real line, which is an improper distribution. And TruncatedFlat represents a truncated version of the Flat distribution.\n\nOnly pdf, logpdf, minimum, and maximum are implemented for these Distributions.\n\nWhen use in a model, the parameters always need to be initialized.\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.Flat","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.Flat","text":"Flat\n\nThe flat distribution mimicking the behavior of the dflat distribution in the BUGS family of softwares.\n\n\n\n\n\n","category":"type"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.TruncatedFlat","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.TruncatedFlat","text":"TruncatedFlat\n\nTruncated version of the Flat distribution.\n\n\n\n\n\n","category":"type"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dexp","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dexp","text":"dexp(λ)\n\nReturns an instance of Exponential with rate frac1λ.\n\np(xλ) = λ e^-λ x\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dchisqr","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dchisqr","text":"dchisqr(k)\n\nReturns an instance of Chi-squared with k degrees of freedom.\n\np(xk) = frac12^k2 Γ(k2) x^k2 - 1 e^-x2\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dweib","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dweib","text":"dweib(a, b)\n\nReturns an instance of Weibull distribution object with shape parameter a and scale parameter frac1b.\n\nThe Weibull distribution is a common model for event times. The hazard or instantaneous risk of the event is abx^a-1. For a 1 the hazard decreases with x; for a 1 it increases. a = 1 results in the exponential distribution with constant hazard.\n\np(xab) = abx^a-1e^-b x^a\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dlnorm","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dlnorm","text":"dlnorm(μ, τ)\n\nReturns an instance of LogNormal with location μ and scale frac1sqrtτ.\n\np(xμτ) = fracsqrtτxsqrt2π e^-τ2 (log(x) - μ)^2\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dgamma","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dgamma","text":"dgamma(a, b)\n\nReturns an instance of Gamma with shape a and scale frac1b.\n\np(xab) = fracb^aΓ(a) x^a-1 e^-bx\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dpar","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dpar","text":"dpar(a, b)\n\nReturns an instance of Pareto with scale parameter b and shape parameter a.\n\np(xab) = fraca b^ax^a+1\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dgev","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dgev","text":"dgev(μ, σ, η)\n\nReturns an instance of GeneralizedExtremeValue with location μ, scale σ, and shape η.\n\np(xμση) = frac1σ left(1 + η fracx - μσright)^-frac1η - 1 e^-left(1 + η fracx - μσright)^-frac1η\n\nwhere fracη(x - μ)σ -1.\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dgpar","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dgpar","text":"dgpar(μ, σ, η)\n\nReturns an instance of GeneralizedPareto with location μ, scale σ, and shape η.\n\np(xμση) = frac1σ (1 + η ((x - μ)σ))^-1η - 1\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.df","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.df","text":"df(n, m, μ=0, τ=1)\n\nReturns an instance of F-distribution object with n and m degrees of freedom, location μ, and scale τ. This function is only valid when μ = 0 and τ = 1,\n\np(xn m μ τ) = fracGammaleft(fracn+m2right)Gammaleft(fracn2right) Gammaleft(fracm2right) left(fracnmright)^fracn2 sqrtτ left(sqrtτ(x - μ)right)^fracn2-1 left(1 + fracn sqrtτ(x-μ)mright)^-fracn+m2\n\nwhere fracn sqrtτ (x - μ)m -1.\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dunif","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dunif","text":"dunif(a, b)\n\nReturns an instance of Uniform with lower bound a and upper bound b.\n\np(xab) = frac1b - a\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dbeta","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dbeta","text":"dbeta(a, b)\n\nReturns an instance of Beta with shape parameters a and b.\n\np(xab) = fracGamma(a + b)Gamma(a)Gamma(b) x^a-1 (1 - x)^b-1\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dmnorm","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dmnorm","text":"dmnorm(μ::AbstractVector, T::AbstractMatrix)\n\nReturns an instance of Multivariate Normal with mean vector μ and covariance matrix T^-1.\n\np(xμT) = (2π)^-k2 T^12 e^-12 (x-μ) T (x-μ)\n\nwhere k is the dimension of x.\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dmt","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dmt","text":"dmt(μ::AbstractVector, T::AbstractMatrix, k)\n\nReturns an instance of Multivariate T with mean vector μ, scale matrix T^-1, and k degrees of freedom.\n\np(xkμΣ) = fracGamma((k+d)2)Gamma(k2) (kpi)^p2 Σ^12 left(1 + frac1k (x-μ)^T Σ^-1 (x-μ)right)^-frack+p2\n\nwhere p is the dimension of x.\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dwish","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dwish","text":"dwish(R::AbstractMatrix, k)\n\nReturns an instance of Wishart with k degrees of freedom and the scale matrix T^-1.\n\np(XRk) = R^k2 X^(k-p-1)2 e^-(12) tr(RX) (2^kp2 Γ_p(k2))\n\nwhere p is the dimension of X, and it should be less than or equal to k. \n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.ddirich","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.ddirich","text":"ddirich(θ::AbstractVector)\n\nReturn an instance of Dirichlet with parameters θ_i.\n\np(xθ) = fracΓ(sum θ) Γ(θ) x_i^θ_i - 1\n\nwhere theta_i 0 x_i in 0 1 sum_i x_i = 1\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dbern","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dbern","text":"dbern(p)\n\nReturn an instance of Bernoulli with success probability p.\n\np(xp) = p^x (1 - p)^1-x\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dbin","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dbin","text":"dbin(p, n)\n\nReturns an instance of Binomial with number of trials n and success probability p.\n\np(xnp) = binomnx p^x (1 - p)^n-x\n\nend\n\nwhere theta in 0 1 n in mathbbZ^+ and x = 0 ldots n.\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dcat","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dcat","text":"dcat(p)\n\nReturns an instance of Categorical with probabilities p.\n\np(xp) = px\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dpois","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dpois","text":"dpois(θ)\n\nReturns an instance of Poisson with mean (and variance) θ.\n\np(xθ) = e^-θ θ^x x\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dgeom","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dgeom","text":"dgeom(θ)\n\nReturns an instance of Geometric with success probability θ.\n\np(xθ) = (1 - θ)^x-1 θ\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dnegbin","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dnegbin","text":"dnegbin(p, r)\n\nReturns an instance of Negative Binomial with number of failures r and success probability p.\n\nP(xrp) = binomx + r - 1x (1 - p)^x p^r\n\nwhere x in mathbbZ^+.\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dbetabin","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dbetabin","text":"dbetabin(a, b, n)\n\nReturns an instance of Beta Binomial with number of trials n and shape parameters a and b.\n\nP(xa b n) = fracbinomnx binoma + b - 1a + x - 1binoma + b + n - 1n\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dhyper","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dhyper","text":"dhyper(n₁, n₂, m₁, ψ=1)\n\nReturns an instance of Hypergeometric. This distribution is used when sampling without replacement from a population consisting of n₁ successes and n₂ failures, with m₁ being the number of trials or the sample size. The function currently only allows for ψ = 1.\n\np(x n₁ n₂ m₁ psi) = fracbinomn₁x binomn₂m₁ - x psi^xsum_i=u_0^u_1 binomn1i binomn2m₁ - i psi^i\n\nwhere u_0 = max(0 m₁-n₂) u_1 = min(n₁m₁) and u_0 leq x leq u_1\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dmulti","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dmulti","text":"dmulti(θ::AbstractVector, n)\n\nReturns an instance Multinomial with number of trials n and success probabilities θ.\n\nP(xnθ) = fracn_r x_r _r θ_r^x_r\n\n\n\n\n\n","category":"function"},{"location":"BUGS_notes/#Miscellaneous-Notes-on-BUGS","page":"Notes on BUGS Implementations","title":"Miscellaneous Notes on BUGS","text":"","category":"section"},{"location":"BUGS_notes/","page":"Notes on BUGS Implementations","title":"Notes on BUGS Implementations","text":"Here are some exert from BUGS Developer Manual and notes on the original BUGS implementations. ","category":"page"},{"location":"BUGS_notes/#Lexing","page":"Notes on BUGS Implementations","title":"Lexing","text":"","category":"section"},{"location":"BUGS_notes/","page":"Notes on BUGS Implementations","title":"Notes on BUGS Implementations","text":"The BUGS language has the convention that if a name is followed immediately by a round bracket, that is by a \"(\", then the names is a reserved name in the BUGS language and does not represent a variable in the model.\nBy scanning the stream of tokens that constitute a BUGS language model the names of all the variables in the model can be found.","category":"page"},{"location":"BUGS_notes/#Table-of-Names","page":"Notes on BUGS Implementations","title":"Table of Names","text":"","category":"section"},{"location":"BUGS_notes/","page":"Notes on BUGS Implementations","title":"Notes on BUGS Implementations","text":"The BUGS language compiler expands all the for loops in the model and records the value of the indices of each use of a tensor on the left hand side of each relation.\nThe range of each index, for a tensor, is set at the maximum value observed value of the index and added to the name table. There is one exception to this procedure for finding index bounds: names that are data, that is in the data source, have the ranges of their indices fixed in the data source.\nEach scalar and each component of a tensor used on the right hand side of a relation must occur either on the left hand side of a relation and or in a data source.","category":"page"},{"location":"BUGS_notes/#Data-Transformations","page":"Notes on BUGS Implementations","title":"Data Transformations","text":"","category":"section"},{"location":"BUGS_notes/","page":"Notes on BUGS Implementations","title":"Notes on BUGS Implementations","text":"If the compiler can prove that a logical assignment can be evaluated to a constant then the assignment is called a data transformation. This occurs if an assignment's right hand side does not depend on any variable quantities. The BUGS language has a general rule that there must only be one assignment statement for each scalar or component of a tensor. This rule is slightly relaxed for data transformations. The language allows a logical assignment and a stochastic assignment to the same scalar or tensor component if and only if the logical assignment is a data transformation. ","category":"page"},{"location":"BUGS_notes/#Generated-Quantities","page":"Notes on BUGS Implementations","title":"Generated Quantities","text":"","category":"section"},{"location":"BUGS_notes/","page":"Notes on BUGS Implementations","title":"Notes on BUGS Implementations","text":"Only need to be evaluated after the inference algorithm has finished its task. ","category":"page"},{"location":"BUGS_notes/","page":"Notes on BUGS Implementations","title":"Notes on BUGS Implementations","text":"Generally, these are leaf nodes that logical variables\nIn the case of stochastic variables that are leaf nodes, do “forward sampling”, also part of the generated Quantities","category":"page"},{"location":"BUGS_notes/#Computation","page":"Notes on BUGS Implementations","title":"Computation","text":"","category":"section"},{"location":"BUGS_notes/","page":"Notes on BUGS Implementations","title":"Notes on BUGS Implementations","text":"All the nodes in the graphical model representing logical relations are placed into an","category":"page"},{"location":"BUGS_notes/","page":"Notes on BUGS Implementations","title":"Notes on BUGS Implementations","text":"array and sorted by their nesting level with the first array entries only depending on quantities defined by stochastic relations. Traversing this array and evaluating nodes gives up to date values to all logical relations.","category":"page"},{"location":"BUGS_notes/#Types","page":"Notes on BUGS Implementations","title":"Types","text":"","category":"section"},{"location":"BUGS_notes/","page":"Notes on BUGS Implementations","title":"Notes on BUGS Implementations","text":"The BUGS compiler uses the properties of the distribution on the right-hand side of a stochastic assignment statement to make deductions about the variable on the left-hand side. For example, r ~ dbin(p, n) implies that r is integer-valued, while x ~ dnorm(mu, tau) implies that x is real-valued.Some distributions are real-valued but have support on a restricted range of the reals. For example, p ~ dbeta(a, b) implies that p is real-valued with support on the unit interval, while x ~ dgamma(r, lambda) implies that x is real-valued but with support on the positive real line.There are two multivariate distributions in the BUGS language, the Dirichlet and the Wishart, that have support on a complex subspace of the reals. The Dirichlet has support on the unit simplex, while the Wishart has support on symmetric positive definite matrices.The BUGS compiler tries to infer if logical relations return an integer value by looking at whether their parents are integer-valued and the operators that combine the values of their parents into the return value. For example, in the cure model example above, the logical relation state1[i] <- state[i] + 1 is integer-valued because state[i] is a Bernoulli variable and therefore integer, the literal 1 is integer, and the sum of two integers is an integer.When the BUGS system reads in data from a data source, it can tag whether the number read is an integer or a real and propagate this information to logical relations. Again, using the cure model as an example, the statement t[i] <- x[i] + y[i] is integer-valued because both x and y are data and are given as integers in the data source.One special type of data is constants: that is just numbers with no associated distribution. Constants have many uses in BUGS language models, but one of the most important is as covariates. A model can contain a large number of constants that are used as covariates. Because of the possible large numbers of these covariate-type constants, they are given special treatment by the BUGS compiler. If a name read in from a data source is only used on the right-hand side of logical relations, no nodes in the graphical model are created to hold its values; they are directly incorporated in the objects that represent the right-hand sides of the logical relations.For example, the large Methadone model contains the regression:mu.indexed[i] <- beta[1] * x1[i] + beta[2] * x2[i] + beta[3] * x3[i] + beta[4] * x4[i] + region.effect[region.indexed[i]] + source.effect[region.indexed[i]] * source.indexed[i] + person.effect[person.indexed[i]]where i ranges from 1 to 240776. Not having to create a node in the graphical model to represent x1, x2, x3, x4, region.indexed, source.index, and person.indexed saves a large amount of space.In the BUGS language, the type information is fine-grained: each component of a tensor can have different type information. This is quite distinct from the situation in STAN and can make it much easier to specify a statistical model. One common case is where some components of a tensor have been observed while other components need to be estimated. The STAN documentation suggests workarounds for these situations, but these are somewhat complex.","category":"page"},{"location":"BUGS_notes/","page":"Notes on BUGS Implementations","title":"Notes on BUGS Implementations","text":"The type propagation is interesting and maybe useful. But we don’t necessarily need to implement a type system. A dirty way to get type information is simply do a dry run with some tricks.","category":"page"},{"location":"BUGS_notes/#Work-flow","page":"Notes on BUGS Implementations","title":"Work flow","text":"","category":"section"},{"location":"BUGS_notes/","page":"Notes on BUGS Implementations","title":"Notes on BUGS Implementations","text":"The statistical model and data are presented to the BUGS system in a series of stages. In the first stage the model text is parsed into a tree and the name table constructed. The data is then loaded and checked against the model. The data can be split over a number of source. Once all the data has been loaded the model is compiled. Compiling builds the graphical model and does a large number of checks on the consistency of the model. Finally initial values can be given or generated for the model.","category":"page"},{"location":"BUGS_notes/","page":"Notes on BUGS Implementations","title":"Notes on BUGS Implementations","text":"The compiler creates a node in the graphical model for each scalar name and each component of a tensor name in the BUGS language model. The compiler checks that only one node is created for each scalar name or component of a tensor name.","category":"page"},{"location":"BUGS_notes/","page":"Notes on BUGS Implementations","title":"Notes on BUGS Implementations","text":"Reading in a data source causes the compiler to create special nodes called constant nodes to hold the values of the data.","category":"page"},{"location":"BUGS_notes/","page":"Notes on BUGS Implementations","title":"Notes on BUGS Implementations","text":"The compiler processes logical relations before stochastic relations. Any logical relations that only have constant nodes on their right hand side become new constant nodes with the appropriate fixed value. Even if a logical relation can not be reduced to a constant some parts of the relation might be reduced to constants.","category":"page"},{"location":"BUGS_notes/","page":"Notes on BUGS Implementations","title":"Notes on BUGS Implementations","text":"Any constant nodes that have an associated stochastic relation become data nodes in the graphical model.","category":"page"},{"location":"BUGS_notes/#Logical-relations-in-the-BUGS-Language","page":"Notes on BUGS Implementations","title":"Logical relations in the BUGS Language","text":"","category":"section"},{"location":"BUGS_notes/","page":"Notes on BUGS Implementations","title":"Notes on BUGS Implementations","text":"The OpenBUGS software compiles a description of a statistical model in the BUGS language into a graph of objects. Each relation in the statistical model gives rise to a node in the graph of objects. Each distinct type of relation in the statistical model is represented by a node of a distinct class. For stochastic relations there is a fixed set of distributions that can be used in the modelling. For logical relations the situation is more complex. The software can use arbitrary logical expressions build out of a fixed set of basic operators and functions. For each distinct logical expression a new software source code module is written to implement a class to represent that logical expression in the graph of objects. The software module is then compiled using the Components Pascal compiler and the executable code merged into the running OpenBUGS software using the run time loading linker.The BUGS language description of a statistical model is parsed into a list of trees. The sub-trees that represent logical relations in the statistical model are first converted into a stack based representation and then into Component Pascal source code. The source code is generated in module BugsCPWrite and the source code is then compiled in module BugsCPCompiler. Usually the generated source code is not displayed. Checking the Verbose option in the Info menu will cause each each source code module generated by the OpenBUGS software to be displayed in a separate window.One advantage of a stack based representation of an expression is that it is straight forward to use it to derive source code that calculates the derivative of the expression with respect to its arguments. This part of the source code generation is carried out in module BugsCPWrite in procedure WriteEvaluateDiffMethod. Each operator in the stack representation of the logical expression causes a snippet of Component Pascal code to be written. These code snippets are generally very simple with those of binary operators slightly more complex than those of unitary operators. Each binary operators can emit three different code snippets: the general case and two special snippets depending on whether the left or right operands are numerical constants. The only complex code snippet is when an operand that is a logical relation in the statistical model is pushed onto the stack – the > case of nested logical relations. In this case the nested logical relation will have its own code to calculate derivatives and these values can be passed up the nesting level.The OpenBUGS software now uses a backward mode scheme to calculate the value of logical nodes in the statistical model. All the logical nodes in the statistical model are held in a global array and sorted according to their nesting level with unnested nodes at the start of the array. To evaluate all the logical nodes in the statistical model this array is then traversed and each logical node evaluated and the value stored in the node. The same scheme is used to calculate derivatives.The graphs derived from the BUGS language representation of statistical models are generally sparse. The OpenBUGS software uses conditional independence arguments to exploit sparsity in the stochastic parts of the model. There is also a sparsity structure in logical relations.Each logical relation will often depend on just a few stochastic parents and derivatives with respect to other stochastic nodes in the model will be structurally zero. Each logical node has an associated array of stochastic parents for which the derivatives are non zero. Moving up the level of nesting the number of parents can grow. Dealing with this issue leads to the complexity in the code snippet for the operator that pushes a logical node onto the stack. These issues can be seen in the non-linear random effects model called Orange trees in volume II of the OpenBUGS examples. In this model eta[i,] is a function of phi[i,1], phi[i,2] and phi[i,3] where the phi are also logical functions of the stochastic theta[i,].One refinement of the backward mode scheme used to calculate the value of logical nodes is to consider separately any logical nodes in the statistical model which are only used for prediction and do not affect the calculation of the joint probability distribution. These nodes need only be evaluated once per iteration of the inference algorithm. Examples of such nodes are sigma[k] and sigma.C in the Orange trees example. There is no need to evaluate the derivatives of these prediction nodes.The workings of the backward mode scheme are easy to visualize when the inference algorithm updates all the stochastic nodes in the statistical model in one block. Local versions of the backward mode scheme can be used when the inference algorithm works on single nodes or when a small blocks of nodes are updated. Each stochastic node is given its own vector of logical nodes that depend on it either directly or via other logical nodes and this vector is sorted by nesting level. Each updater that works on small blocks of nodes contains a vector of logical nodes which is the union of the vectors of dependent logical nodes for each of its components.The idea of the backward mode scheme for evaluating logical nodes can be used with caching in Metropolis Hastings sampling. First the vector of logical nodes depending on the relevant stochastic node(s) is evaluated and their values cached. The log of the conditional distribution is then calculated. Next a new value of the stochastic node is proposed. The vector of logical nodes is re-evaluated and the log of the > conditional distribution calculated. If the proposed value is rejected then the cache is used to set the vector of logical nodes back to its old values.The OpenBUGS software also calculates what class of function each logical node is in terms of its stochastic parents. If the software can prove for example that a logical node is a linear function of its parents more efficient sampling algorithms can be used. If a linear relation can be proved then the calculation of derivatives can also be optimized in some cases because they will be constant and so only need to be calculated once. Generalized linear models are implemented in a way that allows fast calculation of derivatives. The structure of the algorithm to classify the functional form of logical nodes is very similar to that for derivatives and uses a backward mode scheme","category":"page"},{"location":"BUGS_notes/","page":"Notes on BUGS Implementations","title":"Notes on BUGS Implementations","text":"BUGS separates management of logical and stochastic variables, essentially two graphs. Logical variables are stored in an array and values are updated with values in earlier positions of the array.","category":"page"},{"location":"differences/#Differences-Between-BUGS-and-JuliaBUGS","page":"Differences Between BUGS and JuliaBUGS","title":"Differences Between BUGS and JuliaBUGS","text":"","category":"section"},{"location":"differences/#Implicit-Indexing","page":"Differences Between BUGS and JuliaBUGS","title":"Implicit Indexing","text":"","category":"section"},{"location":"differences/","page":"Differences Between BUGS and JuliaBUGS","title":"Differences Between BUGS and JuliaBUGS","text":"In BUGS, x[, ] is used for implicit indexing, which selects all elements from both the first and second dimensions. In JuliaBUGS, users must explicitly use Colon (:) like x[:, :] when using the @bugs macro. The @bugs macro will insert a Colon when given x[], however, the Julia parser will throw an error if given x[, ]. The original BUGS parser will automatically insert a Colon (:) when it encounters x[, ].","category":"page"},{"location":"functions/","page":"Functions","title":"Functions","text":"Most of the functions from BUGS have been implemented. ","category":"page"},{"location":"functions/","page":"Functions","title":"Functions","text":"JuliaBUGS directly utilizes functions from the Julia Standard Library when they share the same names and functionalities. For functions not available in the Julia Standard Library and other popular libraries, we have developed equivalents within JuliaBUGS.BUGSPrimitives.","category":"page"},{"location":"functions/#Function-defined-in-Julia-Standard-Library","page":"Functions","title":"Function defined in Julia Standard Library","text":"","category":"section"},{"location":"functions/","page":"Functions","title":"Functions","text":"warning: Warning\nJuliaBUGS does not allow keyword arguments syntax","category":"page"},{"location":"functions/","page":"Functions","title":"Functions","text":"abs\nexp(x::Real)\nlog(x::Number)\nsqrt(x::Real)\ntrunc\nmin(x::Real, y::Real)\nmax(x::Real, y::Real)\nsum(x::AbstractArray)\nsort(x::AbstractArray)\nsin(x::Real)\ncos(x::Real)\ntan(x::Real)\nasin(x::Real)\nacos(x::Real)\natan(x::Real)\nasinh(x::Real)\nacosh(x::Real)\natanh(x::Real)\nJuliaBUGS.BUGSPrimitives.mean(x::AbstractArray)","category":"page"},{"location":"functions/#Base.abs","page":"Functions","title":"Base.abs","text":"abs(x)\n\nThe absolute value of x.\n\nWhen abs is applied to signed integers, overflow may occur, resulting in the return of a negative value. This overflow occurs only when abs is applied to the minimum representable value of a signed integer. That is, when x == typemin(typeof(x)), abs(x) == x < 0, not -x as might be expected.\n\nSee also: abs2, unsigned, sign.\n\nExamples\n\njulia> abs(-3)\n3\n\njulia> abs(1 + im)\n1.4142135623730951\n\njulia> abs.(Int8[-128 -127 -126 0 126 127]) # overflow at typemin(Int8)\n1×6 Matrix{Int8}:\n -128 127 126 0 126 127\n\njulia> maximum(abs, [1, -2, 3, -4])\n4\n\n\n\n\n\n","category":"function"},{"location":"functions/#Base.exp-Tuple{Real}","page":"Functions","title":"Base.exp","text":"exp(x)\n\nCompute the natural base exponential of x, in other words ℯ^x.\n\nSee also exp2, exp10 and cis.\n\nExamples\n\njulia> exp(1.0)\n2.718281828459045\n\njulia> exp(im * pi) ≈ cis(pi)\ntrue\n\n\n\n\n\n","category":"method"},{"location":"functions/#Base.log-Tuple{Number}","page":"Functions","title":"Base.log","text":"log(x)\n\nCompute the natural logarithm of x. Throws DomainError for negative Real arguments. Use complex negative arguments to obtain complex results.\n\nSee also ℯ, log1p, log2, log10.\n\nExamples\n\njulia> log(2)\n0.6931471805599453\n\njulia> log(-3)\nERROR: DomainError with -3.0:\nlog was called with a negative real argument but will only return a complex result if called with a complex argument. Try log(Complex(x)).\nStacktrace:\n [1] throw_complex_domainerror(::Symbol, ::Float64) at ./math.jl:31\n[...]\n\njulia> log.(exp.(-1:1))\n3-element Vector{Float64}:\n -1.0\n 0.0\n 1.0\n\n\n\n\n\n","category":"method"},{"location":"functions/#Base.sqrt-Tuple{Real}","page":"Functions","title":"Base.sqrt","text":"sqrt(x)\n\nReturn sqrtx. Throws DomainError for negative Real arguments. Use complex negative arguments instead. The prefix operator √ is equivalent to sqrt.\n\nSee also: hypot.\n\nExamples\n\njulia> sqrt(big(81))\n9.0\n\njulia> sqrt(big(-81))\nERROR: DomainError with -81.0:\nNaN result for non-NaN input.\nStacktrace:\n [1] sqrt(::BigFloat) at ./mpfr.jl:501\n[...]\n\njulia> sqrt(big(complex(-81)))\n0.0 + 9.0im\n\njulia> .√(1:4)\n4-element Vector{Float64}:\n 1.0\n 1.4142135623730951\n 1.7320508075688772\n 2.0\n\n\n\n\n\n","category":"method"},{"location":"functions/#Base.trunc","page":"Functions","title":"Base.trunc","text":"trunc([T,] x)\ntrunc(x; digits::Integer= [, base = 10])\ntrunc(x; sigdigits::Integer= [, base = 10])\n\ntrunc(x) returns the nearest integral value of the same type as x whose absolute value is less than or equal to the absolute value of x.\n\ntrunc(T, x) converts the result to type T, throwing an InexactError if the value is not representable.\n\nKeywords digits, sigdigits and base work as for round.\n\nSee also: %, floor, unsigned, unsafe_trunc.\n\nExamples\n\njulia> trunc(2.22)\n2.0\n\njulia> trunc(-2.22, digits=1)\n-2.2\n\njulia> trunc(Int, -2.22)\n-2\n\n\n\n\n\n","category":"function"},{"location":"functions/#Base.min-Tuple{Real, Real}","page":"Functions","title":"Base.min","text":"min(x, y, ...)\n\nReturn the minimum of the arguments (with respect to isless). See also the minimum function to take the minimum element from a collection.\n\nExamples\n\njulia> min(2, 5, 1)\n1\n\n\n\n\n\n","category":"method"},{"location":"functions/#Base.max-Tuple{Real, Real}","page":"Functions","title":"Base.max","text":"max(x, y, ...)\n\nReturn the maximum of the arguments (with respect to isless). See also the maximum function to take the maximum element from a collection.\n\nExamples\n\njulia> max(2, 5, 1)\n5\n\n\n\n\n\n","category":"method"},{"location":"functions/#Base.sum-Tuple{AbstractArray}","page":"Functions","title":"Base.sum","text":"sum(A::AbstractArray; dims)\n\nSum elements of an array over the given dimensions.\n\nExamples\n\njulia> A = [1 2; 3 4]\n2×2 Matrix{Int64}:\n 1 2\n 3 4\n\njulia> sum(A, dims=1)\n1×2 Matrix{Int64}:\n 4 6\n\njulia> sum(A, dims=2)\n2×1 Matrix{Int64}:\n 3\n 7\n\n\n\n\n\n","category":"method"},{"location":"functions/#Base.sort-Tuple{AbstractArray}","page":"Functions","title":"Base.sort","text":"sort(A; dims::Integer, alg::Algorithm=defalg(A), lt=isless, by=identity, rev::Bool=false, order::Ordering=Forward)\n\nSort a multidimensional array A along the given dimension. See sort! for a description of possible keyword arguments.\n\nTo sort slices of an array, refer to sortslices.\n\nExamples\n\njulia> A = [4 3; 1 2]\n2×2 Matrix{Int64}:\n 4 3\n 1 2\n\njulia> sort(A, dims = 1)\n2×2 Matrix{Int64}:\n 1 2\n 4 3\n\njulia> sort(A, dims = 2)\n2×2 Matrix{Int64}:\n 3 4\n 1 2\n\n\n\n\n\n","category":"method"},{"location":"functions/#Base.sin-Tuple{Real}","page":"Functions","title":"Base.sin","text":"sin(x)\n\nCompute sine of x, where x is in radians.\n\nSee also sind, sinpi, sincos, cis, asin.\n\nExamples\n\njulia> round.(sin.(range(0, 2pi, length=9)'), digits=3)\n1×9 Matrix{Float64}:\n 0.0 0.707 1.0 0.707 0.0 -0.707 -1.0 -0.707 -0.0\n\njulia> sind(45)\n0.7071067811865476\n\njulia> sinpi(1/4)\n0.7071067811865475\n\njulia> round.(sincos(pi/6), digits=3)\n(0.5, 0.866)\n\njulia> round(cis(pi/6), digits=3)\n0.866 + 0.5im\n\njulia> round(exp(im*pi/6), digits=3)\n0.866 + 0.5im\n\n\n\n\n\n","category":"method"},{"location":"functions/#Base.cos-Tuple{Real}","page":"Functions","title":"Base.cos","text":"cos(x)\n\nCompute cosine of x, where x is in radians.\n\nSee also cosd, cospi, sincos, cis.\n\n\n\n\n\n","category":"method"},{"location":"functions/#Base.tan-Tuple{Real}","page":"Functions","title":"Base.tan","text":"tan(x)\n\nCompute tangent of x, where x is in radians.\n\n\n\n\n\n","category":"method"},{"location":"functions/#Base.asin-Tuple{Real}","page":"Functions","title":"Base.asin","text":"asin(x)\n\nCompute the inverse sine of x, where the output is in radians.\n\nSee also asind for output in degrees.\n\nExamples\n\njulia> asin.((0, 1/2, 1))\n(0.0, 0.5235987755982989, 1.5707963267948966)\n\njulia> asind.((0, 1/2, 1))\n(0.0, 30.000000000000004, 90.0)\n\n\n\n\n\n","category":"method"},{"location":"functions/#Base.acos-Tuple{Real}","page":"Functions","title":"Base.acos","text":"acos(x)\n\nCompute the inverse cosine of x, where the output is in radians\n\n\n\n\n\n","category":"method"},{"location":"functions/#Base.atan-Tuple{Real}","page":"Functions","title":"Base.atan","text":"atan(y)\natan(y, x)\n\nCompute the inverse tangent of y or y/x, respectively.\n\nFor one argument, this is the angle in radians between the positive x-axis and the point (1, y), returning a value in the interval -pi2 pi2.\n\nFor two arguments, this is the angle in radians between the positive x-axis and the point (x, y), returning a value in the interval -pi pi. This corresponds to a standard atan2 function. Note that by convention atan(0.0,x) is defined as pi and atan(-0.0,x) is defined as -pi when x < 0.\n\nSee also atand for degrees.\n\nExamples\n\njulia> rad2deg(atan(-1/√3))\n-30.000000000000004\n\njulia> rad2deg(atan(-1, √3))\n-30.000000000000004\n\njulia> rad2deg(atan(1, -√3))\n150.0\n\n\n\n\n\n","category":"method"},{"location":"functions/#Base.asinh-Tuple{Real}","page":"Functions","title":"Base.asinh","text":"asinh(x)\n\nCompute the inverse hyperbolic sine of x.\n\n\n\n\n\n","category":"method"},{"location":"functions/#Base.acosh-Tuple{Real}","page":"Functions","title":"Base.acosh","text":"acosh(x)\n\nCompute the inverse hyperbolic cosine of x.\n\n\n\n\n\n","category":"method"},{"location":"functions/#Base.atanh-Tuple{Real}","page":"Functions","title":"Base.atanh","text":"atanh(x)\n\nCompute the inverse hyperbolic tangent of x.\n\n\n\n\n\n","category":"method"},{"location":"functions/#Statistics.mean-Tuple{AbstractArray}","page":"Functions","title":"Statistics.mean","text":"mean(A::AbstractArray; dims)\n\nCompute the mean of an array over the given dimensions.\n\ncompat: Julia 1.1\nmean for empty arrays requires at least Julia 1.1.\n\nExamples\n\njulia> using Statistics\n\njulia> A = [1 2; 3 4]\n2×2 Matrix{Int64}:\n 1 2\n 3 4\n\njulia> mean(A, dims=1)\n1×2 Matrix{Float64}:\n 2.0 3.0\n\njulia> mean(A, dims=2)\n2×1 Matrix{Float64}:\n 1.5\n 3.5\n\n\n\n\n\n","category":"method"},{"location":"functions/#Function-defined-in-[LogExpFunctions](https://github.com/JuliaStats/LogExpFunctions.jl)","page":"Functions","title":"Function defined in LogExpFunctions","text":"","category":"section"},{"location":"functions/","page":"Functions","title":"Functions","text":"cloglog\ncexpexp\nlogit\nlogistic","category":"page"},{"location":"functions/#LogExpFunctions.cloglog","page":"Functions","title":"LogExpFunctions.cloglog","text":"cloglog(x)\n\n\nCompute the complementary log-log, log(-log(1 - x)).\n\n\n\n\n\n","category":"function"},{"location":"functions/#LogExpFunctions.cexpexp","page":"Functions","title":"LogExpFunctions.cexpexp","text":"cexpexp(x)\n\n\nCompute the complementary double exponential, 1 - exp(-exp(x)).\n\n\n\n\n\n","category":"function"},{"location":"functions/#LogExpFunctions.logit","page":"Functions","title":"LogExpFunctions.logit","text":"logit(x)\n\n\nThe logit or log-odds transformation, defined as\n\noperatornamelogit(x) = logleft(fracx1-xright)\n\nfor 0 x 1.\n\nIts inverse is the logistic function.\n\n\n\n\n\n","category":"function"},{"location":"functions/#LogExpFunctions.logistic","page":"Functions","title":"LogExpFunctions.logistic","text":"logistic(x)\n\n\nThe logistic sigmoid function mapping a real number to a value in the interval 01,\n\nsigma(x) = frac1e^-x + 1 = frace^x1+e^x\n\nIts inverse is the logit function.\n\n\n\n\n\n","category":"function"},{"location":"functions/#Function-defined-in-JuliaBUGS.BUGSPrimitives","page":"Functions","title":"Function defined in JuliaBUGS.BUGSPrimitives","text":"","category":"section"},{"location":"functions/","page":"Functions","title":"Functions","text":"JuliaBUGS.BUGSPrimitives.equals\nJuliaBUGS.BUGSPrimitives.inprod\nJuliaBUGS.BUGSPrimitives.inverse\nJuliaBUGS.BUGSPrimitives.logdet\nJuliaBUGS.BUGSPrimitives.logfact\nJuliaBUGS.BUGSPrimitives.loggam\nJuliaBUGS.BUGSPrimitives.icloglog\nJuliaBUGS.BUGSPrimitives.mexp\nJuliaBUGS.BUGSPrimitives.phi\nJuliaBUGS.BUGSPrimitives.pow\nJuliaBUGS.BUGSPrimitives.rank\nJuliaBUGS.BUGSPrimitives.ranked\nJuliaBUGS.BUGSPrimitives.sd\nJuliaBUGS.BUGSPrimitives.softplus\nJuliaBUGS.BUGSPrimitives._step\nJuliaBUGS.BUGSPrimitives.arcsin\nJuliaBUGS.BUGSPrimitives.arcsinh\nJuliaBUGS.BUGSPrimitives.arccos\nJuliaBUGS.BUGSPrimitives.arccosh\nJuliaBUGS.BUGSPrimitives.arctan\nJuliaBUGS.BUGSPrimitives.arctanh","category":"page"},{"location":"functions/#JuliaBUGS.BUGSPrimitives.equals","page":"Functions","title":"JuliaBUGS.BUGSPrimitives.equals","text":"equals(x, y)\n\nReturns 1 if x is equal to y, 0 otherwise.\n\n\n\n\n\n","category":"function"},{"location":"functions/#JuliaBUGS.BUGSPrimitives.inprod","page":"Functions","title":"JuliaBUGS.BUGSPrimitives.inprod","text":"inprod(a, b)\n\nInner product of a and b.\n\n\n\n\n\n","category":"function"},{"location":"functions/#JuliaBUGS.BUGSPrimitives.inverse","page":"Functions","title":"JuliaBUGS.BUGSPrimitives.inverse","text":"inverse(m::AbstractMatrix)\n\nInverse of matrix mathbfm.\n\n\n\n\n\n","category":"function"},{"location":"functions/#JuliaBUGS.BUGSPrimitives.logdet","page":"Functions","title":"JuliaBUGS.BUGSPrimitives.logdet","text":"logdet(::AbstractMatrix)\n\nLogarithm of the determinant of matrix mathbfv.\n\n\n\n\n\n","category":"function"},{"location":"functions/#JuliaBUGS.BUGSPrimitives.logfact","page":"Functions","title":"JuliaBUGS.BUGSPrimitives.logfact","text":"logfact(x)\n\nLogarithm of the factorial of x.\n\n\n\n\n\n","category":"function"},{"location":"functions/#JuliaBUGS.BUGSPrimitives.loggam","page":"Functions","title":"JuliaBUGS.BUGSPrimitives.loggam","text":"loggam(x)\n\nLogarithm of the gamma function of x.\n\n\n\n\n\n","category":"function"},{"location":"functions/#JuliaBUGS.BUGSPrimitives.icloglog","page":"Functions","title":"JuliaBUGS.BUGSPrimitives.icloglog","text":"icloglog(x)\n\nInverse complementary log-log function of x. Alias for cexpexp(x).\n\n\n\n\n\n","category":"function"},{"location":"functions/#JuliaBUGS.BUGSPrimitives.mexp","page":"Functions","title":"JuliaBUGS.BUGSPrimitives.mexp","text":"mexp(x::AbstractMatrix)\n\nMatrix exponential of mathbfx.\n\n\n\n\n\n","category":"function"},{"location":"functions/#JuliaBUGS.BUGSPrimitives.phi","page":"Functions","title":"JuliaBUGS.BUGSPrimitives.phi","text":"phi(x)\n\nCumulative distribution function (CDF) of the standard normal distribution evaluated at x.\n\n\n\n\n\n","category":"function"},{"location":"functions/#JuliaBUGS.BUGSPrimitives.pow","page":"Functions","title":"JuliaBUGS.BUGSPrimitives.pow","text":"pow(a, b)\n\nReturn a raised to the power of b.\n\n\n\n\n\n","category":"function"},{"location":"functions/#JuliaBUGS.BUGSPrimitives.rank","page":"Functions","title":"JuliaBUGS.BUGSPrimitives.rank","text":"rank(v::AbstractVector, i::Integer)\n\nReturn the rank of the i-th element of mathbfv.\n\n\n\n\n\n","category":"function"},{"location":"functions/#JuliaBUGS.BUGSPrimitives.ranked","page":"Functions","title":"JuliaBUGS.BUGSPrimitives.ranked","text":"ranked(v::AbstractVector, i::Integer)\n\nReturn the i-th element of mathbfv sorted in ascending order.\n\n\n\n\n\n","category":"function"},{"location":"functions/#JuliaBUGS.BUGSPrimitives.sd","page":"Functions","title":"JuliaBUGS.BUGSPrimitives.sd","text":"sd(v::AbstractVector)\n\nReturn the standard deviation of the input vector mathbfv.\n\n\n\n\n\n","category":"function"},{"location":"functions/#JuliaBUGS.BUGSPrimitives.softplus","page":"Functions","title":"JuliaBUGS.BUGSPrimitives.softplus","text":"softplus(x)\n\nReturn the softplus function of x, defined as log(1 + exp(x)).\n\n\n\n\n\n","category":"function"},{"location":"functions/#JuliaBUGS.BUGSPrimitives._step","page":"Functions","title":"JuliaBUGS.BUGSPrimitives._step","text":"_step(x)\n\nReturn 1 if x is greater than 0, and 0 otherwise.\n\n\n\n\n\n","category":"function"},{"location":"functions/#JuliaBUGS.BUGSPrimitives.arcsin","page":"Functions","title":"JuliaBUGS.BUGSPrimitives.arcsin","text":"arcsin(x)\n\nReturn the arcsine of x.\n\n\n\n\n\n","category":"function"},{"location":"functions/#JuliaBUGS.BUGSPrimitives.arcsinh","page":"Functions","title":"JuliaBUGS.BUGSPrimitives.arcsinh","text":"arcsinh(x)\n\nReturn the inverse hyperbolic sine of x.\n\n\n\n\n\n","category":"function"},{"location":"functions/#JuliaBUGS.BUGSPrimitives.arccos","page":"Functions","title":"JuliaBUGS.BUGSPrimitives.arccos","text":"arccos(x)\n\nReturn the arccosine of x.\n\n\n\n\n\n","category":"function"},{"location":"functions/#JuliaBUGS.BUGSPrimitives.arccosh","page":"Functions","title":"JuliaBUGS.BUGSPrimitives.arccosh","text":"arccosh(x)\n\nReturn the inverse hyperbolic cosine of x.\n\n\n\n\n\n","category":"function"},{"location":"functions/#JuliaBUGS.BUGSPrimitives.arctan","page":"Functions","title":"JuliaBUGS.BUGSPrimitives.arctan","text":"arctan(x)\n\nReturn the arctangent of x.\n\n\n\n\n\n","category":"function"},{"location":"functions/#JuliaBUGS.BUGSPrimitives.arctanh","page":"Functions","title":"JuliaBUGS.BUGSPrimitives.arctanh","text":"arctanh(x)\n\nSee atanh.\n\n\n\n\n\n","category":"function"},{"location":"user_defined_functions/#Define-and-Use-Your-Own-Functions-and-Distributions","page":"User-defined Functions and Distributions","title":"Define and Use Your Own Functions and Distributions","text":"","category":"section"},{"location":"user_defined_functions/","page":"User-defined Functions and Distributions","title":"User-defined Functions and Distributions","text":"For now, out of the box, JuliaBUGS only allows functions and distributions defined in BUGSPrimitives to be used in the model. With the @register_primitive macro, users can register their own functions and distributions with JuliaBUGS. It is important to ensure that any functions used are pure mathematical functions. This implies that such functions should not alter any external state including but not limited to modifying global variables, writing data to files. (Printing might be okay, but do at discretion.)","category":"page"},{"location":"user_defined_functions/","page":"User-defined Functions and Distributions","title":"User-defined Functions and Distributions","text":"julia> JuliaBUGS.@register_primitive function f(x)\n return x + 1\nend\nf (generic function with 1 method)\n\njulia> JuliaBUGS.f(2)\n3","category":"page"},{"location":"user_defined_functions/","page":"User-defined Functions and Distributions","title":"User-defined Functions and Distributions","text":"Users can also introduce a function into JuliaBUGS, by ","category":"page"},{"location":"user_defined_functions/","page":"User-defined Functions and Distributions","title":"User-defined Functions and Distributions","text":"julia> f(x) = x + 1\nf (generic function with 1 method)\n\njulia> JuliaBUGS.@register_primitive(f);\n\njulia> JuliaBUGS.f(1)\n2","category":"page"},{"location":"user_defined_functions/","page":"User-defined Functions and Distributions","title":"User-defined Functions and Distributions","text":"After registering the function or distributions, they can be used just like any other functions or distributions provided by BUGS.","category":"page"},{"location":"parser/#BUGS-Parser","page":"Parser","title":"BUGS Parser","text":"","category":"section"},{"location":"parser/","page":"Parser","title":"Parser","text":"The macro @bugs produces a Julia Expr object that represents the BUGS model definition.","category":"page"},{"location":"parser/","page":"Parser","title":"Parser","text":"If the input is a String, it's assumed to be a program in the original BUGS language. In this case, the macro will first convert the program to an equivalent Julia program, then use the Julia parser to parse the program into an Expr object.","category":"page"},{"location":"parser/","page":"Parser","title":"Parser","text":"Both model definitions written in Julia and those written in the original BUGS and subsequently parsed are now represented as a Julia Expr object. These objects go through syntax checking and post-processing to create the input for the compile function.","category":"page"},{"location":"parser/","page":"Parser","title":"Parser","text":"Below, we describe how the original BUGS program is translated to an equivalent Julia program and detail the post-processing done to the Expr object.","category":"page"},{"location":"parser/#BUGS-to-Julia-Translation","page":"Parser","title":"BUGS to Julia Translation","text":"","category":"section"},{"location":"parser/","page":"Parser","title":"Parser","text":"In this section, we refer to the translation program as the \"parser\" and the translating process as \"parsing\". Although the parser doesn't produce a syntax tree, it does follow the form of a recursive descent parser, building a Julia program in the form of a vector of tokens rather than a syntax tree.","category":"page"},{"location":"parser/","page":"Parser","title":"Parser","text":"This general implementation is heavily inspired by JuliaSyntax.jl, the official parser for Julia since version 1.10.","category":"page"},{"location":"parser/","page":"Parser","title":"Parser","text":"The BUGS parser implemented here takes a token stream with a recursive descent structure and checks the program's correctness. Here's how it works:","category":"page"},{"location":"parser/","page":"Parser","title":"Parser","text":"Use tokenize to obtain the token vector.\nInspect the tokens and build the Julia version of the program as a vector of tokens.\nPush the token to the Julia version of the program vector when appropriate.\nDetect errors and make necessary alterations to tokens, such as deletion, combination, or replacement.","category":"page"},{"location":"parser/","page":"Parser","title":"Parser","text":"During the recursive descent, BUGS syntax tokens will be translated into Julia syntax tokens. Some tokens will remain as they are, while others will be transformed, removed, or new tokens may be added.","category":"page"},{"location":"parser/","page":"Parser","title":"Parser","text":"The parser will throw an error if it encounters a program that does not adhere to strict BUGS syntax.","category":"page"},{"location":"parser/#Some-Notes-on-Error-Recovery","page":"Parser","title":"Some Notes on Error Recovery","text":"","category":"section"},{"location":"parser/","page":"Parser","title":"Parser","text":"The current error recovery is ad hoc and primarily rudimentary. If the program is correct, it will produce the correct result. If the program is syntactically or semantically incorrect, the token stream will not be pushed forward, resulting in failure.","category":"page"},{"location":"parser/","page":"Parser","title":"Parser","text":"The failure detection mechanism checks if two errors occur with the same \"current token\". If they do, the parser stops and reports the error. This ensures that the parser won't incorrectly parse a flawed program.","category":"page"},{"location":"bugs_lang/#Special-Cases-in-the-BUGS-Language","page":"Special Cases in the BUGS Language","title":"Special Cases in the BUGS Language","text":"","category":"section"},{"location":"bugs_lang/","page":"Special Cases in the BUGS Language","title":"Special Cases in the BUGS Language","text":"Here we record some of the special cases in the BUGS language in the original BUGS softwares, JuliaBUGS may or may not inherent these behaviors.","category":"page"},{"location":"bugs_lang/#Nested-Indexing-on-the-Left-Hand-Side","page":"Special Cases in the BUGS Language","title":"Nested Indexing on the Left-Hand Side","text":"","category":"section"},{"location":"bugs_lang/","page":"Special Cases in the BUGS Language","title":"Special Cases in the BUGS Language","text":"model_def = @bugs x[y[1]] ~ dnorm(0, 1) # this is permitted\ndata = (y=[1,2,3],)","category":"page"},{"location":"bugs_lang/#Function-Application-on-the-Left-Hand-Side-Index-is-Not-Permitted","page":"Special Cases in the BUGS Language","title":"Function Application on the Left-Hand Side Index is Not Permitted","text":"","category":"section"},{"location":"bugs_lang/","page":"Special Cases in the BUGS Language","title":"Special Cases in the BUGS Language","text":"This is identified during the compilation phase, not during parsing.","category":"page"},{"location":"bugs_lang/#The-arguments-to-the-distributions-functions-must-be-\"simple\":-they-are-either-a-variable-or-a-constant,-their-should-be-no-function-application-other-than-the-arithmetic-operators","page":"Special Cases in the BUGS Language","title":"The arguments to the distributions functions must be \"simple\": they are either a variable or a constant, their should be no function application other than the arithmetic operators","text":"","category":"section"},{"location":"bugs_lang/","page":"Special Cases in the BUGS Language","title":"Special Cases in the BUGS Language","text":"model {\n x ~ dnorm(y[1] + 1, 1)\n} # fail at parsing\n\nmodel {\n x ~ dnorm(sum(y[1:2]), 1)\n} # fail at parsing\n\nmodel {\n x ~ dnorm(y[sum(y[1:2])], 1)\n} # pass the parser, but fail at compile\n\nmodel {\n x ~ dnorm(y[y[2]], 1)\n} # this is okay\n\nmodel {\n x ~ dnorm(y[y[2]+1], 1)\n} # this is also okay\n\nlist(y = c(1, 2, 3))","category":"page"},{"location":"bugs_lang/","page":"Special Cases in the BUGS Language","title":"Special Cases in the BUGS Language","text":"The reason behind this is two fold: (1) forcing the user to use ~ to identify the dependency between random variables, instead of expressions, (2) make implementation of automatic differentiation easier.","category":"page"},{"location":"bugs_lang/#Indexing-with-Transformed-Variables-is-Not-Permitted","page":"Special Cases in the BUGS Language","title":"Indexing with Transformed Variables is Not Permitted","text":"","category":"section"},{"location":"bugs_lang/","page":"Special Cases in the BUGS Language","title":"Special Cases in the BUGS Language","text":"This restriction is due to the current compiler implementation.","category":"page"},{"location":"bugs_lang/","page":"Special Cases in the BUGS Language","title":"Special Cases in the BUGS Language","text":"model{\n a <- max(y[2], y[3])\n x[a + 1] ~ dnorm(0, 1)\n}\n\nlist(y=c(1, 2, 3))","category":"page"},{"location":"bugs_lang/","page":"Special Cases in the BUGS Language","title":"Special Cases in the BUGS Language","text":"The reason for this is that the transformed variable and indices are evaluated simultaneously, and the transformed variable is not available at the time of indexing. This is a limitation of the current implementation.","category":"page"},{"location":"api/#API","page":"General","title":"API","text":"","category":"section"},{"location":"api/","page":"General","title":"General","text":"@bugs\ncompile\nBUGSModel\nConcreteNodeInfo\nBUGSGraph","category":"page"},{"location":"api/#JuliaBUGS.Parser.@bugs","page":"General","title":"JuliaBUGS.Parser.@bugs","text":"@bugs(prog::String, replace_period=true, no_enclosure=false)\n\nProduce similar output as @bugs, but takes a string as input. This is useful for parsing original BUGS programs.\n\nArguments\n\nprog::String: The BUGS program code as a string.\nreplace_period::Bool: If true, periods in the BUGS code will be replaced (default true).\nno_enclosure::Bool: If true, the parser will not expect the program to be wrapped between model{ } (default false).\n\n\n\n\n\n","category":"macro"},{"location":"api/#JuliaBUGS.compile","page":"General","title":"JuliaBUGS.compile","text":"compile(model_def[, data, initializations])\n\nCompile a BUGS model into a log density problem.\n\nArguments\n\nmodel_def::Expr: The BUGS model definition.\ndata::NamedTuple or AbstractDict: The data to be used in the model. If none is passed, the data will be assumed to be empty.\ninitializations::NamedTuple or AbstractDict: The initial values for the model parameters. If none is passed, the parameters will be assumed to be initialized to zero.\nis_transformed::Bool=true: If true, the model parameters during inference will be transformed to the unconstrained space. \n\nReturns\n\nA BUGSModel object representing the compiled model.\n\n\n\n\n\n","category":"function"},{"location":"api/#JuliaBUGS.BUGSModel","page":"General","title":"JuliaBUGS.BUGSModel","text":"BUGSModel\n\nThe BUGSModel object is used for inference and represents the output of compilation. It fully implements the LogDensityProblems.jl interface.\n\nFields\n\ntransformed::Bool: Indicates whether the model parameters are in the transformed space.\nuntransformed_param_length::Int: The length of the parameters vector in the original space.\ntransformed_param_length::Int: The length of the parameters vector in the transformed space.\nuntransformed_var_lengths::Dict{VarName,Int}: A dictionary mapping the names of the variables to their lengths in the original space.\ntransformed_var_lengths::Dict{VarName,Int}: A dictionary mapping the names of the variables to their lengths in the transformed space.\nvarinfo::SimpleVarInfo: An instance of DynamicPPL.SimpleVarInfo, which is a dictionary-like data structure that maps both data and values of variables in the model to the corresponding values.\nparameters::Vector{VarName}: A vector containing the names of the parameters in the model, defined as stochastic variables that are not observed. This vector should be consistent with sorted_nodes.\nsorted_nodes::Vector{VarName}: A vector containing the names of all the variables in the model, sorted in topological order. In the case of a conditioned model, sorted_nodes include all the variables in parameters and the variables in the Markov blanket of parameters.\ng::BUGSGraph: An instance of BUGSGraph, representing the dependency graph of the model.\nbase_model::Union{BUGSModel,Nothing}: If not Nothing, the model is a conditioned model; otherwise, it's the model returned by compile.\n\n\n\n\n\n","category":"type"},{"location":"api/#JuliaBUGS.ConcreteNodeInfo","page":"General","title":"JuliaBUGS.ConcreteNodeInfo","text":"ConcreteNodeInfo\n\nDefines the information stored in each node of the BUGS graph, encapsulating the essential characteristics and functions associated with a node within the BUGS model's dependency graph.\n\nFields\n\nnode_type::VariableTypes: Specifies whether the node is a stochastic or logical variable.\nnode_function_expr::Expr: The node function expression.\nnode_args::Vector{VarName}: A vector containing the names of the variables that are arguments to the node function.\n\n\n\n\n\n","category":"type"},{"location":"api/#JuliaBUGS.BUGSGraph","page":"General","title":"JuliaBUGS.BUGSGraph","text":"BUGSGraph\n\nThe BUGSGraph object represents the graph structure for a BUGS model. It is a type alias for MetaGraphsNext.MetaGraph with node type specified to ConcreteNodeInfo.\n\n\n\n\n\n","category":"type"},{"location":"graph_plotting/#Plotting-graphs","page":"Plotting","title":"Plotting graphs","text":"","category":"section"},{"location":"graph_plotting/","page":"Plotting","title":"Plotting","text":"Plotting the graph can be very beneficial for debugging the model.","category":"page"},{"location":"graph_plotting/","page":"Plotting","title":"Plotting","text":"Note Plate notation is not yet supported. Therefore, it's advisable for users to begin with a more streamlined model that contains fewer nodes, allowing for clearer visualization.","category":"page"},{"location":"graph_plotting/","page":"Plotting","title":"Plotting","text":"In Julia, we've set up standard plotting routines using various graphing libraries. You can visualize graphs with three different libraries by employing a common model, as detailed below:","category":"page"},{"location":"graph_plotting/","page":"Plotting","title":"Plotting","text":"model_def = @bugs begin\n a ~ dnorm(f, c)\n f = b - 1\n b ~ dnorm(0, 1)\n c ~ dnorm(l, 1)\n g = a * 2\n d ~ dnorm(g, 1)\n h = g + 2\n e ~ dnorm(h, i)\n i ~ dnorm(0, 1)\n l ~ dnorm(0, 1)\nend\n\ninits = Dict(\n :a => 1.0,\n :b => 2.0,\n :c => 3.0,\n :d => 4.0,\n :e => 5.0,\n\n # :f => 1.0,\n # :g => 2.0,\n # :h => 4.0,\n\n :i => 4.0,\n :l => -2.0,\n)\n\nmodel = compile(model_def, NamedTuple(), inits)","category":"page"},{"location":"graph_plotting/#[TikzGraphs.jl](https://github.com/JuliaTeX/TikzGraphs.jl).","page":"Plotting","title":"TikzGraphs.jl.","text":"","category":"section"},{"location":"graph_plotting/","page":"Plotting","title":"Plotting","text":"using TikzGraphs\nTikzGraphs.plot(model)","category":"page"},{"location":"graph_plotting/","page":"Plotting","title":"Plotting","text":"(Image: TikzGraphs)","category":"page"},{"location":"graph_plotting/#[GraphPlot.jl](https://github.com/JuliaGraphs/GraphPlot.jl)","page":"Plotting","title":"GraphPlot.jl","text":"","category":"section"},{"location":"graph_plotting/","page":"Plotting","title":"Plotting","text":"using GraphPlot\ngplot(model)","category":"page"},{"location":"graph_plotting/","page":"Plotting","title":"Plotting","text":"(Image: GraphPlot)","category":"page"},{"location":"graph_plotting/#[GraphMakie.jl](https://github.com/MakieOrg/GraphMakie.jl)","page":"Plotting","title":"GraphMakie.jl","text":"","category":"section"},{"location":"graph_plotting/","page":"Plotting","title":"Plotting","text":"using GLMakie, GraphMakie\ngraphplot(model)","category":"page"},{"location":"graph_plotting/","page":"Plotting","title":"Plotting","text":"(Image: GraphMakie)","category":"page"},{"location":"bugs_examples/Leuk/#Leuk:-Cox-regression","page":"Leuk: Cox regression","title":"Leuk: Cox regression","text":"","category":"section"},{"location":"bugs_examples/Leuk/#Description","page":"Leuk: Cox regression","title":"Description","text":"","category":"section"},{"location":"bugs_examples/Leuk/","page":"Leuk: Cox regression","title":"Leuk: Cox regression","text":"Several authors have discussed Bayesian inference for censored survival data where the integrated baseline hazard function is to be estimated non-parametrically, including Kalbfleisch (1978), Kalbfleisch and Prentice (1980), Clayton (1991), and Clayton (1994). Clayton (1994) formulates the Cox model using counting process notation introduced by Andersen and Gill (1982) and discusses estimation of the baseline hazard and regression parameters using MCMC methods. Although this approach may seem somewhat contrived, it lays the groundwork for extensions to random effect (frailty) models, time-dependent covariates, smoothed hazards, multiple events, and more. Below is how to implement this formulation of the Cox model in BUGS.","category":"page"},{"location":"bugs_examples/Leuk/","page":"Leuk: Cox regression","title":"Leuk: Cox regression","text":"For subjects i = 1n, we observe processes N_i(t) which count the number of failures up to time t. The corresponding intensity process I_i(t) is given by","category":"page"},{"location":"bugs_examples/Leuk/","page":"Leuk: Cox regression","title":"Leuk: Cox regression","text":"I_i(t) dt = mathbbEdN_i(t) mathcalF_t-","category":"page"},{"location":"bugs_examples/Leuk/","page":"Leuk: Cox regression","title":"Leuk: Cox regression","text":"where dN_i(t) is the increment of N_i over the interval t t+dt), and mathcalF_t- represents the available data just before time t. If subject i fails during this interval, dN_i(t) = 1; otherwise, dN_i(t) = 0. Hence, mathbbE(dN_i(t) mathcalF_t-) corresponds to the probability of subject i failing in the interval t t+dt). As dt to 0, this probability becomes the instantaneous hazard at time t for subject i, assumed to have the form","category":"page"},{"location":"bugs_examples/Leuk/","page":"Leuk: Cox regression","title":"Leuk: Cox regression","text":"I_i(t) = Y_i(t)lambda_0(t) exp(beta z_i)","category":"page"},{"location":"bugs_examples/Leuk/","page":"Leuk: Cox regression","title":"Leuk: Cox regression","text":"where Y_i(t) is an observed process taking the value 1 or 0 according to whether subject i is at risk at time t, and lambda_0(t) exp(beta z_i) is the Cox regression model. Thus, the observed data D = N_i(t) Y_i(t) z_i i = 1n and unknown parameters beta and Lambda_0(t) = int_0^t lambda_0(u) du, the latter estimated non-parametrically.","category":"page"},{"location":"bugs_examples/Leuk/","page":"Leuk: Cox regression","title":"Leuk: Cox regression","text":"The joint posterior distribution is","category":"page"},{"location":"bugs_examples/Leuk/","page":"Leuk: Cox regression","title":"Leuk: Cox regression","text":"P(beta Lambda_0() D) propto P(D beta Lambda_0()) P(beta) P(Lambda_0())","category":"page"},{"location":"bugs_examples/Leuk/","page":"Leuk: Cox regression","title":"Leuk: Cox regression","text":"For BUGS, specify the likelihood P(D beta Lambda_0()) and priors for beta and Lambda_0(). Under non-informative censoring, the data likelihood is","category":"page"},{"location":"bugs_examples/Leuk/","page":"Leuk: Cox regression","title":"Leuk: Cox regression","text":"prod_i=1^n left( prod_t geq 0 I_i(t) dN_i(t) right) exp(- I_i(t) dt)","category":"page"},{"location":"bugs_examples/Leuk/","page":"Leuk: Cox regression","title":"Leuk: Cox regression","text":"Viewing the increments dN_i(t) as independent Poisson variables with means I_i(t)dt:","category":"page"},{"location":"bugs_examples/Leuk/","page":"Leuk: Cox regression","title":"Leuk: Cox regression","text":"dN_i(t) sim textPoisson(I_i(t)dt)","category":"page"},{"location":"bugs_examples/Leuk/","page":"Leuk: Cox regression","title":"Leuk: Cox regression","text":"I_i(t)dt = Y_i(t) exp(beta z_i) dLambda_0(t)","category":"page"},{"location":"bugs_examples/Leuk/","page":"Leuk: Cox regression","title":"Leuk: Cox regression","text":"where (d\\Lambda0(t) = \\Lambda0(t)dt) is the increment or jump in the integrated baseline hazard function occurring during the time interval ([t, t+dt)). Since the conjugate prior for the Poisson mean is the gamma distribution, it would be convenient if (\\Lambda0()) were a process in which the increments (d\\Lambda0(t)) are distributed according to gamma distributions. We assume the conjugate independent increments prior suggested by Kalbfleisch (1978), namely ","category":"page"},{"location":"bugs_examples/Leuk/","page":"Leuk: Cox regression","title":"Leuk: Cox regression","text":"dLambda_0(t) sim textGamma(c cdot dLambda^*_0(t) c)","category":"page"},{"location":"bugs_examples/Leuk/","page":"Leuk: Cox regression","title":"Leuk: Cox regression","text":"Here, dLambda^*_0(t) can be thought of as a prior guess at the unknown hazard function, with c representing the degree of confidence in this guess. Small values of c correspond to weak prior beliefs. In the example below, we set dLambda^*_0(t) = r cdot dt where r is a guess at the failure rate per unit time, and dt is the size of the time interval. ","category":"page"},{"location":"bugs_examples/Leuk/","page":"Leuk: Cox regression","title":"Leuk: Cox regression","text":"The above formulation is appropriate when genuine prior information exists concerning the underlying hazard function. Alternatively, if we wish to reproduce a Cox analysis but with, say, additional hierarchical structure, we may use the Multinomial-Poisson trick described in the BUGS manual. This is equivalent to assuming independent increments in the cumulative non-informative priors. This formulation is also shown below.","category":"page"},{"location":"bugs_examples/Leuk/","page":"Leuk: Cox regression","title":"Leuk: Cox regression","text":"The fixed effect regression coefficients (b) are assigned a vague prior","category":"page"},{"location":"bugs_examples/Leuk/","page":"Leuk: Cox regression","title":"Leuk: Cox regression","text":"b sim textNormal(00 0000001)","category":"page"},{"location":"bugs_examples/Leuk/#BUGS-code-for-the-Leuk-example:","page":"Leuk: Cox regression","title":"BUGS code for the Leuk example:","text":"","category":"section"},{"location":"bugs_examples/Leuk/","page":"Leuk: Cox regression","title":"Leuk: Cox regression","text":"model\n{\n # Set up data\n for(i in 1:N) {\n for(j in 1:T) {\n # risk set = 1 if obs.t >= t\n Y[i,j] <- step(obs.t[i] - t[j] + eps)\n # counting process jump = 1 if obs.t in [ t[j], t[j+1] )\n # i.e. if t[j] <= obs.t < t[j+1]\n dN[i, j] <- Y[i, j] * step(t[j + 1] - obs.t[i] - eps) * fail[i]\n }\n }\n # Model\n for(j in 1:T) {\n for(i in 1:N) {\n dN[i, j] ~ dpois(Idt[i, j]) # Likelihood\n Idt[i, j] <- Y[i, j] * exp(beta * Z[i]) * dL0[j] # Intensity\n }\n dL0[j] ~ dgamma(mu[j], c)\n mu[j] <- dL0.star[j] * c # prior mean hazard\n\n # Survivor function = exp(-Integral{l0(u)du})^exp(beta*z)\n S.treat[j] <- pow(exp(-sum(dL0[1 : j])), exp(beta * -0.5));\n S.placebo[j] <- pow(exp(-sum(dL0[1 : j])), exp(beta * 0.5)); \n }\n # Priors\n c <- 0.001\n r <- 0.1\n for (j in 1 : T) {\n dL0.star[j] <- r * (t[j + 1] - t[j])\n }\n beta ~ dnorm(0.0,0.000001)\n}","category":"page"},{"location":"bugs_examples/Leuk/#Reference-Results","page":"Leuk: Cox regression","title":"Reference Results","text":"","category":"section"},{"location":"bugs_examples/Leuk/","page":"Leuk: Cox regression","title":"Leuk: Cox regression","text":"Variable Mean Median Standard Deviation Monte Carlo Error 2.5% Value 97.5% Value Start Sample ESS\nS.placebo[1] 0.9264 0.9374 0.04989 3.349E-4 0.8029 0.9909 1001 20000 22184\nS.placebo[17] 0.04431 0.03344 0.03909 2.698E-4 0.002478 0.1487 1001 20000 20992\nS.treat[1] 0.9826 0.9863 0.01413 1.074E-4 0.9457 0.9982 1001 20000 17315\nS.treat[17] 0.4767 0.4763 0.1198 0.001009 0.2474 0.7086 1001 20000 14104\nbeta 1.539 1.524 0.4211 0.0034 0.7475 2.388 1001 20000 15340","category":"page"},{"location":"R_interface/#Integrating-R-in-Julia","page":"R Interface","title":"Integrating R in Julia","text":"","category":"section"},{"location":"R_interface/","page":"R Interface","title":"R Interface","text":"Julia offers a seamless interface to the R language. ","category":"page"},{"location":"R_interface/","page":"R Interface","title":"R Interface","text":"The RCall.jl package enables interaction with R functions in Julia.\nThe RData.jl package allows interfacing with R data in Julia.","category":"page"},{"location":"R_interface/#Reading-BUGS-data-and-init-from-R-like-lists","page":"R Interface","title":"Reading BUGS data and init from R like lists","text":"","category":"section"},{"location":"R_interface/","page":"R Interface","title":"R Interface","text":"Warning: The data layout in BUGS assumes that the data is stored in row-major order, while R uses column-major order. This discrepancy can lead to issues. Stan developers have transformed the data and initializations of example BUGS models for R, which can be found here.","category":"page"},{"location":"R_interface/#Reading-the-list-data-structure-from-R","page":"R Interface","title":"Reading the list data structure from R","text":"","category":"section"},{"location":"R_interface/","page":"R Interface","title":"R Interface","text":"The data for Rats is available here. ","category":"page"},{"location":"R_interface/","page":"R Interface","title":"R Interface","text":"In Julia, we can read this data into a Julia dictionary using the RCall.jl package.","category":"page"},{"location":"R_interface/","page":"R Interface","title":"R Interface","text":"julia> using RCall\n\njulia> data = R\"\nlist(\n x = c(8.0, 15.0, 22.0, 29.0, 36.0), xbar = 22, N = 30, T = 5,\n Y = structure(\n .Data = c(\n 151, 199, 246, 283, 320,\n 145, 199, 249, 293, 354,\n 147, 214, 263, 312, 328,\n 155, 200, 237, 272, 297,\n 135, 188, 230, 280, 323,\n 159, 210, 252, 298, 331,\n 141, 189, 231, 275, 305,\n 159, 201, 248, 297, 338,\n 177, 236, 285, 350, 376,\n 134, 182, 220, 260, 296,\n 160, 208, 261, 313, 352,\n 143, 188, 220, 273, 314,\n 154, 200, 244, 289, 325,\n 171, 221, 270, 326, 358,\n 163, 216, 242, 281, 312,\n 160, 207, 248, 288, 324,\n 142, 187, 234, 280, 316,\n 156, 203, 243, 283, 317,\n 157, 212, 259, 307, 336,\n 152, 203, 246, 286, 321,\n 154, 205, 253, 298, 334,\n 139, 190, 225, 267, 302,\n 146, 191, 229, 272, 302,\n 157, 211, 250, 285, 323,\n 132, 185, 237, 286, 331,\n 160, 207, 257, 303, 345,\n 169, 216, 261, 295, 333,\n 157, 205, 248, 289, 316,\n 137, 180, 219, 258, 291,\n 153, 200, 244, 286, 324\n ),\n .Dim = c(30, 5)\n )\n)\n\"\nRObject{VecSxp}\n$x\n[1] 8 15 22 29 36\n\n$xbar\n[1] 22\n\n$N\n[1] 30\n\n$T\n[1] 5\n\n$Y\n [,1] [,2] [,3] [,4] [,5]\n [1,] 151 141 154 157 132\n [2,] 199 189 200 212 185\n [3,] 246 231 244 259 237\n [4,] 283 275 289 307 286\n [5,] 320 305 325 336 331\n [6,] 145 159 171 152 160\n [7,] 199 201 221 203 207\n [8,] 249 248 270 246 257\n [9,] 293 297 326 286 303\n[10,] 354 338 358 321 345\n[11,] 147 177 163 154 169\n[12,] 214 236 216 205 216\n[13,] 263 285 242 253 261\n[14,] 312 350 281 298 295\n[15,] 328 376 312 334 333\n[16,] 155 134 160 139 157\n[17,] 200 182 207 190 205\n[18,] 237 220 248 225 248\n[19,] 272 260 288 267 289\n[20,] 297 296 324 302 316\n[21,] 135 160 142 146 137\n[22,] 188 208 187 191 180\n[23,] 230 261 234 229 219\n[24,] 280 313 280 272 258\n[25,] 323 352 316 302 291\n[26,] 159 143 156 157 153\n[27,] 210 188 203 211 200\n[28,] 252 220 243 250 244\n[29,] 298 273 283 285 286\n[30,] 331 314 317 323 324","category":"page"},{"location":"R_interface/","page":"R Interface","title":"R Interface","text":"alternatively, reval(s::String) will produce the same result in this case.","category":"page"},{"location":"R_interface/","page":"R Interface","title":"R Interface","text":"If the data is stores in a file, user can use function (may require customizing the function to fit specific needs)","category":"page"},{"location":"R_interface/","page":"R Interface","title":"R Interface","text":"function read_rlist_to_dictionary(filepath::String)\n r_data = open(filepath) do f\n s = read(f, String)\n reval(s)\n end\n return rcopy(r_data)\nend","category":"page"},{"location":"R_interface/","page":"R Interface","title":"R Interface","text":", and save the result to a Julia variable and access the data as a Julia dictionary","category":"page"},{"location":"R_interface/","page":"R Interface","title":"R Interface","text":"julia> rcopy(data)\nOrderedDict{Symbol, Any} with 5 entries:\n :x => [8.0, 15.0, 22.0, 29.0, 36.0]\n :xbar => 22.0\n :N => 30.0\n :T => 5.0\n :Y => [151.0 141.0 … 157.0 132.0; 199.0 189.0 … 212.0 185.0; … ; 298.0 273.0 … 285.0 286.0; 331.0 314.0 … 323.0 324.0]","category":"page"},{"location":"R_interface/","page":"R Interface","title":"R Interface","text":"It is worth noting that rcopy will automatically convert data names contains . to _ in Julia. E.g.","category":"page"},{"location":"R_interface/","page":"R Interface","title":"R Interface","text":"julia> rcopy(R\"list(a.b = 1)\")\nOrderedDict{Symbol, Any} with 1 entry:\n :a_b => 1.0","category":"page"},{"location":"R_interface/#Transform-Dta-read-from-R-to-Julia-convention","page":"R Interface","title":"Transform Dta read from R to Julia convention","text":"","category":"section"},{"location":"R_interface/","page":"R Interface","title":"R Interface","text":"If you want to load data using the R interface, but the data source is in the same layout as BUGS, you can process the data in Julia, for instance","category":"page"},{"location":"R_interface/","page":"R Interface","title":"R Interface","text":"# define a row-major reshape function, because Julia's `reshape` is column-major\njulia> function rreshape(v::Vector, dim)\n return permutedims(reshape(v, reverse(dim)), length(dim):-1:1)\n end \nrreshape (generic function with 1 method)\n\njulia> rreshape(vcat(data[:Y]...), (30, 5))\n30×5 Matrix{Float64}:\n 151.0 199.0 246.0 283.0 320.0\n 145.0 199.0 249.0 293.0 354.0\n 147.0 214.0 263.0 312.0 328.0\n 155.0 200.0 237.0 272.0 297.0\n 135.0 188.0 230.0 280.0 323.0\n 159.0 210.0 252.0 298.0 331.0\n 141.0 189.0 231.0 275.0 305.0\n 159.0 201.0 248.0 297.0 338.0\n ⋮ \n 146.0 191.0 229.0 272.0 302.0\n 157.0 211.0 250.0 285.0 323.0\n 132.0 185.0 237.0 286.0 331.0\n 160.0 207.0 257.0 303.0 345.0\n 169.0 216.0 261.0 295.0 333.0\n 157.0 205.0 248.0 289.0 316.0\n 137.0 180.0 219.0 258.0 291.0\n 153.0 200.0 244.0 286.0 324.0","category":"page"},{"location":"R_interface/","page":"R Interface","title":"R Interface","text":"Please always verify the data before using.","category":"page"},{"location":"#Introduction","page":"Introduction","title":"Introduction","text":"","category":"section"},{"location":"","page":"Introduction","title":"Introduction","text":"JuliaBUGS is a graph-based probabilistic programming language and a component of the Turing ecosystem. The package aims to support modelling and inference for probabilistic programs written in the BUGS language. ","category":"page"},{"location":"","page":"Introduction","title":"Introduction","text":"This project is still in its early stage, with many key components needing to be completed. ","category":"page"},{"location":"","page":"Introduction","title":"Introduction","text":"Please refer to the example for usage information and a complete example.","category":"page"},{"location":"#What-is-BUGS?","page":"Introduction","title":"What is BUGS?","text":"","category":"section"},{"location":"","page":"Introduction","title":"Introduction","text":"The BUGS (Bayesian inference Using Gibbs Sampling) system is a probabilistic programming framework designed for specifying directed graphical models. Unlike certain other probabilistic programming languages (PPLs), such as Turing.jl or Pyro, the focus of BUGS is on specifying declarative relationships between nodes in a graph, which can be either logical or stochastic. This means that explicit declarations of variables, inputs, outputs, etc., are not required, and the order of statements is not critical.","category":"page"},{"location":"#The-BUGS-Approach-and-Benefits","page":"Introduction","title":"The BUGS Approach and Benefits","text":"","category":"section"},{"location":"","page":"Introduction","title":"Introduction","text":"Loops in BUGS are essentially a form of \"plate notation,\" offering a concise way to express repetitive statements across many constant indices. Variables in BUGS are either the names of nodes within the program or constant parts of the \"data\" that must be combined with a model for instantiation.","category":"page"},{"location":"","page":"Introduction","title":"Introduction","text":"A BUGS model provides a comprehensive representation of the relationships and dependencies among a set of variables within a Bayesian framework. Our goal is to support BUGS programs as much as possible while also incorporating Julia-specific syntax enhancements.","category":"page"},{"location":"","page":"Introduction","title":"Introduction","text":"The key advantage of utilizing such a graph-based approach is the clarity it provides in understanding the dependencies and relationships within a complex system. These graphical models allow users to explicitly state the conditional dependencies between variables. This makes the model's structure and assumptions transparent, aiding both in the development and interpretation stages. Furthermore, using such a graphical approach makes it easier to apply advanced algorithms for model inference, as it enables more efficient computation by identifying and exploiting the structure of the model.","category":"page"},{"location":"example/#Example:-Logistic-Regression-with-Random-Effects","page":"Example","title":"Example: Logistic Regression with Random Effects","text":"","category":"section"},{"location":"example/","page":"Example","title":"Example","text":"We will use the Seeds model for demonstration. This example concerns the proportion of seeds that germinated on each of 21 plates. Here, we transform the data into a NamedTuple:","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"data = (\n r = [10, 23, 23, 26, 17, 5, 53, 55, 32, 46, 10, 8, 10, 8, 23, 0, 3, 22, 15, 32, 3],\n n = [39, 62, 81, 51, 39, 6, 74, 72, 51, 79, 13, 16, 30, 28, 45, 4, 12, 41, 30, 51, 7],\n x1 = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],\n x2 = [0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1],\n N = 21,\n)","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"where r[i] is the number of germinated seeds and n[i] is the total number of the seeds on the i-th plate. Let p_i be the probability of germination on the i-th plate. Then, the model is defined by:","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"$","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"\\begin{aligned} bi &\\sim \\text{Normal}(0, \\tau) \\\n\\text{logit}(pi) &= \\alpha0 + \\alpha1 x{1 i} + \\alpha2 x{2i} + \\alpha{12} x{1i} x{2i} + b{i} \\\nri &\\sim \\text{Binomial}(pi, ni) \\end{aligned} $","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"where x_1i and x_2i are the seed type and root extract of the i-th plate. The original BUGS program for the model is:","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"model\n{\n for( i in 1 : N ) {\n r[i] ~ dbin(p[i],n[i])\n b[i] ~ dnorm(0.0,tau)\n logit(p[i]) <- alpha0 + alpha1 * x1[i] + alpha2 * x2[i] +\n alpha12 * x1[i] * x2[i] + b[i]\n }\n alpha0 ~ dnorm(0.0, 1.0E-6)\n alpha1 ~ dnorm(0.0, 1.0E-6)\n alpha2 ~ dnorm(0.0, 1.0E-6)\n alpha12 ~ dnorm(0.0, 1.0E-6)\n tau ~ dgamma(0.001, 0.001)\n sigma <- 1 / sqrt(tau)\n}","category":"page"},{"location":"example/#Modeling-Language","page":"Example","title":"Modeling Language","text":"","category":"section"},{"location":"example/#Writing-a-Model-in-BUGS","page":"Example","title":"Writing a Model in BUGS","text":"","category":"section"},{"location":"example/","page":"Example","title":"Example","text":"BUGS language syntax: BNF definition","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"Language References: ","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"MultiBUGS\nOpenBUGS","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"Implementations in C++ and R:","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"JAGS and its user manual\nNimble","category":"page"},{"location":"example/#Writing-a-Model-in-Julia","page":"Example","title":"Writing a Model in Julia","text":"","category":"section"},{"location":"example/","page":"Example","title":"Example","text":"We provide a macro solution which allows users to write down model definitions using Julia:","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"model_def = @bugs begin\n for i in 1:N\n r[i] ~ dbin(p[i], n[i])\n b[i] ~ dnorm(0.0, tau)\n p[i] = logistic(alpha0 + alpha1 * x1[i] + alpha2 * x2[i] + alpha12 * x1[i] * x2[i] + b[i])\n end\n alpha0 ~ dnorm(0.0, 1.0E-6)\n alpha1 ~ dnorm(0.0, 1.0E-6)\n alpha2 ~ dnorm(0.0, 1.0E-6)\n alpha12 ~ dnorm(0.0, 1.0E-6)\n tau ~ dgamma(0.001, 0.001)\n sigma = 1 / sqrt(tau)\nend","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"BUGS syntax carries over almost one-to-one to Julia, with minor exceptions. In general, when basic Julia syntax and BUGS syntax conflict, it is necessary to use Julia syntax. For example, curly braces are replaced with begin ... end blocks, and for loops do not require parentheses. In addition, Julia uses f(x) = ... as a shorthand for function definition, so BUGS' link function syntax can be confusing and ambiguous. Thus, instead of calling the link function, we call the inverse link function from the RHS.","category":"page"},{"location":"example/#Support-for-Legacy-BUGS-Programs","page":"Example","title":"Support for Legacy BUGS Programs","text":"","category":"section"},{"location":"example/","page":"Example","title":"Example","text":"The @bugs macro also works with original (R-like) BUGS syntax:","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"model_def = @bugs(\"\"\"\nmodel{\n for( i in 1 : N ) {\n r[i] ~ dbin(p[i],n[i])\n b[i] ~ dnorm(0.0,tau)\n logit(p[i]) <- alpha0 + alpha1 * x1[i] + alpha2 * x2[i] +\n alpha12 * x1[i] * x2[i] + b[i]\n }\n alpha0 ~ dnorm(0.0,1.0E-6)\n alpha1 ~ dnorm(0.0,1.0E-6)\n alpha2 ~ dnorm(0.0,1.0E-6)\n alpha12 ~ dnorm(0.0,1.0E-6)\n tau ~ dgamma(0.001,0.001)\n sigma <- 1 / sqrt(tau)\n}\n\"\"\", true)","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"By default, @bugs will translate R-style variable names like a.b.c to a_b_c, user can pass false as the second argument to disable this. We still encourage users to write new programs using the Julia-native syntax, because of better debuggability and perks like syntax highlighting. ","category":"page"},{"location":"example/#Compilation","page":"Example","title":"Compilation","text":"","category":"section"},{"location":"example/","page":"Example","title":"Example","text":"For now, the compile function will create a BUGSModel, which implements LogDensityProblems.jl interface.","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"compile(model_def::Expr, data, initializations),","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"The function compile takes three arguments: ","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"the output of @bugs, \nthe data, and\nthe initializations of parameters.","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"initializations = Dict(:alpha => 1, :beta => 1)","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"then we can compile the model with the data and initializations,","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"model = compile(model_def, data, initializations)","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"LogDensityProblemsAD.jl defined some extensions that support automatic differentiation packages. For example, with ReverseDiff.jl","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"using LogDensityProblemsAD, ReverseDiff\n\nad_model = ADgradient(:ReverseDiff, model; compile=Val(true))","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"Here ad_model will also implement all the interfaces of LogDensityProblems.jl. LogDensityProblemsAD.jl will automatically add the interface function logdensity_and_gradient to the model, which will return the log density and gradient of the model. And ad_model can be used in the same way as model in the example below.","category":"page"},{"location":"example/#Inference","page":"Example","title":"Inference","text":"","category":"section"},{"location":"example/","page":"Example","title":"Example","text":"For a differentiable model, we can use AdvancedHMC.jl to perform inference. For instance,","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"using AdvancedHMC, AbstractMCMC, LogDensityProblems, MCMCChains\n\nn_samples, n_adapts = 2000, 1000\n\nD = LogDensityProblems.dimension(model); initial_θ = rand(D)\n\nsamples_and_stats = AbstractMCMC.sample(\n ad_model,\n NUTS(0.8),\n n_samples;\n chain_type = Chains,\n n_adapts = n_adapts,\n init_params = initial_θ,\n discard_initial = n_adapts\n )","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"This will return the MCMC Chain,","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"Chains MCMC chain (2000×40×1 Array{Real, 3}):\n\nIterations = 1001:1:3000\nNumber of chains = 1\nSamples per chain = 2000\nparameters = alpha0, alpha12, alpha1, alpha2, tau, b[16], b[12], b[10], b[14], b[13], b[7], b[6], b[20], b[1], b[4], b[5], b[2], b[18], b[8], b[3], b[9], b[21], b[17], b[15], b[11], b[19], sigma\ninternals = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size, is_adapt\n\nSummary Statistics\n parameters mean std mcse ess_bulk ess_tail rhat ess_per_sec \n Symbol Float64 Float64 Float64 Real Float64 Float64 Missing \n\n alpha0 -0.5642 0.2320 0.0084 766.9305 1022.5211 1.0021 missing\n alpha12 -0.8489 0.5247 0.0170 946.0418 1044.1109 1.0002 missing\n alpha1 0.0587 0.3715 0.0119 966.4367 1233.2257 1.0007 missing\n alpha2 1.3852 0.3410 0.0127 712.2978 974.1566 1.0002 missing\n tau 1.8880 0.7705 0.0447 348.9331 338.3655 1.0030 missing\n b[16] -0.2445 0.4459 0.0132 1528.0578 843.8225 1.0003 missing\n b[12] 0.2050 0.3602 0.0086 1868.6126 1202.1363 0.9996 missing\n b[10] -0.3500 0.2893 0.0090 1047.3119 1245.9358 1.0008 missing\n ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮\n 19 rows omitted\n\nQuantiles\n parameters 2.5% 25.0% 50.0% 75.0% 97.5% \n Symbol Float64 Float64 Float64 Float64 Float64 \n\n alpha0 -1.0143 -0.7143 -0.5590 -0.4100 -0.1185\n alpha12 -1.9063 -1.1812 -0.8296 -0.5153 0.1521\n alpha1 -0.6550 -0.1822 0.0512 0.2885 0.8180\n alpha2 0.7214 1.1663 1.3782 1.5998 2.0986\n tau 0.5461 1.3941 1.8353 2.3115 3.6225\n b[16] -1.2359 -0.4836 -0.1909 0.0345 0.5070\n b[12] -0.4493 -0.0370 0.1910 0.4375 0.9828\n b[10] -0.9570 -0.5264 -0.3331 -0.1514 0.1613\n ⋮ ⋮ ⋮ ⋮ ⋮ ⋮\n 19 rows omitted\n","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"This is consistent with the result in the OpenBUGS seeds example.","category":"page"},{"location":"example/#Parallel-and-Distributed-Sampling-with-AbstractMCMC","page":"Example","title":"Parallel and Distributed Sampling with AbstractMCMC","text":"","category":"section"},{"location":"example/","page":"Example","title":"Example","text":"AbstractMCMC and AdvancedHMC support both parallel and distributed sampling.","category":"page"},{"location":"example/#Parallel-Sampling","page":"Example","title":"Parallel Sampling","text":"","category":"section"},{"location":"example/","page":"Example","title":"Example","text":"To perform multi-threaded sampling of multiple chains, start the Julia session with the -t argument. The model compilation code remains the same, and we can sample multiple chains in parallel as follows:","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"n_chains = 4\nsamples_and_stats = AbstractMCMC.sample(\n ad_model,\n AdvancedHMC.NUTS(0.65),\n AbstractMCMC.MCMCThreads(),\n n_samples,\n n_chains;\n chain_type = Chains,\n n_adapts = n_adapts,\n init_params = [initial_θ for _ = 1:n_chains],\n discard_initial = n_adapts,\n)","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"In this case, we pass two additional arguments to AbstractMCMC.sample:","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"AbstractMCMC.MCMCThreads(): the sampler type, and\nn_chains: the number of chains to sample.","category":"page"},{"location":"example/#Distributed-Sampling","page":"Example","title":"Distributed Sampling","text":"","category":"section"},{"location":"example/","page":"Example","title":"Example","text":"To perform distributed sampling of multiple chains, start the Julia session with the -p argument.","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"In distributed mode, ensure that all functions and modules are available on all processes. Use @everywhere to make the functions and modules available on all processes.","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"For example:","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"@everywhere begin\n using JuliaBUGS, LogDensityProblems, LogDensityProblemsAD, AbstractMCMC, AdvancedHMC, MCMCChains, ReverseDiff # also other packages one may need\n\n # Define the functions to use\n # Use `@register_primitive` to register the functions to use in the model\n\n # Distributed can handle data dependencies in some cases, for more detail, see https://docs.julialang.org/en/v1/manual/distributed-computing/\n\nend\n\nn_chains = nprocs() - 1 # use all the processes except the master process\nsamples_and_stats = AbstractMCMC.sample(\n ad_model,\n AdvancedHMC.NUTS(0.65),\n AbstractMCMC.MCMCDistributed(),\n n_samples,\n n_chains;\n chain_type = Chains,\n n_adapts = n_adapts,\n init_params = [initial_θ for _ = 1:n_chains], # each chain has its own initial parameters\n discard_initial = n_adapts,\n progress = false, # Base.TTY creating problems in distributed setting\n)","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"In this case, we pass two additional arguments to AbstractMCMC.sample:","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"AbstractMCMC.MCMCDistributed(): the sampler type, and\nn_chains: the number of chains to sample.","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"Note that the init_params argument is now a vector of initial parameters for each chain. Sometimes the progress logger can cause problems in distributed setting, so we can disable it by setting progress = false.","category":"page"},{"location":"example/#More-Examples","page":"Example","title":"More Examples","text":"","category":"section"},{"location":"example/","page":"Example","title":"Example","text":"We have transcribed all the examples from the first volume of the BUGS Examples (original and transcribed). All programs and data are included, and can be compiled using the steps described in the tutorial above.","category":"page"}] +[{"location":"pitfalls/#Understanding-Pitfalls-in-Model-Definitions","page":"Pitfalls","title":"Understanding Pitfalls in Model Definitions","text":"","category":"section"},{"location":"pitfalls/#Consequence-of-Observations-on-Model-Parameters","page":"Pitfalls","title":"Consequence of Observations on Model Parameters","text":"","category":"section"},{"location":"pitfalls/","page":"Pitfalls","title":"Pitfalls","text":"When providing observations for the parameters of a model, the dependencies may become disrupted. Consider the following example written in Julia:","category":"page"},{"location":"pitfalls/","page":"Pitfalls","title":"Pitfalls","text":"model_def = @bugs begin\n a ~ Normal(0, 1)\n b ~ Normal(0, 1)\n c ~ Normal(a, b)\nend\n\ndata = (a=1.0, b=2.0)","category":"page"},{"location":"pitfalls/","page":"Pitfalls","title":"Pitfalls","text":"In this scenario, the generated graph will lack the edges a -> c and b -> c, leading the node function of c to become c ~ Normal(1.0, 2.0).","category":"page"},{"location":"pitfalls/#Ambiguity-Between-Data-and-Observed-Stochastic-Variables","page":"Pitfalls","title":"Ambiguity Between Data and Observed Stochastic Variables","text":"","category":"section"},{"location":"pitfalls/","page":"Pitfalls","title":"Pitfalls","text":"A subtle and possibly contentious feature of BUGS syntax is that the observation value of a stochastic variable is treated identically to any model parameters supplied in the data. Here's a legal example in BUGS:","category":"page"},{"location":"pitfalls/","page":"Pitfalls","title":"Pitfalls","text":"model {\n N ~ dcat(p[])\n for (i in 1:N) {\n y[i] ~ dnorm(mu, tau)\n }\n p[1] <- 0.5\n p[2] <- 0.5\n}","category":"page"},{"location":"pitfalls/","page":"Pitfalls","title":"Pitfalls","text":"For a variable to be used as an observation in loop bounds or indexing, it must be part of the provided data, not a transformed variable.","category":"page"},{"location":"pitfalls/","page":"Pitfalls","title":"Pitfalls","text":"This behavior is maintained in the current version of JuliaBUGS, although it was prohibited in the earlier SymbolicPPL.","category":"page"},{"location":"pitfalls/#Possible-Check-Implementation-in-JuliaBUGS","page":"Pitfalls","title":"Possible Check Implementation in JuliaBUGS","text":"","category":"section"},{"location":"pitfalls/","page":"Pitfalls","title":"Pitfalls","text":"Implementing a check for this behavior in JuliaBUGS is feasible. A simplistic approach could be to invalidate (e.g., mark as missing) all observations after the first pass and verify if any are used in loop bounds or indexing. However, there is currently no plan to implement this check.","category":"page"},{"location":"distributions/","page":"Distributions","title":"Distributions","text":"dnorm\ndlogis\ndt\nTDistShiftedScaled\nddexp\ndflat\nFlat\nTruncatedFlat\ndexp\ndchisqr\ndweib\ndlnorm\ndgamma\ndpar\ndgev\ndgpar\ndf\ndunif\ndbeta\ndmnorm\ndmt\ndwish\nddirich\ndbern\ndbin\ndcat\ndpois\ndgeom\ndnegbin\ndbetabin\ndhyper\ndmulti","category":"page"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dnorm","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dnorm","text":"dnorm(μ, τ)\n\nReturns an instance of Normal with mean μ and standard deviation frac1τ. \n\np(xμτ) = sqrtfracτ2π e^-τ frac(x-μ)^22\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dlogis","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dlogis","text":"dlogis(μ, τ)\n\nReturn an instance of Logistic with location parameter μ and scale parameter frac1τ.\n\np(xμτ) = fracsqrtτ e^-sqrtτ(x-μ)(1+e^-sqrtτ(x-μ))^2\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dt","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dt","text":"dt(μ, τ, ν)\n\nIf μ = 0 and σ = 1, the function returns an instance of TDist with ν degrees of freedom, location μ, and scale σ = frac1sqrtτ. Otherwise, it returns an instance of TDistShiftedScaled.\n\np(xνμσ) = fracΓ((ν+1)2)Γ(ν2) sqrtνπσ\nleft(1+frac1νleft(fracx-μσright)^2right)^-fracν+12\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.TDistShiftedScaled","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.TDistShiftedScaled","text":"TDistShiftedScaled(ν, μ, σ)\n\nStudent's t-distribution with ν degrees of freedom, location μ, and scale σ. \n\nThis struct allows for a shift (determined by μ) and a scale (determined by σ) of the standard Student's t-distribution provided by the Distributions.jl package. \n\nOnly pdf and logpdf are implemented for this distribution.\n\nSee Also\n\nTDist\n\n\n\n\n\n","category":"type"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.ddexp","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.ddexp","text":"ddexp(μ, τ)\n\nReturn an instance of Laplace (Double Exponential) with location μ and scale frac1sqrtτ.\n\np(xμτ) = fracsqrtτ2 e^-sqrtτ x-μ\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dflat","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dflat","text":"dflat()\n\nReturns an instance of Flat or TruncatedFlat if truncated.\n\nFlat represents a flat (uniform) prior over the real line, which is an improper distribution. And TruncatedFlat represents a truncated version of the Flat distribution.\n\nOnly pdf, logpdf, minimum, and maximum are implemented for these Distributions.\n\nWhen use in a model, the parameters always need to be initialized.\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.Flat","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.Flat","text":"Flat\n\nThe flat distribution mimicking the behavior of the dflat distribution in the BUGS family of softwares.\n\n\n\n\n\n","category":"type"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.TruncatedFlat","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.TruncatedFlat","text":"TruncatedFlat\n\nTruncated version of the Flat distribution.\n\n\n\n\n\n","category":"type"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dexp","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dexp","text":"dexp(λ)\n\nReturns an instance of Exponential with rate frac1λ.\n\np(xλ) = λ e^-λ x\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dchisqr","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dchisqr","text":"dchisqr(k)\n\nReturns an instance of Chi-squared with k degrees of freedom.\n\np(xk) = frac12^k2 Γ(k2) x^k2 - 1 e^-x2\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dweib","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dweib","text":"dweib(a, b)\n\nReturns an instance of Weibull distribution object with shape parameter a and scale parameter frac1b.\n\nThe Weibull distribution is a common model for event times. The hazard or instantaneous risk of the event is abx^a-1. For a 1 the hazard decreases with x; for a 1 it increases. a = 1 results in the exponential distribution with constant hazard.\n\np(xab) = abx^a-1e^-b x^a\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dlnorm","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dlnorm","text":"dlnorm(μ, τ)\n\nReturns an instance of LogNormal with location μ and scale frac1sqrtτ.\n\np(xμτ) = fracsqrtτxsqrt2π e^-τ2 (log(x) - μ)^2\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dgamma","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dgamma","text":"dgamma(a, b)\n\nReturns an instance of Gamma with shape a and scale frac1b.\n\np(xab) = fracb^aΓ(a) x^a-1 e^-bx\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dpar","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dpar","text":"dpar(a, b)\n\nReturns an instance of Pareto with scale parameter b and shape parameter a.\n\np(xab) = fraca b^ax^a+1\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dgev","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dgev","text":"dgev(μ, σ, η)\n\nReturns an instance of GeneralizedExtremeValue with location μ, scale σ, and shape η.\n\np(xμση) = frac1σ left(1 + η fracx - μσright)^-frac1η - 1 e^-left(1 + η fracx - μσright)^-frac1η\n\nwhere fracη(x - μ)σ -1.\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dgpar","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dgpar","text":"dgpar(μ, σ, η)\n\nReturns an instance of GeneralizedPareto with location μ, scale σ, and shape η.\n\np(xμση) = frac1σ (1 + η ((x - μ)σ))^-1η - 1\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.df","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.df","text":"df(n, m, μ=0, τ=1)\n\nReturns an instance of F-distribution object with n and m degrees of freedom, location μ, and scale τ. This function is only valid when μ = 0 and τ = 1,\n\np(xn m μ τ) = fracGammaleft(fracn+m2right)Gammaleft(fracn2right) Gammaleft(fracm2right) left(fracnmright)^fracn2 sqrtτ left(sqrtτ(x - μ)right)^fracn2-1 left(1 + fracn sqrtτ(x-μ)mright)^-fracn+m2\n\nwhere fracn sqrtτ (x - μ)m -1.\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dunif","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dunif","text":"dunif(a, b)\n\nReturns an instance of Uniform with lower bound a and upper bound b.\n\np(xab) = frac1b - a\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dbeta","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dbeta","text":"dbeta(a, b)\n\nReturns an instance of Beta with shape parameters a and b.\n\np(xab) = fracGamma(a + b)Gamma(a)Gamma(b) x^a-1 (1 - x)^b-1\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dmnorm","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dmnorm","text":"dmnorm(μ::AbstractVector, T::AbstractMatrix)\n\nReturns an instance of Multivariate Normal with mean vector μ and covariance matrix T^-1.\n\np(xμT) = (2π)^-k2 T^12 e^-12 (x-μ) T (x-μ)\n\nwhere k is the dimension of x.\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dmt","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dmt","text":"dmt(μ::AbstractVector, T::AbstractMatrix, k)\n\nReturns an instance of Multivariate T with mean vector μ, scale matrix T^-1, and k degrees of freedom.\n\np(xkμΣ) = fracGamma((k+d)2)Gamma(k2) (kpi)^p2 Σ^12 left(1 + frac1k (x-μ)^T Σ^-1 (x-μ)right)^-frack+p2\n\nwhere p is the dimension of x.\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dwish","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dwish","text":"dwish(R::AbstractMatrix, k)\n\nReturns an instance of Wishart with k degrees of freedom and the scale matrix T^-1.\n\np(XRk) = R^k2 X^(k-p-1)2 e^-(12) tr(RX) (2^kp2 Γ_p(k2))\n\nwhere p is the dimension of X, and it should be less than or equal to k. \n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.ddirich","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.ddirich","text":"ddirich(θ::AbstractVector)\n\nReturn an instance of Dirichlet with parameters θ_i.\n\np(xθ) = fracΓ(sum θ) Γ(θ) x_i^θ_i - 1\n\nwhere theta_i 0 x_i in 0 1 sum_i x_i = 1\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dbern","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dbern","text":"dbern(p)\n\nReturn an instance of Bernoulli with success probability p.\n\np(xp) = p^x (1 - p)^1-x\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dbin","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dbin","text":"dbin(p, n)\n\nReturns an instance of Binomial with number of trials n and success probability p.\n\np(xnp) = binomnx p^x (1 - p)^n-x\n\nend\n\nwhere theta in 0 1 n in mathbbZ^+ and x = 0 ldots n.\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dcat","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dcat","text":"dcat(p)\n\nReturns an instance of Categorical with probabilities p.\n\np(xp) = px\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dpois","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dpois","text":"dpois(θ)\n\nReturns an instance of Poisson with mean (and variance) θ.\n\np(xθ) = e^-θ θ^x x\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dgeom","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dgeom","text":"dgeom(θ)\n\nReturns an instance of Geometric with success probability θ.\n\np(xθ) = (1 - θ)^x-1 θ\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dnegbin","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dnegbin","text":"dnegbin(p, r)\n\nReturns an instance of Negative Binomial with number of failures r and success probability p.\n\nP(xrp) = binomx + r - 1x (1 - p)^x p^r\n\nwhere x in mathbbZ^+.\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dbetabin","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dbetabin","text":"dbetabin(a, b, n)\n\nReturns an instance of Beta Binomial with number of trials n and shape parameters a and b.\n\nP(xa b n) = fracbinomnx binoma + b - 1a + x - 1binoma + b + n - 1n\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dhyper","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dhyper","text":"dhyper(n₁, n₂, m₁, ψ=1)\n\nReturns an instance of Hypergeometric. This distribution is used when sampling without replacement from a population consisting of n₁ successes and n₂ failures, with m₁ being the number of trials or the sample size. The function currently only allows for ψ = 1.\n\np(x n₁ n₂ m₁ psi) = fracbinomn₁x binomn₂m₁ - x psi^xsum_i=u_0^u_1 binomn1i binomn2m₁ - i psi^i\n\nwhere u_0 = max(0 m₁-n₂) u_1 = min(n₁m₁) and u_0 leq x leq u_1\n\n\n\n\n\n","category":"function"},{"location":"distributions/#JuliaBUGS.BUGSPrimitives.dmulti","page":"Distributions","title":"JuliaBUGS.BUGSPrimitives.dmulti","text":"dmulti(θ::AbstractVector, n)\n\nReturns an instance Multinomial with number of trials n and success probabilities θ.\n\nP(xnθ) = fracn_r x_r _r θ_r^x_r\n\n\n\n\n\n","category":"function"},{"location":"BUGS_notes/#Miscellaneous-Notes-on-BUGS","page":"Notes on BUGS Implementations","title":"Miscellaneous Notes on BUGS","text":"","category":"section"},{"location":"BUGS_notes/","page":"Notes on BUGS Implementations","title":"Notes on BUGS Implementations","text":"Here are some exert from BUGS Developer Manual and notes on the original BUGS implementations. ","category":"page"},{"location":"BUGS_notes/#Lexing","page":"Notes on BUGS Implementations","title":"Lexing","text":"","category":"section"},{"location":"BUGS_notes/","page":"Notes on BUGS Implementations","title":"Notes on BUGS Implementations","text":"The BUGS language has the convention that if a name is followed immediately by a round bracket, that is by a \"(\", then the names is a reserved name in the BUGS language and does not represent a variable in the model.\nBy scanning the stream of tokens that constitute a BUGS language model the names of all the variables in the model can be found.","category":"page"},{"location":"BUGS_notes/#Table-of-Names","page":"Notes on BUGS Implementations","title":"Table of Names","text":"","category":"section"},{"location":"BUGS_notes/","page":"Notes on BUGS Implementations","title":"Notes on BUGS Implementations","text":"The BUGS language compiler expands all the for loops in the model and records the value of the indices of each use of a tensor on the left hand side of each relation.\nThe range of each index, for a tensor, is set at the maximum value observed value of the index and added to the name table. There is one exception to this procedure for finding index bounds: names that are data, that is in the data source, have the ranges of their indices fixed in the data source.\nEach scalar and each component of a tensor used on the right hand side of a relation must occur either on the left hand side of a relation and or in a data source.","category":"page"},{"location":"BUGS_notes/#Data-Transformations","page":"Notes on BUGS Implementations","title":"Data Transformations","text":"","category":"section"},{"location":"BUGS_notes/","page":"Notes on BUGS Implementations","title":"Notes on BUGS Implementations","text":"If the compiler can prove that a logical assignment can be evaluated to a constant then the assignment is called a data transformation. This occurs if an assignment's right hand side does not depend on any variable quantities. The BUGS language has a general rule that there must only be one assignment statement for each scalar or component of a tensor. This rule is slightly relaxed for data transformations. The language allows a logical assignment and a stochastic assignment to the same scalar or tensor component if and only if the logical assignment is a data transformation. ","category":"page"},{"location":"BUGS_notes/#Generated-Quantities","page":"Notes on BUGS Implementations","title":"Generated Quantities","text":"","category":"section"},{"location":"BUGS_notes/","page":"Notes on BUGS Implementations","title":"Notes on BUGS Implementations","text":"Only need to be evaluated after the inference algorithm has finished its task. ","category":"page"},{"location":"BUGS_notes/","page":"Notes on BUGS Implementations","title":"Notes on BUGS Implementations","text":"Generally, these are leaf nodes that logical variables\nIn the case of stochastic variables that are leaf nodes, do “forward sampling”, also part of the generated Quantities","category":"page"},{"location":"BUGS_notes/#Computation","page":"Notes on BUGS Implementations","title":"Computation","text":"","category":"section"},{"location":"BUGS_notes/","page":"Notes on BUGS Implementations","title":"Notes on BUGS Implementations","text":"All the nodes in the graphical model representing logical relations are placed into an","category":"page"},{"location":"BUGS_notes/","page":"Notes on BUGS Implementations","title":"Notes on BUGS Implementations","text":"array and sorted by their nesting level with the first array entries only depending on quantities defined by stochastic relations. Traversing this array and evaluating nodes gives up to date values to all logical relations.","category":"page"},{"location":"BUGS_notes/#Types","page":"Notes on BUGS Implementations","title":"Types","text":"","category":"section"},{"location":"BUGS_notes/","page":"Notes on BUGS Implementations","title":"Notes on BUGS Implementations","text":"The BUGS compiler uses the properties of the distribution on the right-hand side of a stochastic assignment statement to make deductions about the variable on the left-hand side. For example, r ~ dbin(p, n) implies that r is integer-valued, while x ~ dnorm(mu, tau) implies that x is real-valued.Some distributions are real-valued but have support on a restricted range of the reals. For example, p ~ dbeta(a, b) implies that p is real-valued with support on the unit interval, while x ~ dgamma(r, lambda) implies that x is real-valued but with support on the positive real line.There are two multivariate distributions in the BUGS language, the Dirichlet and the Wishart, that have support on a complex subspace of the reals. The Dirichlet has support on the unit simplex, while the Wishart has support on symmetric positive definite matrices.The BUGS compiler tries to infer if logical relations return an integer value by looking at whether their parents are integer-valued and the operators that combine the values of their parents into the return value. For example, in the cure model example above, the logical relation state1[i] <- state[i] + 1 is integer-valued because state[i] is a Bernoulli variable and therefore integer, the literal 1 is integer, and the sum of two integers is an integer.When the BUGS system reads in data from a data source, it can tag whether the number read is an integer or a real and propagate this information to logical relations. Again, using the cure model as an example, the statement t[i] <- x[i] + y[i] is integer-valued because both x and y are data and are given as integers in the data source.One special type of data is constants: that is just numbers with no associated distribution. Constants have many uses in BUGS language models, but one of the most important is as covariates. A model can contain a large number of constants that are used as covariates. Because of the possible large numbers of these covariate-type constants, they are given special treatment by the BUGS compiler. If a name read in from a data source is only used on the right-hand side of logical relations, no nodes in the graphical model are created to hold its values; they are directly incorporated in the objects that represent the right-hand sides of the logical relations.For example, the large Methadone model contains the regression:mu.indexed[i] <- beta[1] * x1[i] + beta[2] * x2[i] + beta[3] * x3[i] + beta[4] * x4[i] + region.effect[region.indexed[i]] + source.effect[region.indexed[i]] * source.indexed[i] + person.effect[person.indexed[i]]where i ranges from 1 to 240776. Not having to create a node in the graphical model to represent x1, x2, x3, x4, region.indexed, source.index, and person.indexed saves a large amount of space.In the BUGS language, the type information is fine-grained: each component of a tensor can have different type information. This is quite distinct from the situation in STAN and can make it much easier to specify a statistical model. One common case is where some components of a tensor have been observed while other components need to be estimated. The STAN documentation suggests workarounds for these situations, but these are somewhat complex.","category":"page"},{"location":"BUGS_notes/","page":"Notes on BUGS Implementations","title":"Notes on BUGS Implementations","text":"The type propagation is interesting and maybe useful. But we don’t necessarily need to implement a type system. A dirty way to get type information is simply do a dry run with some tricks.","category":"page"},{"location":"BUGS_notes/#Work-flow","page":"Notes on BUGS Implementations","title":"Work flow","text":"","category":"section"},{"location":"BUGS_notes/","page":"Notes on BUGS Implementations","title":"Notes on BUGS Implementations","text":"The statistical model and data are presented to the BUGS system in a series of stages. In the first stage the model text is parsed into a tree and the name table constructed. The data is then loaded and checked against the model. The data can be split over a number of source. Once all the data has been loaded the model is compiled. Compiling builds the graphical model and does a large number of checks on the consistency of the model. Finally initial values can be given or generated for the model.","category":"page"},{"location":"BUGS_notes/","page":"Notes on BUGS Implementations","title":"Notes on BUGS Implementations","text":"The compiler creates a node in the graphical model for each scalar name and each component of a tensor name in the BUGS language model. The compiler checks that only one node is created for each scalar name or component of a tensor name.","category":"page"},{"location":"BUGS_notes/","page":"Notes on BUGS Implementations","title":"Notes on BUGS Implementations","text":"Reading in a data source causes the compiler to create special nodes called constant nodes to hold the values of the data.","category":"page"},{"location":"BUGS_notes/","page":"Notes on BUGS Implementations","title":"Notes on BUGS Implementations","text":"The compiler processes logical relations before stochastic relations. Any logical relations that only have constant nodes on their right hand side become new constant nodes with the appropriate fixed value. Even if a logical relation can not be reduced to a constant some parts of the relation might be reduced to constants.","category":"page"},{"location":"BUGS_notes/","page":"Notes on BUGS Implementations","title":"Notes on BUGS Implementations","text":"Any constant nodes that have an associated stochastic relation become data nodes in the graphical model.","category":"page"},{"location":"BUGS_notes/#Logical-relations-in-the-BUGS-Language","page":"Notes on BUGS Implementations","title":"Logical relations in the BUGS Language","text":"","category":"section"},{"location":"BUGS_notes/","page":"Notes on BUGS Implementations","title":"Notes on BUGS Implementations","text":"The OpenBUGS software compiles a description of a statistical model in the BUGS language into a graph of objects. Each relation in the statistical model gives rise to a node in the graph of objects. Each distinct type of relation in the statistical model is represented by a node of a distinct class. For stochastic relations there is a fixed set of distributions that can be used in the modelling. For logical relations the situation is more complex. The software can use arbitrary logical expressions build out of a fixed set of basic operators and functions. For each distinct logical expression a new software source code module is written to implement a class to represent that logical expression in the graph of objects. The software module is then compiled using the Components Pascal compiler and the executable code merged into the running OpenBUGS software using the run time loading linker.The BUGS language description of a statistical model is parsed into a list of trees. The sub-trees that represent logical relations in the statistical model are first converted into a stack based representation and then into Component Pascal source code. The source code is generated in module BugsCPWrite and the source code is then compiled in module BugsCPCompiler. Usually the generated source code is not displayed. Checking the Verbose option in the Info menu will cause each each source code module generated by the OpenBUGS software to be displayed in a separate window.One advantage of a stack based representation of an expression is that it is straight forward to use it to derive source code that calculates the derivative of the expression with respect to its arguments. This part of the source code generation is carried out in module BugsCPWrite in procedure WriteEvaluateDiffMethod. Each operator in the stack representation of the logical expression causes a snippet of Component Pascal code to be written. These code snippets are generally very simple with those of binary operators slightly more complex than those of unitary operators. Each binary operators can emit three different code snippets: the general case and two special snippets depending on whether the left or right operands are numerical constants. The only complex code snippet is when an operand that is a logical relation in the statistical model is pushed onto the stack – the > case of nested logical relations. In this case the nested logical relation will have its own code to calculate derivatives and these values can be passed up the nesting level.The OpenBUGS software now uses a backward mode scheme to calculate the value of logical nodes in the statistical model. All the logical nodes in the statistical model are held in a global array and sorted according to their nesting level with unnested nodes at the start of the array. To evaluate all the logical nodes in the statistical model this array is then traversed and each logical node evaluated and the value stored in the node. The same scheme is used to calculate derivatives.The graphs derived from the BUGS language representation of statistical models are generally sparse. The OpenBUGS software uses conditional independence arguments to exploit sparsity in the stochastic parts of the model. There is also a sparsity structure in logical relations.Each logical relation will often depend on just a few stochastic parents and derivatives with respect to other stochastic nodes in the model will be structurally zero. Each logical node has an associated array of stochastic parents for which the derivatives are non zero. Moving up the level of nesting the number of parents can grow. Dealing with this issue leads to the complexity in the code snippet for the operator that pushes a logical node onto the stack. These issues can be seen in the non-linear random effects model called Orange trees in volume II of the OpenBUGS examples. In this model eta[i,] is a function of phi[i,1], phi[i,2] and phi[i,3] where the phi are also logical functions of the stochastic theta[i,].One refinement of the backward mode scheme used to calculate the value of logical nodes is to consider separately any logical nodes in the statistical model which are only used for prediction and do not affect the calculation of the joint probability distribution. These nodes need only be evaluated once per iteration of the inference algorithm. Examples of such nodes are sigma[k] and sigma.C in the Orange trees example. There is no need to evaluate the derivatives of these prediction nodes.The workings of the backward mode scheme are easy to visualize when the inference algorithm updates all the stochastic nodes in the statistical model in one block. Local versions of the backward mode scheme can be used when the inference algorithm works on single nodes or when a small blocks of nodes are updated. Each stochastic node is given its own vector of logical nodes that depend on it either directly or via other logical nodes and this vector is sorted by nesting level. Each updater that works on small blocks of nodes contains a vector of logical nodes which is the union of the vectors of dependent logical nodes for each of its components.The idea of the backward mode scheme for evaluating logical nodes can be used with caching in Metropolis Hastings sampling. First the vector of logical nodes depending on the relevant stochastic node(s) is evaluated and their values cached. The log of the conditional distribution is then calculated. Next a new value of the stochastic node is proposed. The vector of logical nodes is re-evaluated and the log of the > conditional distribution calculated. If the proposed value is rejected then the cache is used to set the vector of logical nodes back to its old values.The OpenBUGS software also calculates what class of function each logical node is in terms of its stochastic parents. If the software can prove for example that a logical node is a linear function of its parents more efficient sampling algorithms can be used. If a linear relation can be proved then the calculation of derivatives can also be optimized in some cases because they will be constant and so only need to be calculated once. Generalized linear models are implemented in a way that allows fast calculation of derivatives. The structure of the algorithm to classify the functional form of logical nodes is very similar to that for derivatives and uses a backward mode scheme","category":"page"},{"location":"BUGS_notes/","page":"Notes on BUGS Implementations","title":"Notes on BUGS Implementations","text":"BUGS separates management of logical and stochastic variables, essentially two graphs. Logical variables are stored in an array and values are updated with values in earlier positions of the array.","category":"page"},{"location":"differences/#Differences-Between-BUGS-and-JuliaBUGS","page":"Differences Between BUGS and JuliaBUGS","title":"Differences Between BUGS and JuliaBUGS","text":"","category":"section"},{"location":"differences/#Implicit-Indexing","page":"Differences Between BUGS and JuliaBUGS","title":"Implicit Indexing","text":"","category":"section"},{"location":"differences/","page":"Differences Between BUGS and JuliaBUGS","title":"Differences Between BUGS and JuliaBUGS","text":"In BUGS, x[, ] is used for implicit indexing, which selects all elements from both the first and second dimensions. In JuliaBUGS, users must explicitly use Colon (:) like x[:, :] when using the @bugs macro. The @bugs macro will insert a Colon when given x[], however, the Julia parser will throw an error if given x[, ]. The original BUGS parser will automatically insert a Colon (:) when it encounters x[, ].","category":"page"},{"location":"functions/","page":"Functions","title":"Functions","text":"Most of the functions from BUGS have been implemented. ","category":"page"},{"location":"functions/","page":"Functions","title":"Functions","text":"JuliaBUGS directly utilizes functions from the Julia Standard Library when they share the same names and functionalities. For functions not available in the Julia Standard Library and other popular libraries, we have developed equivalents within JuliaBUGS.BUGSPrimitives.","category":"page"},{"location":"functions/#Function-defined-in-Julia-Standard-Library","page":"Functions","title":"Function defined in Julia Standard Library","text":"","category":"section"},{"location":"functions/","page":"Functions","title":"Functions","text":"abs(x::Real)\nexp(x::Real)\nlog(x::Number)\nsqrt(x::Real)\ntrunc\nmin(x::Real, y::Real)\nmax(x::Real, y::Real)\nsum(x::AbstractArray)\nsort(x::AbstractArray)\n!!! warning\n `JuliaBUGS` does not allow keyword arguments syntax\nsin(x::Real)\ncos(x::Real)\ntan(x::Real)\nasin(x::Real)\nacos(x::Real)\natan(x::Real)\nasinh(x::Real)\nacosh(x::Real)\natanh(x::Real)\nJuliaBUGS.BUGSPrimitives.mean(x::AbstractArray)","category":"page"},{"location":"functions/#Base.exp-Tuple{Real}","page":"Functions","title":"Base.exp","text":"exp(x)\n\nCompute the natural base exponential of x, in other words ℯ^x.\n\nSee also exp2, exp10 and cis.\n\nExamples\n\njulia> exp(1.0)\n2.718281828459045\n\njulia> exp(im * pi) ≈ cis(pi)\ntrue\n\n\n\n\n\n","category":"method"},{"location":"functions/#Base.log-Tuple{Number}","page":"Functions","title":"Base.log","text":"log(x)\n\nCompute the natural logarithm of x. Throws DomainError for negative Real arguments. Use complex negative arguments to obtain complex results.\n\nSee also ℯ, log1p, log2, log10.\n\nExamples\n\njulia> log(2)\n0.6931471805599453\n\njulia> log(-3)\nERROR: DomainError with -3.0:\nlog was called with a negative real argument but will only return a complex result if called with a complex argument. Try log(Complex(x)).\nStacktrace:\n [1] throw_complex_domainerror(::Symbol, ::Float64) at ./math.jl:31\n[...]\n\njulia> log.(exp.(-1:1))\n3-element Vector{Float64}:\n -1.0\n 0.0\n 1.0\n\n\n\n\n\n","category":"method"},{"location":"functions/#Base.sqrt-Tuple{Real}","page":"Functions","title":"Base.sqrt","text":"sqrt(x)\n\nReturn sqrtx. Throws DomainError for negative Real arguments. Use complex negative arguments instead. The prefix operator √ is equivalent to sqrt.\n\nSee also: hypot.\n\nExamples\n\njulia> sqrt(big(81))\n9.0\n\njulia> sqrt(big(-81))\nERROR: DomainError with -81.0:\nNaN result for non-NaN input.\nStacktrace:\n [1] sqrt(::BigFloat) at ./mpfr.jl:501\n[...]\n\njulia> sqrt(big(complex(-81)))\n0.0 + 9.0im\n\njulia> .√(1:4)\n4-element Vector{Float64}:\n 1.0\n 1.4142135623730951\n 1.7320508075688772\n 2.0\n\n\n\n\n\n","category":"method"},{"location":"functions/#Base.trunc","page":"Functions","title":"Base.trunc","text":"trunc([T,] x)\ntrunc(x; digits::Integer= [, base = 10])\ntrunc(x; sigdigits::Integer= [, base = 10])\n\ntrunc(x) returns the nearest integral value of the same type as x whose absolute value is less than or equal to the absolute value of x.\n\ntrunc(T, x) converts the result to type T, throwing an InexactError if the value is not representable.\n\nKeywords digits, sigdigits and base work as for round.\n\nSee also: %, floor, unsigned, unsafe_trunc.\n\nExamples\n\njulia> trunc(2.22)\n2.0\n\njulia> trunc(-2.22, digits=1)\n-2.2\n\njulia> trunc(Int, -2.22)\n-2\n\n\n\n\n\n","category":"function"},{"location":"functions/#Base.min-Tuple{Real, Real}","page":"Functions","title":"Base.min","text":"min(x, y, ...)\n\nReturn the minimum of the arguments (with respect to isless). See also the minimum function to take the minimum element from a collection.\n\nExamples\n\njulia> min(2, 5, 1)\n1\n\n\n\n\n\n","category":"method"},{"location":"functions/#Base.max-Tuple{Real, Real}","page":"Functions","title":"Base.max","text":"max(x, y, ...)\n\nReturn the maximum of the arguments (with respect to isless). See also the maximum function to take the maximum element from a collection.\n\nExamples\n\njulia> max(2, 5, 1)\n5\n\n\n\n\n\n","category":"method"},{"location":"functions/#Base.sum-Tuple{AbstractArray}","page":"Functions","title":"Base.sum","text":"sum(A::AbstractArray; dims)\n\nSum elements of an array over the given dimensions.\n\nExamples\n\njulia> A = [1 2; 3 4]\n2×2 Matrix{Int64}:\n 1 2\n 3 4\n\njulia> sum(A, dims=1)\n1×2 Matrix{Int64}:\n 4 6\n\njulia> sum(A, dims=2)\n2×1 Matrix{Int64}:\n 3\n 7\n\n\n\n\n\n","category":"method"},{"location":"functions/#Base.sort-Tuple{AbstractArray}","page":"Functions","title":"Base.sort","text":"sort(A; dims::Integer, alg::Algorithm=defalg(A), lt=isless, by=identity, rev::Bool=false, order::Ordering=Forward)\n\nSort a multidimensional array A along the given dimension. See sort! for a description of possible keyword arguments.\n\nTo sort slices of an array, refer to sortslices.\n\nExamples\n\njulia> A = [4 3; 1 2]\n2×2 Matrix{Int64}:\n 4 3\n 1 2\n\njulia> sort(A, dims = 1)\n2×2 Matrix{Int64}:\n 1 2\n 4 3\n\njulia> sort(A, dims = 2)\n2×2 Matrix{Int64}:\n 3 4\n 1 2\n\n\n\n\n\n","category":"method"},{"location":"functions/#Function-defined-in-[LogExpFunctions](https://github.com/JuliaStats/LogExpFunctions.jl)","page":"Functions","title":"Function defined in LogExpFunctions","text":"","category":"section"},{"location":"functions/","page":"Functions","title":"Functions","text":"cloglog\ncexpexp\nlogit\nlogistic","category":"page"},{"location":"functions/#LogExpFunctions.cloglog","page":"Functions","title":"LogExpFunctions.cloglog","text":"cloglog(x)\n\n\nCompute the complementary log-log, log(-log(1 - x)).\n\n\n\n\n\n","category":"function"},{"location":"functions/#LogExpFunctions.cexpexp","page":"Functions","title":"LogExpFunctions.cexpexp","text":"cexpexp(x)\n\n\nCompute the complementary double exponential, 1 - exp(-exp(x)).\n\n\n\n\n\n","category":"function"},{"location":"functions/#LogExpFunctions.logit","page":"Functions","title":"LogExpFunctions.logit","text":"logit(x)\n\n\nThe logit or log-odds transformation, defined as\n\noperatornamelogit(x) = logleft(fracx1-xright)\n\nfor 0 x 1.\n\nIts inverse is the logistic function.\n\n\n\n\n\n","category":"function"},{"location":"functions/#LogExpFunctions.logistic","page":"Functions","title":"LogExpFunctions.logistic","text":"logistic(x)\n\n\nThe logistic sigmoid function mapping a real number to a value in the interval 01,\n\nsigma(x) = frac1e^-x + 1 = frace^x1+e^x\n\nIts inverse is the logit function.\n\n\n\n\n\n","category":"function"},{"location":"functions/#Function-defined-in-JuliaBUGS.BUGSPrimitives","page":"Functions","title":"Function defined in JuliaBUGS.BUGSPrimitives","text":"","category":"section"},{"location":"functions/","page":"Functions","title":"Functions","text":"JuliaBUGS.BUGSPrimitives.equals\nJuliaBUGS.BUGSPrimitives.inprod\nJuliaBUGS.BUGSPrimitives.inverse\nJuliaBUGS.BUGSPrimitives.logdet\nJuliaBUGS.BUGSPrimitives.logfact\nJuliaBUGS.BUGSPrimitives.loggam\nJuliaBUGS.BUGSPrimitives.icloglog\nJuliaBUGS.BUGSPrimitives.mexp\nJuliaBUGS.BUGSPrimitives.phi\nJuliaBUGS.BUGSPrimitives.pow\nJuliaBUGS.BUGSPrimitives.rank\nJuliaBUGS.BUGSPrimitives.ranked\nJuliaBUGS.BUGSPrimitives.sd\nJuliaBUGS.BUGSPrimitives.softplus\nJuliaBUGS.BUGSPrimitives._step\nJuliaBUGS.BUGSPrimitives.arcsin\nJuliaBUGS.BUGSPrimitives.arcsinh\nJuliaBUGS.BUGSPrimitives.arccos\nJuliaBUGS.BUGSPrimitives.arccosh\nJuliaBUGS.BUGSPrimitives.arctan\nJuliaBUGS.BUGSPrimitives.arctanh","category":"page"},{"location":"functions/#JuliaBUGS.BUGSPrimitives.equals","page":"Functions","title":"JuliaBUGS.BUGSPrimitives.equals","text":"equals(x, y)\n\nReturns 1 if x is equal to y, 0 otherwise.\n\n\n\n\n\n","category":"function"},{"location":"functions/#JuliaBUGS.BUGSPrimitives.inprod","page":"Functions","title":"JuliaBUGS.BUGSPrimitives.inprod","text":"inprod(a, b)\n\nInner product of a and b.\n\n\n\n\n\n","category":"function"},{"location":"functions/#JuliaBUGS.BUGSPrimitives.inverse","page":"Functions","title":"JuliaBUGS.BUGSPrimitives.inverse","text":"inverse(m::AbstractMatrix)\n\nInverse of matrix mathbfm.\n\n\n\n\n\n","category":"function"},{"location":"functions/#JuliaBUGS.BUGSPrimitives.logdet","page":"Functions","title":"JuliaBUGS.BUGSPrimitives.logdet","text":"logdet(::AbstractMatrix)\n\nLogarithm of the determinant of matrix mathbfv.\n\n\n\n\n\n","category":"function"},{"location":"functions/#JuliaBUGS.BUGSPrimitives.logfact","page":"Functions","title":"JuliaBUGS.BUGSPrimitives.logfact","text":"logfact(x)\n\nLogarithm of the factorial of x.\n\n\n\n\n\n","category":"function"},{"location":"functions/#JuliaBUGS.BUGSPrimitives.loggam","page":"Functions","title":"JuliaBUGS.BUGSPrimitives.loggam","text":"loggam(x)\n\nLogarithm of the gamma function of x.\n\n\n\n\n\n","category":"function"},{"location":"functions/#JuliaBUGS.BUGSPrimitives.icloglog","page":"Functions","title":"JuliaBUGS.BUGSPrimitives.icloglog","text":"icloglog(x)\n\nInverse complementary log-log function of x. Alias for cexpexp(x).\n\n\n\n\n\n","category":"function"},{"location":"functions/#JuliaBUGS.BUGSPrimitives.mexp","page":"Functions","title":"JuliaBUGS.BUGSPrimitives.mexp","text":"mexp(x::AbstractMatrix)\n\nMatrix exponential of mathbfx.\n\n\n\n\n\n","category":"function"},{"location":"functions/#JuliaBUGS.BUGSPrimitives.phi","page":"Functions","title":"JuliaBUGS.BUGSPrimitives.phi","text":"phi(x)\n\nCumulative distribution function (CDF) of the standard normal distribution evaluated at x.\n\n\n\n\n\n","category":"function"},{"location":"functions/#JuliaBUGS.BUGSPrimitives.pow","page":"Functions","title":"JuliaBUGS.BUGSPrimitives.pow","text":"pow(a, b)\n\nReturn a raised to the power of b.\n\n\n\n\n\n","category":"function"},{"location":"functions/#JuliaBUGS.BUGSPrimitives.rank","page":"Functions","title":"JuliaBUGS.BUGSPrimitives.rank","text":"rank(v::AbstractVector, i::Integer)\n\nReturn the rank of the i-th element of mathbfv.\n\n\n\n\n\n","category":"function"},{"location":"functions/#JuliaBUGS.BUGSPrimitives.ranked","page":"Functions","title":"JuliaBUGS.BUGSPrimitives.ranked","text":"ranked(v::AbstractVector, i::Integer)\n\nReturn the i-th element of mathbfv sorted in ascending order.\n\n\n\n\n\n","category":"function"},{"location":"functions/#JuliaBUGS.BUGSPrimitives.sd","page":"Functions","title":"JuliaBUGS.BUGSPrimitives.sd","text":"sd(v::AbstractVector)\n\nReturn the standard deviation of the input vector mathbfv.\n\n\n\n\n\n","category":"function"},{"location":"functions/#JuliaBUGS.BUGSPrimitives.softplus","page":"Functions","title":"JuliaBUGS.BUGSPrimitives.softplus","text":"softplus(x)\n\nReturn the softplus function of x, defined as log(1 + exp(x)).\n\n\n\n\n\n","category":"function"},{"location":"functions/#JuliaBUGS.BUGSPrimitives._step","page":"Functions","title":"JuliaBUGS.BUGSPrimitives._step","text":"_step(x)\n\nReturn 1 if x is greater than 0, and 0 otherwise.\n\n\n\n\n\n","category":"function"},{"location":"functions/#JuliaBUGS.BUGSPrimitives.arcsin","page":"Functions","title":"JuliaBUGS.BUGSPrimitives.arcsin","text":"arcsin(x)\n\nReturn the arcsine of x.\n\n\n\n\n\n","category":"function"},{"location":"functions/#JuliaBUGS.BUGSPrimitives.arcsinh","page":"Functions","title":"JuliaBUGS.BUGSPrimitives.arcsinh","text":"arcsinh(x)\n\nReturn the inverse hyperbolic sine of x.\n\n\n\n\n\n","category":"function"},{"location":"functions/#JuliaBUGS.BUGSPrimitives.arccos","page":"Functions","title":"JuliaBUGS.BUGSPrimitives.arccos","text":"arccos(x)\n\nReturn the arccosine of x.\n\n\n\n\n\n","category":"function"},{"location":"functions/#JuliaBUGS.BUGSPrimitives.arccosh","page":"Functions","title":"JuliaBUGS.BUGSPrimitives.arccosh","text":"arccosh(x)\n\nReturn the inverse hyperbolic cosine of x.\n\n\n\n\n\n","category":"function"},{"location":"functions/#JuliaBUGS.BUGSPrimitives.arctan","page":"Functions","title":"JuliaBUGS.BUGSPrimitives.arctan","text":"arctan(x)\n\nReturn the arctangent of x.\n\n\n\n\n\n","category":"function"},{"location":"functions/#JuliaBUGS.BUGSPrimitives.arctanh","page":"Functions","title":"JuliaBUGS.BUGSPrimitives.arctanh","text":"arctanh(x)\n\nSee atanh.\n\n\n\n\n\n","category":"function"},{"location":"user_defined_functions/#Define-and-Use-Your-Own-Functions-and-Distributions","page":"User-defined Functions and Distributions","title":"Define and Use Your Own Functions and Distributions","text":"","category":"section"},{"location":"user_defined_functions/","page":"User-defined Functions and Distributions","title":"User-defined Functions and Distributions","text":"For now, out of the box, JuliaBUGS only allows functions and distributions defined in BUGSPrimitives to be used in the model. With the @register_primitive macro, users can register their own functions and distributions with JuliaBUGS. It is important to ensure that any functions used are pure mathematical functions. This implies that such functions should not alter any external state including but not limited to modifying global variables, writing data to files. (Printing might be okay, but do at discretion.)","category":"page"},{"location":"user_defined_functions/","page":"User-defined Functions and Distributions","title":"User-defined Functions and Distributions","text":"julia> JuliaBUGS.@register_primitive function f(x)\n return x + 1\nend\nf (generic function with 1 method)\n\njulia> JuliaBUGS.f(2)\n3","category":"page"},{"location":"user_defined_functions/","page":"User-defined Functions and Distributions","title":"User-defined Functions and Distributions","text":"Users can also introduce a function into JuliaBUGS, by ","category":"page"},{"location":"user_defined_functions/","page":"User-defined Functions and Distributions","title":"User-defined Functions and Distributions","text":"julia> f(x) = x + 1\nf (generic function with 1 method)\n\njulia> JuliaBUGS.@register_primitive(f);\n\njulia> JuliaBUGS.f(1)\n2","category":"page"},{"location":"user_defined_functions/","page":"User-defined Functions and Distributions","title":"User-defined Functions and Distributions","text":"After registering the function or distributions, they can be used just like any other functions or distributions provided by BUGS.","category":"page"},{"location":"parser/#BUGS-Parser","page":"Parser","title":"BUGS Parser","text":"","category":"section"},{"location":"parser/","page":"Parser","title":"Parser","text":"The macro @bugs produces a Julia Expr object that represents the BUGS model definition.","category":"page"},{"location":"parser/","page":"Parser","title":"Parser","text":"If the input is a String, it's assumed to be a program in the original BUGS language. In this case, the macro will first convert the program to an equivalent Julia program, then use the Julia parser to parse the program into an Expr object.","category":"page"},{"location":"parser/","page":"Parser","title":"Parser","text":"Both model definitions written in Julia and those written in the original BUGS and subsequently parsed are now represented as a Julia Expr object. These objects go through syntax checking and post-processing to create the input for the compile function.","category":"page"},{"location":"parser/","page":"Parser","title":"Parser","text":"Below, we describe how the original BUGS program is translated to an equivalent Julia program and detail the post-processing done to the Expr object.","category":"page"},{"location":"parser/#BUGS-to-Julia-Translation","page":"Parser","title":"BUGS to Julia Translation","text":"","category":"section"},{"location":"parser/","page":"Parser","title":"Parser","text":"In this section, we refer to the translation program as the \"parser\" and the translating process as \"parsing\". Although the parser doesn't produce a syntax tree, it does follow the form of a recursive descent parser, building a Julia program in the form of a vector of tokens rather than a syntax tree.","category":"page"},{"location":"parser/","page":"Parser","title":"Parser","text":"This general implementation is heavily inspired by JuliaSyntax.jl, the official parser for Julia since version 1.10.","category":"page"},{"location":"parser/","page":"Parser","title":"Parser","text":"The BUGS parser implemented here takes a token stream with a recursive descent structure and checks the program's correctness. Here's how it works:","category":"page"},{"location":"parser/","page":"Parser","title":"Parser","text":"Use tokenize to obtain the token vector.\nInspect the tokens and build the Julia version of the program as a vector of tokens.\nPush the token to the Julia version of the program vector when appropriate.\nDetect errors and make necessary alterations to tokens, such as deletion, combination, or replacement.","category":"page"},{"location":"parser/","page":"Parser","title":"Parser","text":"During the recursive descent, BUGS syntax tokens will be translated into Julia syntax tokens. Some tokens will remain as they are, while others will be transformed, removed, or new tokens may be added.","category":"page"},{"location":"parser/","page":"Parser","title":"Parser","text":"The parser will throw an error if it encounters a program that does not adhere to strict BUGS syntax.","category":"page"},{"location":"parser/#Some-Notes-on-Error-Recovery","page":"Parser","title":"Some Notes on Error Recovery","text":"","category":"section"},{"location":"parser/","page":"Parser","title":"Parser","text":"The current error recovery is ad hoc and primarily rudimentary. If the program is correct, it will produce the correct result. If the program is syntactically or semantically incorrect, the token stream will not be pushed forward, resulting in failure.","category":"page"},{"location":"parser/","page":"Parser","title":"Parser","text":"The failure detection mechanism checks if two errors occur with the same \"current token\". If they do, the parser stops and reports the error. This ensures that the parser won't incorrectly parse a flawed program.","category":"page"},{"location":"bugs_lang/#Special-Cases-in-the-BUGS-Language","page":"Special Cases in the BUGS Language","title":"Special Cases in the BUGS Language","text":"","category":"section"},{"location":"bugs_lang/","page":"Special Cases in the BUGS Language","title":"Special Cases in the BUGS Language","text":"Here we record some of the special cases in the BUGS language in the original BUGS softwares, JuliaBUGS may or may not inherent these behaviors.","category":"page"},{"location":"bugs_lang/#Nested-Indexing-on-the-Left-Hand-Side","page":"Special Cases in the BUGS Language","title":"Nested Indexing on the Left-Hand Side","text":"","category":"section"},{"location":"bugs_lang/","page":"Special Cases in the BUGS Language","title":"Special Cases in the BUGS Language","text":"model_def = @bugs x[y[1]] ~ dnorm(0, 1) # this is permitted\ndata = (y=[1,2,3],)","category":"page"},{"location":"bugs_lang/#Function-Application-on-the-Left-Hand-Side-Index-is-Not-Permitted","page":"Special Cases in the BUGS Language","title":"Function Application on the Left-Hand Side Index is Not Permitted","text":"","category":"section"},{"location":"bugs_lang/","page":"Special Cases in the BUGS Language","title":"Special Cases in the BUGS Language","text":"This is identified during the compilation phase, not during parsing.","category":"page"},{"location":"bugs_lang/#The-arguments-to-the-distributions-functions-must-be-\"simple\":-they-are-either-a-variable-or-a-constant,-their-should-be-no-function-application-other-than-the-arithmetic-operators","page":"Special Cases in the BUGS Language","title":"The arguments to the distributions functions must be \"simple\": they are either a variable or a constant, their should be no function application other than the arithmetic operators","text":"","category":"section"},{"location":"bugs_lang/","page":"Special Cases in the BUGS Language","title":"Special Cases in the BUGS Language","text":"model {\n x ~ dnorm(y[1] + 1, 1)\n} # fail at parsing\n\nmodel {\n x ~ dnorm(sum(y[1:2]), 1)\n} # fail at parsing\n\nmodel {\n x ~ dnorm(y[sum(y[1:2])], 1)\n} # pass the parser, but fail at compile\n\nmodel {\n x ~ dnorm(y[y[2]], 1)\n} # this is okay\n\nmodel {\n x ~ dnorm(y[y[2]+1], 1)\n} # this is also okay\n\nlist(y = c(1, 2, 3))","category":"page"},{"location":"bugs_lang/","page":"Special Cases in the BUGS Language","title":"Special Cases in the BUGS Language","text":"The reason behind this is two fold: (1) forcing the user to use ~ to identify the dependency between random variables, instead of expressions, (2) make implementation of automatic differentiation easier.","category":"page"},{"location":"bugs_lang/#Indexing-with-Transformed-Variables-is-Not-Permitted","page":"Special Cases in the BUGS Language","title":"Indexing with Transformed Variables is Not Permitted","text":"","category":"section"},{"location":"bugs_lang/","page":"Special Cases in the BUGS Language","title":"Special Cases in the BUGS Language","text":"This restriction is due to the current compiler implementation.","category":"page"},{"location":"bugs_lang/","page":"Special Cases in the BUGS Language","title":"Special Cases in the BUGS Language","text":"model{\n a <- max(y[2], y[3])\n x[a + 1] ~ dnorm(0, 1)\n}\n\nlist(y=c(1, 2, 3))","category":"page"},{"location":"bugs_lang/","page":"Special Cases in the BUGS Language","title":"Special Cases in the BUGS Language","text":"The reason for this is that the transformed variable and indices are evaluated simultaneously, and the transformed variable is not available at the time of indexing. This is a limitation of the current implementation.","category":"page"},{"location":"api/#API","page":"General","title":"API","text":"","category":"section"},{"location":"api/","page":"General","title":"General","text":"@bugs\ncompile\nBUGSModel\nConcreteNodeInfo\nBUGSGraph","category":"page"},{"location":"api/#JuliaBUGS.Parser.@bugs","page":"General","title":"JuliaBUGS.Parser.@bugs","text":"@bugs(prog::String, replace_period=true, no_enclosure=false)\n\nProduce similar output as @bugs, but takes a string as input. This is useful for parsing original BUGS programs.\n\nArguments\n\nprog::String: The BUGS program code as a string.\nreplace_period::Bool: If true, periods in the BUGS code will be replaced (default true).\nno_enclosure::Bool: If true, the parser will not expect the program to be wrapped between model{ } (default false).\n\n\n\n\n\n","category":"macro"},{"location":"api/#JuliaBUGS.compile","page":"General","title":"JuliaBUGS.compile","text":"compile(model_def[, data, initializations])\n\nCompile a BUGS model into a log density problem.\n\nArguments\n\nmodel_def::Expr: The BUGS model definition.\ndata::NamedTuple or AbstractDict: The data to be used in the model. If none is passed, the data will be assumed to be empty.\ninitializations::NamedTuple or AbstractDict: The initial values for the model parameters. If none is passed, the parameters will be assumed to be initialized to zero.\nis_transformed::Bool=true: If true, the model parameters during inference will be transformed to the unconstrained space. \n\nReturns\n\nA BUGSModel object representing the compiled model.\n\n\n\n\n\n","category":"function"},{"location":"api/#JuliaBUGS.BUGSModel","page":"General","title":"JuliaBUGS.BUGSModel","text":"BUGSModel\n\nThe BUGSModel object is used for inference and represents the output of compilation. It fully implements the LogDensityProblems.jl interface.\n\nFields\n\ntransformed::Bool: Indicates whether the model parameters are in the transformed space.\nuntransformed_param_length::Int: The length of the parameters vector in the original space.\ntransformed_param_length::Int: The length of the parameters vector in the transformed space.\nuntransformed_var_lengths::Dict{VarName,Int}: A dictionary mapping the names of the variables to their lengths in the original space.\ntransformed_var_lengths::Dict{VarName,Int}: A dictionary mapping the names of the variables to their lengths in the transformed space.\nvarinfo::SimpleVarInfo: An instance of DynamicPPL.SimpleVarInfo, which is a dictionary-like data structure that maps both data and values of variables in the model to the corresponding values.\nparameters::Vector{VarName}: A vector containing the names of the parameters in the model, defined as stochastic variables that are not observed. This vector should be consistent with sorted_nodes.\nsorted_nodes::Vector{VarName}: A vector containing the names of all the variables in the model, sorted in topological order. In the case of a conditioned model, sorted_nodes include all the variables in parameters and the variables in the Markov blanket of parameters.\ng::BUGSGraph: An instance of BUGSGraph, representing the dependency graph of the model.\nbase_model::Union{BUGSModel,Nothing}: If not Nothing, the model is a conditioned model; otherwise, it's the model returned by compile.\n\n\n\n\n\n","category":"type"},{"location":"api/#JuliaBUGS.ConcreteNodeInfo","page":"General","title":"JuliaBUGS.ConcreteNodeInfo","text":"ConcreteNodeInfo\n\nDefines the information stored in each node of the BUGS graph, encapsulating the essential characteristics and functions associated with a node within the BUGS model's dependency graph.\n\nFields\n\nnode_type::VariableTypes: Specifies whether the node is a stochastic or logical variable.\nnode_function_expr::Expr: The node function expression.\nnode_args::Vector{VarName}: A vector containing the names of the variables that are arguments to the node function.\n\n\n\n\n\n","category":"type"},{"location":"api/#JuliaBUGS.BUGSGraph","page":"General","title":"JuliaBUGS.BUGSGraph","text":"BUGSGraph\n\nThe BUGSGraph object represents the graph structure for a BUGS model. It is a type alias for MetaGraphsNext.MetaGraph with node type specified to ConcreteNodeInfo.\n\n\n\n\n\n","category":"type"},{"location":"graph_plotting/#Plotting-graphs","page":"Plotting","title":"Plotting graphs","text":"","category":"section"},{"location":"graph_plotting/","page":"Plotting","title":"Plotting","text":"Plotting the graph can be very beneficial for debugging the model.","category":"page"},{"location":"graph_plotting/","page":"Plotting","title":"Plotting","text":"Note Plate notation is not yet supported. Therefore, it's advisable for users to begin with a more streamlined model that contains fewer nodes, allowing for clearer visualization.","category":"page"},{"location":"graph_plotting/","page":"Plotting","title":"Plotting","text":"In Julia, we've set up standard plotting routines using various graphing libraries. You can visualize graphs with three different libraries by employing a common model, as detailed below:","category":"page"},{"location":"graph_plotting/","page":"Plotting","title":"Plotting","text":"model_def = @bugs begin\n a ~ dnorm(f, c)\n f = b - 1\n b ~ dnorm(0, 1)\n c ~ dnorm(l, 1)\n g = a * 2\n d ~ dnorm(g, 1)\n h = g + 2\n e ~ dnorm(h, i)\n i ~ dnorm(0, 1)\n l ~ dnorm(0, 1)\nend\n\ninits = Dict(\n :a => 1.0,\n :b => 2.0,\n :c => 3.0,\n :d => 4.0,\n :e => 5.0,\n\n # :f => 1.0,\n # :g => 2.0,\n # :h => 4.0,\n\n :i => 4.0,\n :l => -2.0,\n)\n\nmodel = compile(model_def, NamedTuple(), inits)","category":"page"},{"location":"graph_plotting/#[TikzGraphs.jl](https://github.com/JuliaTeX/TikzGraphs.jl).","page":"Plotting","title":"TikzGraphs.jl.","text":"","category":"section"},{"location":"graph_plotting/","page":"Plotting","title":"Plotting","text":"using TikzGraphs\nTikzGraphs.plot(model)","category":"page"},{"location":"graph_plotting/","page":"Plotting","title":"Plotting","text":"(Image: TikzGraphs)","category":"page"},{"location":"graph_plotting/#[GraphPlot.jl](https://github.com/JuliaGraphs/GraphPlot.jl)","page":"Plotting","title":"GraphPlot.jl","text":"","category":"section"},{"location":"graph_plotting/","page":"Plotting","title":"Plotting","text":"using GraphPlot\ngplot(model)","category":"page"},{"location":"graph_plotting/","page":"Plotting","title":"Plotting","text":"(Image: GraphPlot)","category":"page"},{"location":"graph_plotting/#[GraphMakie.jl](https://github.com/MakieOrg/GraphMakie.jl)","page":"Plotting","title":"GraphMakie.jl","text":"","category":"section"},{"location":"graph_plotting/","page":"Plotting","title":"Plotting","text":"using GLMakie, GraphMakie\ngraphplot(model)","category":"page"},{"location":"graph_plotting/","page":"Plotting","title":"Plotting","text":"(Image: GraphMakie)","category":"page"},{"location":"bugs_examples/Leuk/#Leuk:-Cox-regression","page":"Leuk: Cox regression","title":"Leuk: Cox regression","text":"","category":"section"},{"location":"bugs_examples/Leuk/#Description","page":"Leuk: Cox regression","title":"Description","text":"","category":"section"},{"location":"bugs_examples/Leuk/","page":"Leuk: Cox regression","title":"Leuk: Cox regression","text":"Several authors have discussed Bayesian inference for censored survival data where the integrated baseline hazard function is to be estimated non-parametrically, including Kalbfleisch (1978), Kalbfleisch and Prentice (1980), Clayton (1991), and Clayton (1994). Clayton (1994) formulates the Cox model using counting process notation introduced by Andersen and Gill (1982) and discusses estimation of the baseline hazard and regression parameters using MCMC methods. Although this approach may seem somewhat contrived, it lays the groundwork for extensions to random effect (frailty) models, time-dependent covariates, smoothed hazards, multiple events, and more. Below is how to implement this formulation of the Cox model in BUGS.","category":"page"},{"location":"bugs_examples/Leuk/","page":"Leuk: Cox regression","title":"Leuk: Cox regression","text":"For subjects i = 1n, we observe processes N_i(t) which count the number of failures up to time t. The corresponding intensity process I_i(t) is given by","category":"page"},{"location":"bugs_examples/Leuk/","page":"Leuk: Cox regression","title":"Leuk: Cox regression","text":"I_i(t) dt = mathbbEdN_i(t) mathcalF_t-","category":"page"},{"location":"bugs_examples/Leuk/","page":"Leuk: Cox regression","title":"Leuk: Cox regression","text":"where dN_i(t) is the increment of N_i over the interval t t+dt), and mathcalF_t- represents the available data just before time t. If subject i fails during this interval, dN_i(t) = 1; otherwise, dN_i(t) = 0. Hence, mathbbE(dN_i(t) mathcalF_t-) corresponds to the probability of subject i failing in the interval t t+dt). As dt to 0, this probability becomes the instantaneous hazard at time t for subject i, assumed to have the form","category":"page"},{"location":"bugs_examples/Leuk/","page":"Leuk: Cox regression","title":"Leuk: Cox regression","text":"I_i(t) = Y_i(t)lambda_0(t) exp(beta z_i)","category":"page"},{"location":"bugs_examples/Leuk/","page":"Leuk: Cox regression","title":"Leuk: Cox regression","text":"where Y_i(t) is an observed process taking the value 1 or 0 according to whether subject i is at risk at time t, and lambda_0(t) exp(beta z_i) is the Cox regression model. Thus, the observed data D = N_i(t) Y_i(t) z_i i = 1n and unknown parameters beta and Lambda_0(t) = int_0^t lambda_0(u) du, the latter estimated non-parametrically.","category":"page"},{"location":"bugs_examples/Leuk/","page":"Leuk: Cox regression","title":"Leuk: Cox regression","text":"The joint posterior distribution is","category":"page"},{"location":"bugs_examples/Leuk/","page":"Leuk: Cox regression","title":"Leuk: Cox regression","text":"P(beta Lambda_0() D) propto P(D beta Lambda_0()) P(beta) P(Lambda_0())","category":"page"},{"location":"bugs_examples/Leuk/","page":"Leuk: Cox regression","title":"Leuk: Cox regression","text":"For BUGS, specify the likelihood P(D beta Lambda_0()) and priors for beta and Lambda_0(). Under non-informative censoring, the data likelihood is","category":"page"},{"location":"bugs_examples/Leuk/","page":"Leuk: Cox regression","title":"Leuk: Cox regression","text":"prod_i=1^n left( prod_t geq 0 I_i(t) dN_i(t) right) exp(- I_i(t) dt)","category":"page"},{"location":"bugs_examples/Leuk/","page":"Leuk: Cox regression","title":"Leuk: Cox regression","text":"Viewing the increments dN_i(t) as independent Poisson variables with means I_i(t)dt:","category":"page"},{"location":"bugs_examples/Leuk/","page":"Leuk: Cox regression","title":"Leuk: Cox regression","text":"dN_i(t) sim textPoisson(I_i(t)dt)","category":"page"},{"location":"bugs_examples/Leuk/","page":"Leuk: Cox regression","title":"Leuk: Cox regression","text":"I_i(t)dt = Y_i(t) exp(beta z_i) dLambda_0(t)","category":"page"},{"location":"bugs_examples/Leuk/","page":"Leuk: Cox regression","title":"Leuk: Cox regression","text":"where (d\\Lambda0(t) = \\Lambda0(t)dt) is the increment or jump in the integrated baseline hazard function occurring during the time interval ([t, t+dt)). Since the conjugate prior for the Poisson mean is the gamma distribution, it would be convenient if (\\Lambda0()) were a process in which the increments (d\\Lambda0(t)) are distributed according to gamma distributions. We assume the conjugate independent increments prior suggested by Kalbfleisch (1978), namely ","category":"page"},{"location":"bugs_examples/Leuk/","page":"Leuk: Cox regression","title":"Leuk: Cox regression","text":"dLambda_0(t) sim textGamma(c cdot dLambda^*_0(t) c)","category":"page"},{"location":"bugs_examples/Leuk/","page":"Leuk: Cox regression","title":"Leuk: Cox regression","text":"Here, dLambda^*_0(t) can be thought of as a prior guess at the unknown hazard function, with c representing the degree of confidence in this guess. Small values of c correspond to weak prior beliefs. In the example below, we set dLambda^*_0(t) = r cdot dt where r is a guess at the failure rate per unit time, and dt is the size of the time interval. ","category":"page"},{"location":"bugs_examples/Leuk/","page":"Leuk: Cox regression","title":"Leuk: Cox regression","text":"The above formulation is appropriate when genuine prior information exists concerning the underlying hazard function. Alternatively, if we wish to reproduce a Cox analysis but with, say, additional hierarchical structure, we may use the Multinomial-Poisson trick described in the BUGS manual. This is equivalent to assuming independent increments in the cumulative non-informative priors. This formulation is also shown below.","category":"page"},{"location":"bugs_examples/Leuk/","page":"Leuk: Cox regression","title":"Leuk: Cox regression","text":"The fixed effect regression coefficients (b) are assigned a vague prior","category":"page"},{"location":"bugs_examples/Leuk/","page":"Leuk: Cox regression","title":"Leuk: Cox regression","text":"b sim textNormal(00 0000001)","category":"page"},{"location":"bugs_examples/Leuk/#BUGS-code-for-the-Leuk-example:","page":"Leuk: Cox regression","title":"BUGS code for the Leuk example:","text":"","category":"section"},{"location":"bugs_examples/Leuk/","page":"Leuk: Cox regression","title":"Leuk: Cox regression","text":"model\n{\n # Set up data\n for(i in 1:N) {\n for(j in 1:T) {\n # risk set = 1 if obs.t >= t\n Y[i,j] <- step(obs.t[i] - t[j] + eps)\n # counting process jump = 1 if obs.t in [ t[j], t[j+1] )\n # i.e. if t[j] <= obs.t < t[j+1]\n dN[i, j] <- Y[i, j] * step(t[j + 1] - obs.t[i] - eps) * fail[i]\n }\n }\n # Model\n for(j in 1:T) {\n for(i in 1:N) {\n dN[i, j] ~ dpois(Idt[i, j]) # Likelihood\n Idt[i, j] <- Y[i, j] * exp(beta * Z[i]) * dL0[j] # Intensity\n }\n dL0[j] ~ dgamma(mu[j], c)\n mu[j] <- dL0.star[j] * c # prior mean hazard\n\n # Survivor function = exp(-Integral{l0(u)du})^exp(beta*z)\n S.treat[j] <- pow(exp(-sum(dL0[1 : j])), exp(beta * -0.5));\n S.placebo[j] <- pow(exp(-sum(dL0[1 : j])), exp(beta * 0.5)); \n }\n # Priors\n c <- 0.001\n r <- 0.1\n for (j in 1 : T) {\n dL0.star[j] <- r * (t[j + 1] - t[j])\n }\n beta ~ dnorm(0.0,0.000001)\n}","category":"page"},{"location":"bugs_examples/Leuk/#Reference-Results","page":"Leuk: Cox regression","title":"Reference Results","text":"","category":"section"},{"location":"bugs_examples/Leuk/","page":"Leuk: Cox regression","title":"Leuk: Cox regression","text":"Variable Mean Median Standard Deviation Monte Carlo Error 2.5% Value 97.5% Value Start Sample ESS\nS.placebo[1] 0.9264 0.9374 0.04989 3.349E-4 0.8029 0.9909 1001 20000 22184\nS.placebo[17] 0.04431 0.03344 0.03909 2.698E-4 0.002478 0.1487 1001 20000 20992\nS.treat[1] 0.9826 0.9863 0.01413 1.074E-4 0.9457 0.9982 1001 20000 17315\nS.treat[17] 0.4767 0.4763 0.1198 0.001009 0.2474 0.7086 1001 20000 14104\nbeta 1.539 1.524 0.4211 0.0034 0.7475 2.388 1001 20000 15340","category":"page"},{"location":"R_interface/#Integrating-R-in-Julia","page":"R Interface","title":"Integrating R in Julia","text":"","category":"section"},{"location":"R_interface/","page":"R Interface","title":"R Interface","text":"Julia offers a seamless interface to the R language. ","category":"page"},{"location":"R_interface/","page":"R Interface","title":"R Interface","text":"The RCall.jl package enables interaction with R functions in Julia.\nThe RData.jl package allows interfacing with R data in Julia.","category":"page"},{"location":"R_interface/#Reading-BUGS-data-and-init-from-R-like-lists","page":"R Interface","title":"Reading BUGS data and init from R like lists","text":"","category":"section"},{"location":"R_interface/","page":"R Interface","title":"R Interface","text":"Warning: The data layout in BUGS assumes that the data is stored in row-major order, while R uses column-major order. This discrepancy can lead to issues. Stan developers have transformed the data and initializations of example BUGS models for R, which can be found here.","category":"page"},{"location":"R_interface/#Reading-the-list-data-structure-from-R","page":"R Interface","title":"Reading the list data structure from R","text":"","category":"section"},{"location":"R_interface/","page":"R Interface","title":"R Interface","text":"The data for Rats is available here. ","category":"page"},{"location":"R_interface/","page":"R Interface","title":"R Interface","text":"In Julia, we can read this data into a Julia dictionary using the RCall.jl package.","category":"page"},{"location":"R_interface/","page":"R Interface","title":"R Interface","text":"julia> using RCall\n\njulia> data = R\"\nlist(\n x = c(8.0, 15.0, 22.0, 29.0, 36.0), xbar = 22, N = 30, T = 5,\n Y = structure(\n .Data = c(\n 151, 199, 246, 283, 320,\n 145, 199, 249, 293, 354,\n 147, 214, 263, 312, 328,\n 155, 200, 237, 272, 297,\n 135, 188, 230, 280, 323,\n 159, 210, 252, 298, 331,\n 141, 189, 231, 275, 305,\n 159, 201, 248, 297, 338,\n 177, 236, 285, 350, 376,\n 134, 182, 220, 260, 296,\n 160, 208, 261, 313, 352,\n 143, 188, 220, 273, 314,\n 154, 200, 244, 289, 325,\n 171, 221, 270, 326, 358,\n 163, 216, 242, 281, 312,\n 160, 207, 248, 288, 324,\n 142, 187, 234, 280, 316,\n 156, 203, 243, 283, 317,\n 157, 212, 259, 307, 336,\n 152, 203, 246, 286, 321,\n 154, 205, 253, 298, 334,\n 139, 190, 225, 267, 302,\n 146, 191, 229, 272, 302,\n 157, 211, 250, 285, 323,\n 132, 185, 237, 286, 331,\n 160, 207, 257, 303, 345,\n 169, 216, 261, 295, 333,\n 157, 205, 248, 289, 316,\n 137, 180, 219, 258, 291,\n 153, 200, 244, 286, 324\n ),\n .Dim = c(30, 5)\n )\n)\n\"\nRObject{VecSxp}\n$x\n[1] 8 15 22 29 36\n\n$xbar\n[1] 22\n\n$N\n[1] 30\n\n$T\n[1] 5\n\n$Y\n [,1] [,2] [,3] [,4] [,5]\n [1,] 151 141 154 157 132\n [2,] 199 189 200 212 185\n [3,] 246 231 244 259 237\n [4,] 283 275 289 307 286\n [5,] 320 305 325 336 331\n [6,] 145 159 171 152 160\n [7,] 199 201 221 203 207\n [8,] 249 248 270 246 257\n [9,] 293 297 326 286 303\n[10,] 354 338 358 321 345\n[11,] 147 177 163 154 169\n[12,] 214 236 216 205 216\n[13,] 263 285 242 253 261\n[14,] 312 350 281 298 295\n[15,] 328 376 312 334 333\n[16,] 155 134 160 139 157\n[17,] 200 182 207 190 205\n[18,] 237 220 248 225 248\n[19,] 272 260 288 267 289\n[20,] 297 296 324 302 316\n[21,] 135 160 142 146 137\n[22,] 188 208 187 191 180\n[23,] 230 261 234 229 219\n[24,] 280 313 280 272 258\n[25,] 323 352 316 302 291\n[26,] 159 143 156 157 153\n[27,] 210 188 203 211 200\n[28,] 252 220 243 250 244\n[29,] 298 273 283 285 286\n[30,] 331 314 317 323 324","category":"page"},{"location":"R_interface/","page":"R Interface","title":"R Interface","text":"alternatively, reval(s::String) will produce the same result in this case.","category":"page"},{"location":"R_interface/","page":"R Interface","title":"R Interface","text":"If the data is stores in a file, user can use function (may require customizing the function to fit specific needs)","category":"page"},{"location":"R_interface/","page":"R Interface","title":"R Interface","text":"function read_rlist_to_dictionary(filepath::String)\n r_data = open(filepath) do f\n s = read(f, String)\n reval(s)\n end\n return rcopy(r_data)\nend","category":"page"},{"location":"R_interface/","page":"R Interface","title":"R Interface","text":", and save the result to a Julia variable and access the data as a Julia dictionary","category":"page"},{"location":"R_interface/","page":"R Interface","title":"R Interface","text":"julia> rcopy(data)\nOrderedDict{Symbol, Any} with 5 entries:\n :x => [8.0, 15.0, 22.0, 29.0, 36.0]\n :xbar => 22.0\n :N => 30.0\n :T => 5.0\n :Y => [151.0 141.0 … 157.0 132.0; 199.0 189.0 … 212.0 185.0; … ; 298.0 273.0 … 285.0 286.0; 331.0 314.0 … 323.0 324.0]","category":"page"},{"location":"R_interface/","page":"R Interface","title":"R Interface","text":"It is worth noting that rcopy will automatically convert data names contains . to _ in Julia. E.g.","category":"page"},{"location":"R_interface/","page":"R Interface","title":"R Interface","text":"julia> rcopy(R\"list(a.b = 1)\")\nOrderedDict{Symbol, Any} with 1 entry:\n :a_b => 1.0","category":"page"},{"location":"R_interface/#Transform-Dta-read-from-R-to-Julia-convention","page":"R Interface","title":"Transform Dta read from R to Julia convention","text":"","category":"section"},{"location":"R_interface/","page":"R Interface","title":"R Interface","text":"If you want to load data using the R interface, but the data source is in the same layout as BUGS, you can process the data in Julia, for instance","category":"page"},{"location":"R_interface/","page":"R Interface","title":"R Interface","text":"# define a row-major reshape function, because Julia's `reshape` is column-major\njulia> function rreshape(v::Vector, dim)\n return permutedims(reshape(v, reverse(dim)), length(dim):-1:1)\n end \nrreshape (generic function with 1 method)\n\njulia> rreshape(vcat(data[:Y]...), (30, 5))\n30×5 Matrix{Float64}:\n 151.0 199.0 246.0 283.0 320.0\n 145.0 199.0 249.0 293.0 354.0\n 147.0 214.0 263.0 312.0 328.0\n 155.0 200.0 237.0 272.0 297.0\n 135.0 188.0 230.0 280.0 323.0\n 159.0 210.0 252.0 298.0 331.0\n 141.0 189.0 231.0 275.0 305.0\n 159.0 201.0 248.0 297.0 338.0\n ⋮ \n 146.0 191.0 229.0 272.0 302.0\n 157.0 211.0 250.0 285.0 323.0\n 132.0 185.0 237.0 286.0 331.0\n 160.0 207.0 257.0 303.0 345.0\n 169.0 216.0 261.0 295.0 333.0\n 157.0 205.0 248.0 289.0 316.0\n 137.0 180.0 219.0 258.0 291.0\n 153.0 200.0 244.0 286.0 324.0","category":"page"},{"location":"R_interface/","page":"R Interface","title":"R Interface","text":"Please always verify the data before using.","category":"page"},{"location":"#Introduction","page":"Introduction","title":"Introduction","text":"","category":"section"},{"location":"","page":"Introduction","title":"Introduction","text":"JuliaBUGS is a graph-based probabilistic programming language and a component of the Turing ecosystem. The package aims to support modelling and inference for probabilistic programs written in the BUGS language. ","category":"page"},{"location":"","page":"Introduction","title":"Introduction","text":"This project is still in its early stage, with many key components needing to be completed. ","category":"page"},{"location":"","page":"Introduction","title":"Introduction","text":"Please refer to the example for usage information and a complete example.","category":"page"},{"location":"#What-is-BUGS?","page":"Introduction","title":"What is BUGS?","text":"","category":"section"},{"location":"","page":"Introduction","title":"Introduction","text":"The BUGS (Bayesian inference Using Gibbs Sampling) system is a probabilistic programming framework designed for specifying directed graphical models. Unlike certain other probabilistic programming languages (PPLs), such as Turing.jl or Pyro, the focus of BUGS is on specifying declarative relationships between nodes in a graph, which can be either logical or stochastic. This means that explicit declarations of variables, inputs, outputs, etc., are not required, and the order of statements is not critical.","category":"page"},{"location":"#The-BUGS-Approach-and-Benefits","page":"Introduction","title":"The BUGS Approach and Benefits","text":"","category":"section"},{"location":"","page":"Introduction","title":"Introduction","text":"Loops in BUGS are essentially a form of \"plate notation,\" offering a concise way to express repetitive statements across many constant indices. Variables in BUGS are either the names of nodes within the program or constant parts of the \"data\" that must be combined with a model for instantiation.","category":"page"},{"location":"","page":"Introduction","title":"Introduction","text":"A BUGS model provides a comprehensive representation of the relationships and dependencies among a set of variables within a Bayesian framework. Our goal is to support BUGS programs as much as possible while also incorporating Julia-specific syntax enhancements.","category":"page"},{"location":"","page":"Introduction","title":"Introduction","text":"The key advantage of utilizing such a graph-based approach is the clarity it provides in understanding the dependencies and relationships within a complex system. These graphical models allow users to explicitly state the conditional dependencies between variables. This makes the model's structure and assumptions transparent, aiding both in the development and interpretation stages. Furthermore, using such a graphical approach makes it easier to apply advanced algorithms for model inference, as it enables more efficient computation by identifying and exploiting the structure of the model.","category":"page"},{"location":"example/#Example:-Logistic-Regression-with-Random-Effects","page":"Example","title":"Example: Logistic Regression with Random Effects","text":"","category":"section"},{"location":"example/","page":"Example","title":"Example","text":"We will use the Seeds model for demonstration. This example concerns the proportion of seeds that germinated on each of 21 plates. Here, we transform the data into a NamedTuple:","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"data = (\n r = [10, 23, 23, 26, 17, 5, 53, 55, 32, 46, 10, 8, 10, 8, 23, 0, 3, 22, 15, 32, 3],\n n = [39, 62, 81, 51, 39, 6, 74, 72, 51, 79, 13, 16, 30, 28, 45, 4, 12, 41, 30, 51, 7],\n x1 = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],\n x2 = [0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1],\n N = 21,\n)","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"where r[i] is the number of germinated seeds and n[i] is the total number of the seeds on the i-th plate. Let p_i be the probability of germination on the i-th plate. Then, the model is defined by:","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"$","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"\\begin{aligned} bi &\\sim \\text{Normal}(0, \\tau) \\\n\\text{logit}(pi) &= \\alpha0 + \\alpha1 x{1 i} + \\alpha2 x{2i} + \\alpha{12} x{1i} x{2i} + b{i} \\\nri &\\sim \\text{Binomial}(pi, ni) \\end{aligned} $","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"where x_1i and x_2i are the seed type and root extract of the i-th plate. The original BUGS program for the model is:","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"model\n{\n for( i in 1 : N ) {\n r[i] ~ dbin(p[i],n[i])\n b[i] ~ dnorm(0.0,tau)\n logit(p[i]) <- alpha0 + alpha1 * x1[i] + alpha2 * x2[i] +\n alpha12 * x1[i] * x2[i] + b[i]\n }\n alpha0 ~ dnorm(0.0, 1.0E-6)\n alpha1 ~ dnorm(0.0, 1.0E-6)\n alpha2 ~ dnorm(0.0, 1.0E-6)\n alpha12 ~ dnorm(0.0, 1.0E-6)\n tau ~ dgamma(0.001, 0.001)\n sigma <- 1 / sqrt(tau)\n}","category":"page"},{"location":"example/#Modeling-Language","page":"Example","title":"Modeling Language","text":"","category":"section"},{"location":"example/#Writing-a-Model-in-BUGS","page":"Example","title":"Writing a Model in BUGS","text":"","category":"section"},{"location":"example/","page":"Example","title":"Example","text":"BUGS language syntax: BNF definition","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"Language References: ","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"MultiBUGS\nOpenBUGS","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"Implementations in C++ and R:","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"JAGS and its user manual\nNimble","category":"page"},{"location":"example/#Writing-a-Model-in-Julia","page":"Example","title":"Writing a Model in Julia","text":"","category":"section"},{"location":"example/","page":"Example","title":"Example","text":"We provide a macro solution which allows users to write down model definitions using Julia:","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"model_def = @bugs begin\n for i in 1:N\n r[i] ~ dbin(p[i], n[i])\n b[i] ~ dnorm(0.0, tau)\n p[i] = logistic(alpha0 + alpha1 * x1[i] + alpha2 * x2[i] + alpha12 * x1[i] * x2[i] + b[i])\n end\n alpha0 ~ dnorm(0.0, 1.0E-6)\n alpha1 ~ dnorm(0.0, 1.0E-6)\n alpha2 ~ dnorm(0.0, 1.0E-6)\n alpha12 ~ dnorm(0.0, 1.0E-6)\n tau ~ dgamma(0.001, 0.001)\n sigma = 1 / sqrt(tau)\nend","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"BUGS syntax carries over almost one-to-one to Julia, with minor exceptions. In general, when basic Julia syntax and BUGS syntax conflict, it is necessary to use Julia syntax. For example, curly braces are replaced with begin ... end blocks, and for loops do not require parentheses. In addition, Julia uses f(x) = ... as a shorthand for function definition, so BUGS' link function syntax can be confusing and ambiguous. Thus, instead of calling the link function, we call the inverse link function from the RHS.","category":"page"},{"location":"example/#Support-for-Legacy-BUGS-Programs","page":"Example","title":"Support for Legacy BUGS Programs","text":"","category":"section"},{"location":"example/","page":"Example","title":"Example","text":"The @bugs macro also works with original (R-like) BUGS syntax:","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"model_def = @bugs(\"\"\"\nmodel{\n for( i in 1 : N ) {\n r[i] ~ dbin(p[i],n[i])\n b[i] ~ dnorm(0.0,tau)\n logit(p[i]) <- alpha0 + alpha1 * x1[i] + alpha2 * x2[i] +\n alpha12 * x1[i] * x2[i] + b[i]\n }\n alpha0 ~ dnorm(0.0,1.0E-6)\n alpha1 ~ dnorm(0.0,1.0E-6)\n alpha2 ~ dnorm(0.0,1.0E-6)\n alpha12 ~ dnorm(0.0,1.0E-6)\n tau ~ dgamma(0.001,0.001)\n sigma <- 1 / sqrt(tau)\n}\n\"\"\", true)","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"By default, @bugs will translate R-style variable names like a.b.c to a_b_c, user can pass false as the second argument to disable this. We still encourage users to write new programs using the Julia-native syntax, because of better debuggability and perks like syntax highlighting. ","category":"page"},{"location":"example/#Compilation","page":"Example","title":"Compilation","text":"","category":"section"},{"location":"example/","page":"Example","title":"Example","text":"For now, the compile function will create a BUGSModel, which implements LogDensityProblems.jl interface.","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"compile(model_def::Expr, data, initializations),","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"The function compile takes three arguments: ","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"the output of @bugs, \nthe data, and\nthe initializations of parameters.","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"initializations = Dict(:alpha => 1, :beta => 1)","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"then we can compile the model with the data and initializations,","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"model = compile(model_def, data, initializations)","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"LogDensityProblemsAD.jl defined some extensions that support automatic differentiation packages. For example, with ReverseDiff.jl","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"using LogDensityProblemsAD, ReverseDiff\n\nad_model = ADgradient(:ReverseDiff, model; compile=Val(true))","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"Here ad_model will also implement all the interfaces of LogDensityProblems.jl. LogDensityProblemsAD.jl will automatically add the interface function logdensity_and_gradient to the model, which will return the log density and gradient of the model. And ad_model can be used in the same way as model in the example below.","category":"page"},{"location":"example/#Inference","page":"Example","title":"Inference","text":"","category":"section"},{"location":"example/","page":"Example","title":"Example","text":"For a differentiable model, we can use AdvancedHMC.jl to perform inference. For instance,","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"using AdvancedHMC, AbstractMCMC, LogDensityProblems, MCMCChains\n\nn_samples, n_adapts = 2000, 1000\n\nD = LogDensityProblems.dimension(model); initial_θ = rand(D)\n\nsamples_and_stats = AbstractMCMC.sample(\n ad_model,\n NUTS(0.8),\n n_samples;\n chain_type = Chains,\n n_adapts = n_adapts,\n init_params = initial_θ,\n discard_initial = n_adapts\n )","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"This will return the MCMC Chain,","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"Chains MCMC chain (2000×40×1 Array{Real, 3}):\n\nIterations = 1001:1:3000\nNumber of chains = 1\nSamples per chain = 2000\nparameters = alpha0, alpha12, alpha1, alpha2, tau, b[16], b[12], b[10], b[14], b[13], b[7], b[6], b[20], b[1], b[4], b[5], b[2], b[18], b[8], b[3], b[9], b[21], b[17], b[15], b[11], b[19], sigma\ninternals = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size, is_adapt\n\nSummary Statistics\n parameters mean std mcse ess_bulk ess_tail rhat ess_per_sec \n Symbol Float64 Float64 Float64 Real Float64 Float64 Missing \n\n alpha0 -0.5642 0.2320 0.0084 766.9305 1022.5211 1.0021 missing\n alpha12 -0.8489 0.5247 0.0170 946.0418 1044.1109 1.0002 missing\n alpha1 0.0587 0.3715 0.0119 966.4367 1233.2257 1.0007 missing\n alpha2 1.3852 0.3410 0.0127 712.2978 974.1566 1.0002 missing\n tau 1.8880 0.7705 0.0447 348.9331 338.3655 1.0030 missing\n b[16] -0.2445 0.4459 0.0132 1528.0578 843.8225 1.0003 missing\n b[12] 0.2050 0.3602 0.0086 1868.6126 1202.1363 0.9996 missing\n b[10] -0.3500 0.2893 0.0090 1047.3119 1245.9358 1.0008 missing\n ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮\n 19 rows omitted\n\nQuantiles\n parameters 2.5% 25.0% 50.0% 75.0% 97.5% \n Symbol Float64 Float64 Float64 Float64 Float64 \n\n alpha0 -1.0143 -0.7143 -0.5590 -0.4100 -0.1185\n alpha12 -1.9063 -1.1812 -0.8296 -0.5153 0.1521\n alpha1 -0.6550 -0.1822 0.0512 0.2885 0.8180\n alpha2 0.7214 1.1663 1.3782 1.5998 2.0986\n tau 0.5461 1.3941 1.8353 2.3115 3.6225\n b[16] -1.2359 -0.4836 -0.1909 0.0345 0.5070\n b[12] -0.4493 -0.0370 0.1910 0.4375 0.9828\n b[10] -0.9570 -0.5264 -0.3331 -0.1514 0.1613\n ⋮ ⋮ ⋮ ⋮ ⋮ ⋮\n 19 rows omitted\n","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"This is consistent with the result in the OpenBUGS seeds example.","category":"page"},{"location":"example/#Parallel-and-Distributed-Sampling-with-AbstractMCMC","page":"Example","title":"Parallel and Distributed Sampling with AbstractMCMC","text":"","category":"section"},{"location":"example/","page":"Example","title":"Example","text":"AbstractMCMC and AdvancedHMC support both parallel and distributed sampling.","category":"page"},{"location":"example/#Parallel-Sampling","page":"Example","title":"Parallel Sampling","text":"","category":"section"},{"location":"example/","page":"Example","title":"Example","text":"To perform multi-threaded sampling of multiple chains, start the Julia session with the -t argument. The model compilation code remains the same, and we can sample multiple chains in parallel as follows:","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"n_chains = 4\nsamples_and_stats = AbstractMCMC.sample(\n ad_model,\n AdvancedHMC.NUTS(0.65),\n AbstractMCMC.MCMCThreads(),\n n_samples,\n n_chains;\n chain_type = Chains,\n n_adapts = n_adapts,\n init_params = [initial_θ for _ = 1:n_chains],\n discard_initial = n_adapts,\n)","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"In this case, we pass two additional arguments to AbstractMCMC.sample:","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"AbstractMCMC.MCMCThreads(): the sampler type, and\nn_chains: the number of chains to sample.","category":"page"},{"location":"example/#Distributed-Sampling","page":"Example","title":"Distributed Sampling","text":"","category":"section"},{"location":"example/","page":"Example","title":"Example","text":"To perform distributed sampling of multiple chains, start the Julia session with the -p argument.","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"In distributed mode, ensure that all functions and modules are available on all processes. Use @everywhere to make the functions and modules available on all processes.","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"For example:","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"@everywhere begin\n using JuliaBUGS, LogDensityProblems, LogDensityProblemsAD, AbstractMCMC, AdvancedHMC, MCMCChains, ReverseDiff # also other packages one may need\n\n # Define the functions to use\n # Use `@register_primitive` to register the functions to use in the model\n\n # Distributed can handle data dependencies in some cases, for more detail, see https://docs.julialang.org/en/v1/manual/distributed-computing/\n\nend\n\nn_chains = nprocs() - 1 # use all the processes except the master process\nsamples_and_stats = AbstractMCMC.sample(\n ad_model,\n AdvancedHMC.NUTS(0.65),\n AbstractMCMC.MCMCDistributed(),\n n_samples,\n n_chains;\n chain_type = Chains,\n n_adapts = n_adapts,\n init_params = [initial_θ for _ = 1:n_chains], # each chain has its own initial parameters\n discard_initial = n_adapts,\n progress = false, # Base.TTY creating problems in distributed setting\n)","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"In this case, we pass two additional arguments to AbstractMCMC.sample:","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"AbstractMCMC.MCMCDistributed(): the sampler type, and\nn_chains: the number of chains to sample.","category":"page"},{"location":"example/","page":"Example","title":"Example","text":"Note that the init_params argument is now a vector of initial parameters for each chain. Sometimes the progress logger can cause problems in distributed setting, so we can disable it by setting progress = false.","category":"page"},{"location":"example/#More-Examples","page":"Example","title":"More Examples","text":"","category":"section"},{"location":"example/","page":"Example","title":"Example","text":"We have transcribed all the examples from the first volume of the BUGS Examples (original and transcribed). All programs and data are included, and can be compiled using the steps described in the tutorial above.","category":"page"}] } diff --git a/previews/PR150/user_defined_functions/index.html b/previews/PR150/user_defined_functions/index.html index 36cbcbde9..f3148d2a2 100644 --- a/previews/PR150/user_defined_functions/index.html +++ b/previews/PR150/user_defined_functions/index.html @@ -11,4 +11,4 @@ julia> JuliaBUGS.@register_primitive(f); julia> JuliaBUGS.f(1) -2

      After registering the function or distributions, they can be used just like any other functions or distributions provided by BUGS.

      +2

      After registering the function or distributions, they can be used just like any other functions or distributions provided by BUGS.