MCMC sampling (NUTS) MatlabStan

[Hamed Nikbakht]

Hi all,

I'm a newcomer in Stan. I wanted to implement the Stan on a toy example.
I want to make sure whether I'm on the right track.

This is my example: My problem's dimension is 1000 and my target distribution is a standard normal distribution and I have a vector x (size 1000x1) which i set it as the initial values of my problem for the HMC NUTS sampling.

I also attach the figure to show you what i mean by saying the initial values for each chain.

 

(Lately, I have done the whole process with MCMC metropolis in Matlab but now I want to use this feature of Stan to make it more efficient with the NUTS.)

I would gratefully appreciate it if you could let me know whether I'm doing it correctly or not.

N = dimension of my problem, and ux = Initial state of each dimension for the Markov Chain. I'm trying to implement the HMC NUTS sampling technique
and the density (target distribution) will be the standard normal distribution as you can see in the following code.

load('Ufile.mat');  // Load my initial states which is the vector of ux
N = length(ux);
MCMC_dat = struct('N',N);
initial = struct('x',ux');
fit = stan('file','MC.stan','data',

MCMC_dat,'chains',1,'iter',10,'init',initial,'algorithm','NUTS','sample_file','output-6.csv');
%print(fit) 


This is my Stan code:
data {
  int<lower=1> N;

}
parameters {
  real x[N];
}
model {
for (n in 1:N)
  x[n] ~ normal(0, 1);
}

Thank you very much in advance for your help.
Regards,
Hamed

[Jonah Gabry]

Hi Hamed, 

Some comments below:

On Friday, March 3, 2017 at 12:17:53 PM UTC-5, Hamed Nikbakht wrote:

Hi all,

I'm a newcomer in Stan. I wanted to implement the Stan on a toy example.
I want to make sure whether I'm on the right track.

Welcome to Stan! It's definitely smart to try it out on toy examples before diving into more complicated models. 

 

This is my Stan code:
data {
  int<lower=1> N;

}
parameters {
  real x[N];
}
model {
for (n in 1:N)
  x[n] ~ normal(0, 1);
}

Yes, that is a correct Stan program if the intent is to have an array x with N >= 1 elements that are independent and each is standard normal. 

It probably isn't essential for this toy example, but you can make it faster by replacing the loop in the model block with the vectorized form x ~ normal(0,1).

load('Ufile.mat');  // Load my initial states which is the vector of ux

N = length(ux);
MCMC_dat = struct('N',N);
initial = struct('x',ux');
fit = stan('file','MC.stan','data',

MCMC_dat,'chains',1,'iter',10,'init',initial,'algorithm','NUTS','sample_file','output-6.csv');
%print(fit)  

I'm not actually familiar with the syntax for specifying these things using MatLabStan (we have a bunch of interfaces to Stan and I haven't actually used the one for MatLab). If you're trying to run 1 Markov chain for 10 iterations (I assume only 1 chain and 10 iter just as a test) with initial values for that 1 chain given by the values in ux then it looks ok to me! But hopefully someone who knows the MatLab interface can let you know if that's actually the right way to specify that.

Best,

Jonah

[Bob Carpenter]

I didn't understand what you meant by

"which i set it as the initial values of my problem for the HMC NUTS sampling."

You can initialize parameters, but you almost never need to in Stan.

The best way to write that program is:

data {
int<lower=0> N;
}
model { }
generated quantities {
real y[N];
for (n in 1:N) y[n] = normal_rng(0, 1);
}

Then each iteration gives you a random draw from a normal.
It's a Monte Carlo approach rather than a Markov chain
Monte Carlo (MCMC) approach, which is always better if you
can pull it off. But it's not going to generalize to cases
where you are fitting data.

- Bob

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Hamed Nikbakht]

Hi Jonah,

Thank you very much for your points.

They were helpful a lot for verification.

Hamed

[Hamed Nikbakht]

It is the part of a technique to estimate the small failure probabilities in high dimensions.

This method expresses probability of failure (Pf) as a product of conditional probabilities that are significantly larger than Pf. These conditional probabilities are estimated by application of Markov Chain Monte Carlo (MCMC) sampling.

These samples are used as “seeds” for generating additional samples conditional on F1 (based on figure) at simulation level 1.

The seeds -which are the initial values as I mentioned before - for generating conditional samples at Level i are discarded after use. This provides some convenience in sample accounting and analysis of the failure
probability estimate. It also reduces the correlation between the samples at different simulation levels.

I attach the figure for more clarification about the technique.

My toy example is very similar to example 4.1.1 250-dimensional multivariate normal (MVN) in the paper of Hoffman & Gelman: The No-U-Turn Sampler:…..

But with this difference that I have the initial values as I mentioned above and my problem's dimension is 1000.

So, based on above explanation, Do you think my Matlab and Stan code is correct?

Thank you very much for your consideration.

Regards,

Hamed

[Hamed Nikbakht]

And i still have the following issue:

'Inference for Stan model: MC_model'
    '1 chains: each with iter=(5); warmup=(0); thin=(1); 5 iterations saved.'
    'Warmup took (0.023) seconds, 0.023 seconds total'
    'Sampling took (0.012) seconds, 0.012 seconds total'


Which is zero warm up and i run it for 10 iterations but it just save 5 iterations.

What should i do for overcoming this problem?

I highly appreciate your help.

Thank you

On Saturday, March 4, 2017 at 4:57:23 PM UTC-5, Bob Carpenter wrote:

[Hamed Nikbakht]

I have another issue about running the Stan in Matlab.

I set the warmup period equal to 10 but it still shows me zero warmup up period.

I would appreciate it if you could help me out in this regard.

I don't know what the problem is.

fit = stan('file','MC.stan','data',MCMC_dat,'chains',1,'iter',20,'warmup',10,'init',initial,'algorithm','NUTS','sample_file','output-6.csv');

Results:

Inference for Stan model: MC_model

1 chains: each with iter=(20); warmup=(0); thin=(1); 20 iterations saved.
Warmup took (0.016) seconds, 0.016 seconds total
Sampling took (0.048) seconds, 0.048 seconds total

Thank you very much.

Hamed


On Friday, March 3, 2017 at 7:26:58 PM UTC-5, Jonah Gabry wrote:

[Bob Carpenter]

It's probably doing the warmup but not saving the iterations. It
does seem to be reporting time for warmup.

I don't know if iter in the MATLAB implementation is the number
of sampling iterations or number of total iterations.

- Bob

[Hamed Nikbakht]

Thank you very much for the explanation.

I have another quick question, when i run it for 10 iterations, Is the last row of "fit.sim.samples.x" for 10th iteration?

(fit.sim.samples.x give me a 10*1000 output. i want to make sure that whether the 10th row is for the last iteration or not.) How does it save the iterations?

Thanks a lot.
Hamed

[Bob Carpenter]

It just calls CmdStan, so you might want to look at the doc for
that package. It saves iterations as a CSV file with header for
parameter names and one row per draw. There's other info tucked
into comments. Just look at the header---it's human readable.

- Bob

[Bob Carpenter]

I meant look at the header for a toy problem (or just
look in the manual).

- Bob

> On Mar 12, 2017, at 7:39 PM, Hamed Nikbakht <hamednik...@gmail.com> wrote:
>

[Hamed Nikbakht]

Thank you very much for your help, Dr.Carpenter.

Regards,
Hamed