Random variable constraints

[mrquincle]

Is it possible to add constraints to the walkers?

For example I'm implementing fitting a line segment at https://github.com/mrquincle/gwmcmc/blob/master/examples/linesegmentfit.m

Currently, how I handle the fact that sigma > 0, is by manually adding the following to the logprior function:

    if (sigma < 0)
        result = -Inf;
        return;

Thanks!

[grinsted]

That is also how i usually do it. Or I use something similar to this: logPrior=@(sigma,theta)-1/( (sigma>0)&(theta<-5) )

if the prior is entirely a boolean expression.

There is currently no other way of doing it. I dont really want to complicate the inner loop in the sampling with anything that might slow it down. But i could experiment with either an additional optional "Constraints" parameter that is then inserted into logPfuns behind the scenes. Or i could check if the output of logPfuns is a logical type and treat it as a constraint if it is. A simple hack would be to detect logical output after calculating the logP state of minit.

-Alternatively i often transform my parameters to a space with no limits. Using logarithm for sigmas, and logit for probabilities. But that is usually because i think it is easier to formulate a prior in the transformed space. Something like that could also be built into constraints, but it would be slower, and more complicates and probably not as flexible.

I'd be happy to have any input on this.

[grinsted]

I added a simple experimental implementation that converts logical functions to functions that return -inf. It works on this example but i have not tested further than that:

mu = [5;-3;6]; C = [.5 -.4 0;-.4 .5 0; 0 0 1]; iC=pinv(C); constraints=@(m)m(1)>4; logL=@(m)-0.5_sum((m-mu)'_iC*(m-mu)); logPfuns={constraints logL};

%make a set of starting points for the entire ensemble of walkers minit=rand(length(mu),length(mu)*2)+10;

%Apply the MCMC hammer [models,logP]=gwmcmc(minit,logPfuns,100000); models(:,:,1:floor(size(models,3)*.2))=[]; %remove 20% as burn-in models=models(:,:)'; %reshape matrix to collapse the ensemble member dimension scatter(models(:,1),models(:,2)) prctile(models,[5 50 95])

I still consider it an experimental feature and so have not documented this fully yet.