Sampling from a PDF

[axsk]

Is there any way to sample from any non-standard distribution P (i.e. not contained in Distributions.jl) which is given by its probability density function (PDF), or any proportional?

I wanted to try Mamba.jl on a simple Bayesian inference example from my class: Given by the PDF where f is an arbitrary map representing some measurement, y the measured value, and ϕ the PDF of a normal distribution representing noise in the measurement.

As far as I understood MCMC all it needs is a PDF to produce the samples. Still, reading the manual, I could not find any way to sample from the above example.

One possibility would be specifying a new Distribution, but according to Distribtions.jl's manual this would need a rand(d::D) method, sampling from the distribution, but this is what I wanted to use Mamba for, so I feel stuck.

Is there any way to solve this with Mamba, and if not, shouldn't there be?

[brian-j-smith]

Sure. Probably the most straightforward approach would be to use one of the stand-alone samplers, like that for slice sampling, as follows.

P.S. The code below requires julia 0.4.

## Data
y = 5.0

## Map
f(x) = x

## Normalized or unnormalized prior density function
P(x) = 1.0

## Normalized or unnormalized log-density function whose single argument x is a
## draw from the distribution.
logf = function(x::DenseVector)
  log(P(x)) - 0.5 * sumabs2(y - f(x))
end

## Slice sampling from the distribution
n = 10000
sim = Chains(n, 1, names = ["x"])
x = SliceVariate([0.0])
width = fill(1.0, length(x))
for i in 1:n
  slice!(x, width, logf)
  sim[i,:,1] = x
end
describe(sim)

[axsk]

Thank you for your answer.

I hoped to be able to use the general Model() approach. Digging further in the manual I found how to specify User-Defined Univariate Distributions, so I guess this would be the canonical approach? I think I could define a Posterior <: ContinuousUnivariateDistribution Distribution corresponding to your logf above, but could I split this up into Posterior, Prior and Likelihood nodes, the first being dependent on the latter ones?

Also looking at all those eval tricks I wonder what makes them necessary. Wouldn't it be possible to pass the node code directly as function with either the ::Model argument or directly with its input nodes, e.g. y = Stochastic(1, (mu, s2) -> MvNormal(mu, sqrt(s2)))?

In that case the functions inside the definition should be resolved in the current namespace without any needs for the eval hack, which I think is much cleaner.

[brian-j-smith]

A Model() approach can definitely be used:

• y would be a stochastic node having a Normal distribution with mean f(x) and standard deviation 1.0 - the likelihood.
• x would also be stochastic, with its non-standard distribution P(x) being user-defined - the prior.

As for functions... accepting them in the node specifications would make for neat syntax. My preference for such functions is the signature used in your example (input nodes as the arguments). To be honest, I started this package as my first julia program back in julia 0.2, and how to do that function approach wasn't clear to me back then. Now that it is, I am considering it. The implementation would require working directly with function AST representations to extract the input node symbols and link the supplied function to the model dictionary of nodes.

Accepting functions with a ::Model argument is something that I did consider, but decided against because the syntax would be even messier, more redundant (same function signature for each node), and more prone to function misspecifications (invariably followed by issues being opened by confused users).

[axsk]

So the solution is that easy :)

I am happy that you like the idea with the input nodes as arguments. It also seems to be a rather local change so I guess I will give it a shot the coming days.

I am specifically interested in this to simplify specifying Stochastic nodes. Given the above syntax it should be straightforward to define a function bypdf(f::Function) returning a Distribution with the desired logdensity (maybe this should be a feature of Distributions.jl). The example above could then be written like this:

## simple uniform prior defined via pdf
x = Stochastic(bypdf(x -> 1))
## the likelihood, dependent on the prior with user-specified density (here proportional to the normal pdf)
y = Stochastic(x -> bypdf(y -> -exp(sumabs2(x-y)))

Excited to get this to work :)

[brian-j-smith]

I now have an implementation, running locally, that allows functions in the specification of Stochastic and Logical nodes. A preview of the syntax is shown below for part of the line example from the tutorial. I'll include this as an experimental (until the documentation is updated) feature in the next version of Mamba, planned for release this weekend. The current (expression/macro) syntax will continue to work - the function syntax is just an alternative way to do the specifications.

