tuneSample is unpredictably giving errors during parameter tuning

Hi,

I am trying to use a hamiltonian MCMC to fit exponential line to multinomial data. This means there are several rows to my data set but I treat each row individually (iterating over the rows and treating the data as binomial).

As I am sure you are all aware, MCMC algorithms tend to misbehave if they are not tuned. As my task involves many iterations of fitting, looking at different subsets of the data each time, I cannot observe the behaviour of each MCMC fit. Below I give test data with expontential lines. The data set I want to use this fitting method on is much larger (don't worry I have access to high performance computer).

Test data that is I use when I get the errors:

td = round([normrnd(2,1, [5,2]),...
    [2*ones(1,3);zeros(4,3)],...
    normrnd(2,1, [5,5]),...
    exp([0.5;0.6;0.7;0.8;0.9]*[0,1,2,3,4,5,6,7,8,9,10]),...
    repmat(exp([0.5;0.6;0.7;0.8;0.9]*10), [1,10])+normrnd(0,25,[5,10])]);

My real data is count data so I do correct for negative values that slip in when generating this data.

The following code is conducted within three nested loops that vary line length, data subset. I can provide all the code if that is wanted for people to recreate the errors.

t = a vector of varible length in a loop not included here

l = a vairble I iterate over in a loop not included here

n = sum(td(:,t),1)
k = td(l,t)

logpdf = @(alpha) BetaGrad(k,n,alpha,t) %see below for function
hmc = hmcSampler(logpdf, 0);
hmc = tuneSampler(hmc);

I try to use estimateMAP for a starting position but it sometimes (and I don't know what the pattern is) will give results far too large (so above I have just used 0. Which is reasonable). I have noticed the rididulous results agree with the sampler when it is not properly tuned. As I said in the title somtimes the tuner works, sometimes it doesn't.

It is on the lower values of the data that the tuner doesn't work.

The most common error sequence is:

Error using rand
Complex or sparse input following 'like' is not supported.

Error in unifrnd (line 34)
r = mu + sig .* (2*rand(sizeOut,'like',mu)-1);

Error in stats.mcmc.impl.HamiltonianImpl.hmcParameters (line 241)
                        stepsize = unifrnd(0,maxstepsize,1,1);

Error in stats.mcmc.impl.HamiltonianImpl.doHMC (line 561)
                [stepsize,numsteps] =
                stats.mcmc.impl.HamiltonianImpl.hmcParameters(maxstepsize,maxnumsteps,jittercode);

Error in stats.mcmc.impl.HamiltonianImpl/drawSamples (line 66)
            [xsmpl,accratio,stepSizeTuningInfo,state] =
            stats.mcmc.impl.HamiltonianImpl.doHMC(start,numsamples,this.LogPDF,burnin,stepsize,numsteps,massvec,...

Error in stats.mcmc.impl.HamiltonianImpl/tuneMassVectorIterativeSampling (line 154)
                [xsmpl, finalState] = drawSamples(this, chainInfo, samplerParamsNew);

Error in stats.mcmc.impl.HamiltonianImpl/tuneMassVector (line 99)
                [samplerParams, massVectorTuningInfo] =
                tuneMassVectorIterativeSampling(this,chainInfo,samplerParams);

Error in stats.mcmc.impl.HamiltonianImpl/tuneSampler (line 76)
            [samplerParams, massVectorTuningInfo] =
            tuneMassVector(this,chainInfo,samplerParams);

Error in stats.mcmc.HamiltonianSampler/tuneSampler (line 568)
            [this.Impl,tuningInfo] = tuneSampler(this.Impl,chainInfo,samplerParams);

As you can see the root error is very deep in the tuning process so I am doing somthing fundementally wrong, but I have no idea what. As I will reitterate, this only happens sometimes and, from what I can tell, on lower data values exclusively.

I also get warnings about colons usage within the tuneSampler function

Warning: Colon operands must be real scalars. 
> In stats.mcmc.impl/HamiltonianImpl/doLeapFrog (line 284)
In stats.mcmc.impl/HamiltonianImpl/doHMCProposal (line 605)
In stats.mcmc.impl/HamiltonianImpl/doHMCOneStep (line 639)
In stats.mcmc.impl/HamiltonianImpl/tuneHMCDualAveraging (line 393)
In stats.mcmc.impl/HamiltonianImpl/doHMC (line 533)
In stats.mcmc.impl/HamiltonianImpl/drawSamples (line 66)
In stats.mcmc.impl/HamiltonianImpl/tuneMassVectorIterativeSampling (line 154)
In stats.mcmc.impl/HamiltonianImpl/tuneMassVector (line 99)
In stats.mcmc.impl/HamiltonianImpl/tuneSampler (line 76)
In stats.mcmc/HamiltonianSampler/tuneSampler (line 568)

Here is the logpdf and its gradient as requested by the sampler.

the pdf is a beta distribution sterBeta being the sterling approximation of the log of nchoosek becuase some of my factoirals are massive.

function [lpdf, glpdf] = BetaGrad(k,n,alpha,t)

% parameters of the prior
a = 1;
b = 2;


lpdf = sum(sterBeta(n,k,a,b) + alpha*t.*(a+k-1) - (a+k-1).*log(n)...
            + (b+n-k-1).*log(1-exp(alpha*t)./n));
        
% The derivative of a sum of functions is a sum of their derivatives        
glpdf = sum(t.*(a+k-1) - (b+n-k-1).*t.*exp(alpha*t)./(n-exp(alpha*t)));


end