Just make sure that the functions/expressions in your stochastic nodes return Distribution type objects, like those defined in the Distributions package. For example, syntax for a simple flat (unconstrained uniform) prior would look like x = Stochastic(() -> Flat()); a normal with mapping f(x) would be y = Stochastic(x -> Normal(f(x), 1.0); or other package or user-defined distributions could be returned by the node functions.

using Mamba

## Model Specification

model = Model(

  y = Stochastic(1,
    (mu, s2) -> MvNormal(mu, sqrt(s2)),
    false
  ),

  mu = Logical(1,
    (xmat, beta) -> xmat * beta,
    false
  ),

  beta = Stochastic(1,
    () -> MvNormal(2, sqrt(1000))
  ),

  s2 = Stochastic(
    () -> InverseGamma(0.001, 0.001)
  )

)

## User-Defined Samplers

Gibbs_beta = Sampler([:beta],
  (beta, s2, xmat, y) ->
    begin
      beta_mean = mean(beta.distr)
      beta_invcov = invcov(beta.distr)
      Sigma = inv(xmat' * xmat / s2 + beta_invcov)
      mu = Sigma * (xmat' * y / s2 + beta_invcov * beta_mean)
      rand(MvNormal(mu, Sigma))
    end
)

Gibbs_s2 = Sampler([:s2],
  (mu, s2, y) ->
    begin
      a = length(y) / 2.0 + shape(s2.distr)
      b = sumabs2(y - mu) / 2.0 + scale(s2.distr)
      rand(InverseGamma(a, b))
    end
)

## User-Defined Sampling Scheme
scheme3 = [Gibbs_beta, Gibbs_s2]

## Data
line = Dict{Symbol,Any}(
  :x => [1, 2, 3, 4, 5],
  :y => [1, 3, 3, 3, 5]
)
line[:xmat] = [ones(5) line[:x]]


## Set Random Number Generator Seed
srand(123)


## Initial Values
inits = [
  Dict{Symbol,Any}(
    :y => line[:y],
    :beta => rand(Normal(0, 1), 2),
    :s2 => rand(Gamma(1, 1))
  )
  for i in 1:3
]


## MCMC Simulations

setsamplers!(model, scheme3)
sim3 = mcmc(model, line, inits, 10000, burnin=250, thin=2, chains=3)

sim = mcmc(sim3, 5000)
describe(sim)

[brian-j-smith]

Function syntax for the specification of Model nodes is officially supported now in the latest release of Mamba (0.7.0), along with full documentation and updated tutorial and examples.

[mikejacktzen]

Hi Brian and axsk,

I think this is a good opportunity to check in. Awhile ago, I wanted to implement a custom many to one function to use as a univariate data likelihood. The caveat was, the many to one function needed to make use of every entry of a data vector for the evaluation at each univariate scalar parameter. I think I found a solution (it works)! I get expected behavior using my custom likelihood (without making use of mamba eg only relying on Distributions) AND successful posterior output (making use of mamba).

I was hoping Brian might be able to chime in on my maybe successful example. In essence, I worked around the very stringent Distributions syntax to define the custom likelihood. The key trick was to modify the function signature in logpdf() below. The arithmatic in the body of the custom logpdf() is removed for clarity.

Since my main question is if the modification of the function signature of the logpdf() is correctly implemented (distinguish if what I'm exploiting is a feature or a bug)

 ## Type declaration

  type CustomUnivarDist <: ContinuousUnivariateDistribution
        theta::Float64

    ## Constructor
        function CustomUnivarDist(theta)
            new(convert(Float64, theta))
    end
  end

  ## The following method functions must be implemented

  #################################
  ## Minimum and maximum support values
  #################################
  minimum(d::CustomUnivarDist) = -Inf
  maximum(d::CustomUnivarDist) = Inf

#################################
# logpdf
#################################

# old/unmodified mamba example
# function logpdf(d::CustomUnivarDist, y::Real)

# successful modification of function signature
function logpdf{T<:Real}(d::CustomUnivarDist, y::Vector{T})

# arbitrary computation of 
# many-to one mapping eg  input y to real
# goes here
# removed for clarity

    end
  end

export CustomUnivarDist

[brian-j-smith]

Hi @mikejacktzen,

Congratulations on the progress! Have you tried your modified logpdf approach in recent versions of julia (0.4) and Mamba (0.7)? The package is a little more stringent now about the logpdf signature. I would suggest trying a multivariate approach (see skeleton code below), as a check on and/or alternative to the univariate approach.

  using Distributions
  import Distributions: length, insupport, _logpdf
  export CustomMultivarDist

  type CustomMultivarDist <: ContinuousMultivariateDistribution
        n::Int
        theta::Float64

    ## Constructor
        function CustomMultivarDist(n::Int, theta::Real)
            new(n, convert(Float64, theta))
    end
  end

  ## The following method functions must be implemented

  ## Dimension of the distribution
  length(d::CustomMultivarDist) = d.n

  ## Logical indicating whether x is in the support
  function insupport{T<:Real}(d::CustomMultivarDist, x::Vector{T})
    length(d) == length(x)
  end

  ## Normalized or unnormalized log-density value
  function _logpdf{T<:Real}(d::CustomMultivarDist, x::Vector{T})
    ## Insert log-density calculations here.
  end

[axsk]

I created a UnivariateDensityDistribution(similar to your code above) and sampling it works fine using the Slice sampler, but when attempting to use NUTS I receive a method error:

ERROR: MethodError: `start` has no method matching start(::UnivariateDensityDistribution)
 in mapfoldl at reduce.jl:56
 in link at C:\Users\Alex\.julia\v0.4\Mamba\src\distributions/transformdistribution.jl:7
 in link_sub at C:\Users\Alex\.julia\v0.4\Mamba\src\distributions/distributionstruct.jl:86
 in anonymous at C:\Users\Alex\.julia\v0.4\Mamba\src\distributions/distributionstruct.jl:90
 in map! at abstractarray.jl:1251
 in link_sub at C:\Users\Alex\.julia\v0.4\Mamba\src\distributions/distributionstruct.jl:90
 in unlist at C:\Users\Alex\.julia\v0.4\Mamba\src\model/dependent.jl:194
 in unlist at C:\Users\Alex\.julia\v0.4\Mamba\src\model/dependent.jl:190
 in anonymous at C:\Users\Alex\.julia\v0.4\Mamba\src\model/simulation.jl:155
 in map at abstractarray.jl:1279
 in unlist at C:\Users\Alex\.julia\v0.4\Mamba\src\model/simulation.jl:155
 in unlist at C:\Users\Alex\.julia\v0.4\Mamba\src\model/simulation.jl:142
 [inlined code] from C:\Users\Alex\.julia\v0.4\Mamba\src\samplers/nuts.jl:52
 in anonymous at no file:51
 in simulate! at C:\Users\Alex\.julia\v0.4\Mamba\src\model/simulation.jl:132
 in simulate! at C:\Users\Alex\.julia\v0.4\Mamba\src\model/simulation.jl:124
 in mcmc_worker! at C:\Users\Alex\.julia\v0.4\Mamba\src\model/mcmc.jl:68
 in map at abstractarray.jl:1279
 in mcmc_master! at C:\Users\Alex\.julia\v0.4\Mamba\src\model/mcmc.jl:47
 in mcmc at C:\Users\Alex\.julia\v0.4\Mamba\src\model/mcmc.jl:32
 in ...

So I guess one still needs to supply the start and end values? They are not denoted as necessary in the docs yet :(

Edit: Excuse the fuss. Actually I missed that custom univar. distr. need to supply the minimum and maximum methods, as correctly denoted in the docs. I was not providing these, so that Julia fell back to the Base ones, which of course make no sense on a distribution.

[bdeonovic]

Glad you figured out your problem @axsk. Let us know if you have any more questions. By the way did you see this issue: https://github.com/brian-j-smith/Mamba.jl/issues/42, there were lots of helpful comments on their about setting up a custom distribution.

[mikejacktzen]

@brian-j-smith

I have confirmed that updating julia to 0.4 and Mamba to 0.7 has broken my original 'maybe-working' solution of a customized density that computes on a vector of scalar observations for a scalar parameter.

I'm going to try out the modifications you suggested, and report back.

[mikejacktzen]

@brian-j-smith

Reporting back! the multivariate solution with the newer signature restrictions via Mamba 0.7 work!

I think the new mamba signature restrictions is clearer. In previous mamba pre 0.7, i resorted to the univariate. Since, the multivariate method obfuscated the way to specify the dimensions of the multivariate distribution as the number of data observations (this was giving me incomprehensible errors when crosswalking to the dag model specification)

Maybe it was present pre 0.7, but it was hard to figure out with the work in progress documentation.

Brian's suggestion of

Dimension of the distribution

length(d::CustomMultivarDist) = d.n

took care of that

The new, "functions in node syntax" is great too!

[brian-j-smith]

Awesome! So glad to hear that, and thanks for the updates along the way.

[axsk]

Is there any advantage of using Univariate Dist's over 1d Multivariates?

Edit: I just saw that Univariates support a "transform" bool, but don't understand what it is doing, nor find anything about it in the manual. What is it doing?

[brian-j-smith]

Some of the samplers (those listed as applying transformations in this table) are designed for model parameters that are defined on or can be transformed to an unconstrained space. Mamba has built-in transformations that will work in general for Univariates. However, transformations for Multivariates tend to be distribution specific and not generally implementable. The "transform" bool allows the aforementioned samplers to apply transformations that have been defined for parameters being sampled.

If a user-defined Multivariate is created with a constrained space, can be transformed to an unconstrained space, and it is desirable to do so for compatibility with samplers that operate on unconstrained spaces; then link/invlink transformation functions can be defined for the distribution. An example of this can be seen in the code for positive definite matrices.

Apart from the transformation issue, there may be computational advantages to one over the other. Those will primarily depend on the performance of summing logpdf over individual Univariate elements versus a single Multivariate logpdf calculation. For instance, independent Normals could be represented by a multivariate Normal or by a vector of Univariate Normals. The multivariate calculation may involve matrix inversion and be slower than separate Univariate calculations. Which approach is better depends on the application and logpdf implementation more so than on any overhead differences coming from Mamba.

[mikejacktzen]

Hi @brian-j-smith and @axsk ,

I'm following up this thread by trying to generalize my successful example of a vector containing univariate outcome with logpdf evaluation of a scalar parameter into a more general regression setting (vector containing univariate outcome with logpdf evaluation over a vector of regression parameters).

The custom logpdf needs be evaluated over a vector of regression parameters using the entire data vector in it's function body computation.

I'm making some iterative progress but hitting errors. I think the issue touches on the difference / benefits of 'univariate' versus 'multivariate'

First, my data of outcomes y and input covariates x are columnwise concatenated into a single array called 'yx_mat' where the first column is the outcome y, and the following columns are the input covariates x. rows index the observation points

I think this may be related to 'arrays of multivariate distributions' but i can't think through it clearly.

I've tried modifying the examples above with the following logpdf modifications

The constructor is slightly modified to a vector of thetas (for the parameters of covariates x)

    type CustomMultivarDist2 <: ContinuousMultivariateDistribution
        n::Int
        theta::Vector{Float64}

        ## Constructor
        function CustomMultivarDist2(n::Int, theta::AbstractVector)
            new(n, convert(Vector{Float64}, theta))
        end
    end

I believe the above is correct for the constructor, the tricky part is modifying the logpdf below Sequentially, i'ved tried out three signatures for _logpdf{}()

    function _logpdf{T<:Real}(d::CustomMultivarDist2, yx_mat::Vector{T})
    function _logpdf{T<:Real}(d::CustomMultivarDist2, yx_mat::Array{T})
    function _logpdf{T<:Real}(d::CustomMultivarDist2, yx_mat::AbstractArray{T})

which was swapped into the signature of the 'body' of the logpdf function below

    # issue: figure out how to feed in output y and input x variables    
    # proper? signature for yx_mat as an array, y in first column, x in rest

    function _logpdf{T<:Real}(d::CustomMultivarDist2, yx_mat::AbstractArray{T})

    ## Insert log-density calculations here.

    ## split up yx_mat into seperate outcome y and input x

        y = yx_mat[:,1]  # first column is outcome y
        x = yx_mat[:,2:end]  # second till end is input x

    # toy example using matrix algebra

    # for pure illustration: inner product of residuals demonstrates requirement of using ALL data observations
     dot(y - x * d.theta, y - x * d.theta)

  end

The last signature spits out the least amount of errors (yx_mat::AbstractArray{T}).

Would you all have some suggestions on how to proceed? I'm thinking there might be a better way to feed in 'yx_mat' or the correct way of 'typing' it that I am overlooking.

[bdeonovic]

It looks like yx_mat is a matrix, which could be why yx_mat::AbstractArray{T} or yx_mat::Vector{T} would throw errors. In this case you are trying to define a Matrix distribution (similar to the Wishart Distribution). Is this correct? Check out Distributions.jl source code on the Wishart distribution to get an idea of how to write up an appropriate logpdf: https://github.com/JuliaStats/Distributions.jl/blob/master/src/matrix/wishart.jl

Without seeing what the rest of your logpdf is suppose to calculate I can't tell you much more than this. You might want to ask yourself, is this distribution really a matrix distribution? If for example you simply need the data X and y in the logpdf to calculate the residuals perhaps the distribution should just be defined on real values instead.

Also you need to post the errors you are getting, can't really help if we don't know what the error is.

Since this isn't a Mamba specific error I'd also like to refer you to some great reference material for Julia which you might not be aware of:

Julia Manual : Skimming through all of these sections can be very enlightening on how Julia works. I would especially recommend the chapter on Types since it looks like you are having some issues with those.

Julia Tutorial : This tutorial is quite short and excellent for getting the fundamentals down.

[mikejacktzen]

@bdeonovic

Thanks for the wishart.jl example

I was under the impression that in pure julia, everything is an array, even if 'matrices' are operated in matrix fashion. But the wishart.jl source code appears to explicitly declare matrix types.

In regards to residuals, it is for purely illustrative purposes (a mechanism for us statisticians to quickly get on the same page). The illustrative and end-use goal is to have a single overarching logpdf() computation that needs to involve all data points, parameters, outcomes, and inputs, into a single logpdf() evaluation.

The criteria above definately muddys the boundary between what's univariate, multivariate and what should be scalars, vectors, or matrices.

I think the situation here is integrating the transition between "Distributions" syntax and "Mamba" syntax.

[bdeonovic]

1. Julia does differentiate between "vectors" and matrix of 1 column, although most functions will know how to handle both, e.g. note:

julia> [1,2,3]                                                                                                                                                                                                                            
3-element Array{Int64,1}:
 1
 2
 3

julia> reshape([1,2,3],3,1)                                                                                                                                                                                                               
3x1 Array{Int64,2}:
 1
 2
 3

2. That said both are instances of AbstractArray so your logpdf signature should work with AbstractArray{T,2} and Array{T,2}, (but not Vector{T} which is an alias for AbstractArray{T,1}).

Again without knowing what errors you have I can't really help, but it sounds like the issues you are having are with programming up the distribution, in which case you might want to post at Distributions.jl.

[axsk]

I had the same problem (creating a matrix-variate distribution given some logpdf function). This is what I came up with, guess it might help:

type MultivariateDensityDistribution <: ContinuousMultivariateDistribution
  lpdf::Function
  dim
  insupport::Function
end

logpdf(d::MultivariateDensityDistribution, x::DenseMatrix, transform::Bool=true) = d.lpdf(x)
insupport(d::MultivariateDensityDistribution, x::DenseMatrix) = d.insupport(x)
length(d::MultivariateDensityDistribution)               = d.dim
size(d::MultivariateDensityDistribution)                 = d.dim

# Constructor for a DensityDistribution with logpdf: R^(2x2) -> R, x |-> ||x||²
MultivariateDensityDistribution(x -> sumabs2(x), (2,2), (x) -> true)

Working out the right method signatures was a pain though. There is so many different logpdf's flying around in Mamba.jl and Distributions.jl, it realy made me see multiple dispatch in a different light, as it was hard to work out which is responsible for what and called by whom.

Regarding length and size, if I remember correctly one should give the scalar dimension while the other gave each ranks dimension (i.e. 4 vs. 2x2). Probably only one is needed in the DenseMatrix case, but I didn't dare to further debug this, just leaving it like this.

The type should probably be called MatrixvariateDensityDistribution (named it before I sorted out the difference..)

In the end it would probably be nice to have a single DensityDistribution(dim::Union{Int, Tuple}, pdf::Function; log::Bool=true, insupport::Function=(x->true)) constructor, allowing the end-user to forget about the interna.

Maybe we could also let setinits! set the dim field to the size of the matrix, as the information in the constructor is redundant in that case. I have mixed feelings about this though as it would be mutating the state :/ What do you think about that?

If there's any interest I could create a PR for this.

[mikejacktzen]

@bdeonovic

Just to report back, my particular case was resolved with the correct typing of AbstractArray{T,2}

@axsk I like the idea of a abstract density distribtuion where the user could specify the dimension.

For your proposal of

DensityDistribution(dim::Union{Int, Tuple}

What if it also allows for scalar, eg dimension 0 (i think julia understands scalars as an array of dimension 0).

[mikejacktzen]

Following up this discussion: after some months of collecting dust, I think I've gotten the 'general' case to work via the matrix-variate way of @axsk The mamba sampler is running right now, and i'll followup with a confirmation.

I just wanted to +1 the suggestion of @axsk

Motivated by the discussions above, I think having a general api for the end user to specify the dimensions of the "computational input data" and the dimensions of the "mamba-ready output statistical distribution", would be helpful.

the context is: the user wants to feed all of the 'computational data' in one-shot as input into a custom statistically-formal univariate/multivariate/matrix-variate distribution (to be used as a stochastic node in mamba). Explicitly, the underlying gray-zone is matching the 'computational' data type-dimension with the 'intended/statistical' type-dimension recognized by mamba

From my iterative development process, I went down the following path

1. vector input data of scalars to model a univariate stiatistical model

• type CustomUnivarDist <: ContinuousUnivariateDistribution (in Mamba+Distributions)
• note: this was successful all the way thru sampling

2. matrix input data to model a univariate statistical model (regression onto covariate)

• type CustomMultivariateDist <: ContinuousMultivariateDistribution (in Mamba+Distributions)
• note: the evaluation of _logpdf was successful, but did not work well in a stochastic node of a mamba model (typing is the culprit)

3. matrix input data to model a univariate statistical model (regression onto covariate) - same as above

• type CustomMatrixDist <: ContinuousMatrixDistribution (in Mamba+Distributions)
• note: this is the hopefully successful implementation, as it's sampling right now

So the pattern was, to incrementally try higher dimensional 'arrays' on the Distributions+Mamba side in order to ingest computational data in matrix form.

I think we see leftover 'artifacts' of this as evidenced by how the mamba documentation is pretty 'heavy' for custom uni-variate distributions, but 'light' on the custom matrix-variate distributions

[bdeonovic]

I'm glad you got your code working. Without being able to see any code I can't really give any helpful comments, but are you aware that not all your data have to be nodes in mamba model? If it doesn't have a distribution (such as regression covariate information) then you can just put it in the inits dict. Mamba will have access to those variables.

For example see the rats example note that variables like the covariate X matrix and the centered covariates Xm are stored in the rats data dict. These are then referenced in the stochastic nodes, e.g.

y = Stochastic(1,
    (alpha, beta, rat, Xm, s2_c) ->
      begin
        mu = alpha[rat] + beta[rat] .* Xm
        MvNormal(mu, sqrt(s2_c))
      end,
    false
  ),

Let me know if this helps, or if I didn't understand your particular use case.

[mikejacktzen]

Hi @bdeonovic

Thanks for your ongoing helpful comments. Good news is, I think my current general implementation worked. The posterior samples seem correct.

I'm preparing a first draft write up of the results, and once ready I'll definitely share my code. (Scooping purposes).

That's a nifty trick I was not aware of. I think putting 'non distributional data' as inits may be an elegant way to achieve what I need to do. Specifically, combining that with the ability to 'pass in' a user-defined function into the stochastic node may be reasonable

[bdeonovic]

I'm glad everything is working well. It's not really a trick but how the package was designed to work. Reading through the tutorial and the examples should be quite illuminating. Good luck with your project.