Optimization and Variable Model Components

[ralphe...@gmail.com]

I've created the code (below) to fit a series of time-varying peaks to a surface, the parameters of which each evolve using random walks. It works very well, however, it runs extremely slowly (5-8 hours on an Nx=50  x Ny=100 matrix).

Is there anything obvious that I can do to speed this up? The exponents in my model are time-varying I can't precompute them, and I can't vectorize the ^ operation. I ultimately want to be able to have more than 3 peaks, and right now the model fit time is impractical for use.

Assuming this can run in a reasonable time, I would also like to be able to determine the optimal number of peaks at each time, such that
       surface(t)=sum(In * Pn(params)) where n=1:number of peaks
meaning that the surface at time t is the sum of several peaks (P1...PN), each of which has an indicator function (Ip) dictating whether they are on or off.

Even better would be letting the indicator functions switch on/off based on some estimated markov probability. I know stan doesn't estimate discrete parameters like this, and I had a few thoughts of how to approach this problem. One way might be to create some transition probability parameter for each component and perform a test against a pseudorandom number at each time to determine if a switch occurs. Alternatively, I could try all possible combinations indicator functions at a given time and choose based on maximum likelihood, which I believe is similar to marginalizing out the discrete parameter.

Any thoughts?

Thanks!!

data {
    int<lower=0>         Nx; //number of xvals
    int<lower=0>         Ny; //number of yvals

    vector[Ny]           y_vals; //y values
    matrix[Nx, Ny]       surface; //surface
}

parameters{
    //Peak 1
    real<lower=8, upper=35>     peak1_Ft[Nx];
    real<lower=0, upper=20>     peak1_At[Nx];
    real<lower=5, upper=100>    peak1_Bt[Nx];

    real<lower=.8, upper=5>     peak1_sig2Ft;
    real<lower=0>               peak1_sig2At;
    real<lower=0>               peak1_sig2Bt;

    //Peak 2
    real<lower=0, upper=1>      peak2_Ft[Nx];
    real<lower=0>               peak2_At[Nx];
    real<lower=0, upper=2>      peak2_Bt[Nx];

    real<lower=0, upper= 1>     peak2_sig2Ft;
    real<lower=0>               peak2_sig2At;
    real<lower=0>               peak2_sig2Bt;

    //Peak 3
    real<lower=1, upper=3>      peak3_Ft[Nx];
    real<lower=0>               peak3_At[Nx];
    real<lower=1, upper=10>     peak3_Bt[Nx];

    real<lower=0, upper= 1>     peak3_sig2Ft;
    real<lower=0>               peak3_sig2At;
    real<lower=0>               peak3_sig2Bt;

    //Baseline
    real<lower=0>               baseline_At[Nx];
    real<lower=0, upper=1>      baseline_Rt[Nx];
    real                        baseline_Ot[Nx];

    real<lower=0>               baseline_sig2At;
    real<lower=0>               baseline_sig2Rt;
    real<lower=0>               baseline_sig2Ot;

    //Observation Noise
    real<lower=0, upper=10>      sig2obs;
}


model {
    //Model variables
    matrix[Nx, Ny] mu; //Holds the mean

    //Gamma parameters
    real peak1_k;
    real peak1_theta;
    real peak2_k;
    real peak2_theta;
    real peak3_k;
    real peak3_theta;
   
    //Priors
    peak1_Ft[1] ~ normal(20,5);
    peak2_Ft[1] ~ normal(.1,.5);
    peak2_Bt[1] ~ normal(.1,.05);
    peak3_Ft[1] ~ normal(1,2);
    peak3_Bt[1] ~ normal(2,.5);

    baseline_Ot[1] ~ normal(-25,5);
    baseline_Rt[1] ~ normal(.01,.001);
    baseline_At[1] ~ normal(50,10);

    sig2obs ~ cauchy(0,.5);    
 
    for (x in 1:Nx){
       if(x>1)
        {
            //Create Random Walk
            peak1_Ft[x] ~ normal(peak1_Ft[x-1], peak1_sig2Ft);
            peak1_At[x] ~ normal(peak1_At[x-1], peak1_sig2At);
            peak1_Bt[x] ~ normal(peak1_Bt[x-1], peak1_sig2Bt);

            peak2_Ft[x] ~ normal(peak2_Ft[x-1], peak2_sig2Ft);
            peak2_At[x] ~ normal(peak2_At[x-1], peak2_sig2At);
            peak2_Bt[x] ~ normal(peak2_Bt[x-1], peak2_sig2Bt);

            peak3_Ft[x] ~ normal(peak3_Ft[x-1], peak3_sig2Ft);
            peak3_At[x] ~ normal(peak3_At[x-1], peak3_sig2At);
            peak3_Bt[x] ~ normal(peak3_Bt[x-1], peak3_sig2Bt);

            baseline_At[x] ~ normal(baseline_At[x-1], baseline_sig2At);
            baseline_Rt[x] ~ normal(baseline_Rt[x-1], baseline_sig2Rt);
            baseline_Ot[x] ~ normal(baseline_Ot[x-1], baseline_sig2Ot);
        }

        //Compute Gamma Terms
        peak2_k=(sqrt(peak2_Ft[x]^2*(4*peak2_Bt[x]+peak2_Ft[x]^2))+2*peak2_Bt[x]+peak2_Ft[x]^2)/(2*peak2_Bt[x]);
        peak2_theta=(sqrt(peak2_Ft[x]^2*(4*peak2_Bt[x]+peak2_Ft[x]^2))-peak2_Ft[x]^2)/(2*peak2_Ft[x]);

        peak3_k=(sqrt(peak3_Ft[x]^2*(4*peak3_Bt[x]+peak3_Ft[x]^2))+2*peak3_Bt[x]+peak3_Ft[x]^2)/(2*peak3_Bt[x]);
        peak3_theta=(sqrt(peak3_Ft[x]^2*(4*peak3_Bt[x]+peak3_Ft[x]^2))-peak3_Ft[x]^2)/(2*peak3_Ft[x]);

        peak1_k=(sqrt(peak1_Ft[x]^2*(4*peak1_Bt[x]+peak1_Ft[x]^2))+2*peak1_Bt[x]+peak1_Ft[x]^2)/(2*peak1_Bt[x]);
        peak1_theta=(sqrt(peak1_Ft[x]^2*(4*peak1_Bt[x]+peak1_Ft[x]^2))-peak1_Ft[x]^2)/(2*peak1_Ft[x]);

        for (y in 1:Ny){
            //Add baseline
            mu[x,y] = baseline_At[x]*(1-baseline_Rt[x])^y_vals[y]+baseline_Ot[x];
   
            //Add Peak 1
            mu[x,y] = mu[x,y] + peak1_At[x]/(exp(1-peak1_k)*((peak1_k-1)*peak1_theta)^(peak1_k-1))*y_vals[y]^(peak1_k-1)*exp(-y_vals[y]/peak1_theta);

            //Add Peak 2
            mu[x,y] = mu[x,y] + peak2_At[x]/(exp(1-peak2_k)*((peak2_k-1)*peak2_theta)^(peak2_k-1))*y_vals[y]^(peak2_k-1)*exp(-y_vals[y]/peak2_theta);
   
            //Add Peak 3
            mu[x,y] = mu[x,y] + peak3_At[x]/(exp(1-peak3_k)*((peak3_k-1)*peak3_theta)^(peak3_k-1))*y_vals[y]^(peak3_k-1)*exp(-y_vals[y]/peak3_theta);
        }
    }
 
      to_vector(surface) ~ normal(to_vector(mu), sig2obs);
}

[Bob Carpenter]

> On Sep 30, 2016, at 10:24 AM, ralphe...@gmail.com wrote:
>
> I've created the code (below) to fit a series of time-varying peaks to a surface, the parameters of which each evolve using random walks. It works very well, however, it runs extremely slowly (5-8 hours on an Nx=50 x Ny=100 matrix).
>
> Is there anything obvious that I can do to speed this up?

Scaling so everything is unit scale can be a huge help, especially
in finding a good adaptation quickly.

> peak1_Ft[1] ~ normal(20,5);
> peak2_Ft[1] ~ normal(.1,.5);
> peak2_Bt[1] ~ normal(.1,.05);
> peak3_Ft[1] ~ normal(1,2);
> peak3_Bt[1] ~ normal(2,.5);
>
> baseline_Ot[1] ~ normal(-25,5);
> baseline_Rt[1] ~ normal(.01,.001);
> baseline_At[1] ~ normal(50,10);

If you can't rescale the data, make these

peak_Ft_raw[1] ~ normal(0, 1);

peak_Ft[1] = 20 + 5 * peak_Ft_raw;

There's also a time-series example in the manual that suggests
replacing something like

> for (x in 1:Nx){
> if(x>1)
> {
> //Create Random Walk
> peak1_Ft[x] ~ normal(peak1_Ft[x-1], peak1_sig2Ft);

with

peak1_Ft[2:Nx] ~ normal(peak1_Ft[1:(Nx - 1)], peak1_sig2Ft);

this should provide a substantial speedup.

I didn't understand the rest of the model, but if there are
any shared terms in all those computations, saving them and
reusing is a big gain.

I'd strongly recommend writing that cut-and-paste formula in
the following as a function to avoid errors (assuming these really
are all the same formula---that's one of the reasons to use
a function---to highlight the repetition):

> peak2_k=(sqrt(peak2_Ft[x]^2*(4*peak2_Bt[x]+peak2_Ft[x]^2))+2*peak2_Bt[x]+peak2_Ft[x]^2)/(2*peak2_Bt[x]);
> peak2_theta=(sqrt(peak2_Ft[x]^2*(4*peak2_Bt[x]+peak2_Ft[x]^2))-peak2_Ft[x]^2)/(2*peak2_Ft[x]);
>
> peak3_k=(sqrt(peak3_Ft[x]^2*(4*peak3_Bt[x]+peak3_Ft[x]^2))+2*peak3_Bt[x]+peak3_Ft[x]^2)/(2*peak3_Bt[x]);
> peak3_theta=(sqrt(peak3_Ft[x]^2*(4*peak3_Bt[x]+peak3_Ft[x]^2))-peak3_Ft[x]^2)/(2*peak3_Ft[x]);
>
> peak1_k=(sqrt(peak1_Ft[x]^2*(4*peak1_Bt[x]+peak1_Ft[x]^2))+2*peak1_Bt[x]+peak1_Ft[x]^2)/(2*peak1_Bt[x]);
> peak1_theta=(sqrt(peak1_Ft[x]^2*(4*peak1_Bt[x]+peak1_Ft[x]^2))-peak1_Ft[x]^2)/(2*peak1_Ft[x]);

same for these:

>
> //Add Peak 1
> mu[x,y] = mu[x,y] + peak1_At[x]/(exp(1-peak1_k)*((peak1_k-1)*peak1_theta)^(peak1_k-1))*y_vals[y]^(peak1_k-1)*exp(-y_vals[y]/peak1_theta);
>
> //Add Peak 2
> mu[x,y] = mu[x,y] + peak2_At[x]/(exp(1-peak2_k)*((peak2_k-1)*peak2_theta)^(peak2_k-1))*y_vals[y]^(peak2_k-1)*exp(-y_vals[y]/peak2_theta);
>
> //Add Peak 3
> mu[x,y] = mu[x,y] + peak3_At[x]/(exp(1-peak3_k)*((peak3_k-1)*peak3_theta)^(peak3_k-1))*y_vals[y]^(peak3_k-1)*exp(-y_vals[y]/peak3_theta);

And in that second one, there's clearly a lot of repeated
computation of 1 - peak1_k (and its negation). Save these and reuse.

> The exponents in my model are time-varying I can't precompute them, and I can't vectorize the ^ operation.

Vectorizing ^ wouldn't add any speed --- loops are fast in Stan.

> I ultimately want to be able to have more than 3 peaks, and right now the model fit time is impractical for use.
>
> Assuming this can run in a reasonable time, I would also like to be able to determine the optimal number of peaks at each time, such that
> surface(t)=sum(In * Pn(params)) where n=1:number of peaks
> meaning that the surface at time t is the sum of several peaks (P1...PN), each of which has an indicator function (Ip) dictating whether they are on or off.

You can also fit models with different numbers of peaks and
compare. Stan's not really set up to do this kind of nonparametrics
(with varying numbers of parameters across MCMC iterations). You
could think about it as a Dirichlet process or something, but at
best you'd be able to approximate that in Stan.

> Even better would be letting the indicator functions switch on/off based on some estimated markov probability. I know stan doesn't estimate discrete parameters like this, and I had a few thoughts of how to approach this problem. One way might be to create some transition probability parameter for each component and perform a test against a pseudorandom number at each time to determine if a switch occurs.

That won't work with Stan's samplers as it introduces
a huge discontinuity that'll block the Hamiltonian Evolution.

> Alternatively, I could try all possible combinations indicator functions at a given time and choose based on maximum likelihood, which I believe is similar to marginalizing out the discrete parameter.

You can't do that choosing within the model for the same reason
as the previous comment.

You can actually just marginalize it out if the combinatorics isn't too bad
(it's 2^N in number of variables).

- 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.

[ralphe...@gmail.com]

Thanks so much! This is very helpful!

If you can't rescale the data, make these

  peak_Ft_raw[1] ~ normal(0, 1);

  peak_Ft[1] = 20 + 5 * peak_Ft_raw;

The data unfortunately can't be rescaled in an obvious way. Is the idea to make every variable unit normal and then fit on scaled version of it so as to make the gradient equal in all dimensions for the MCMC?

Can this be done directly in the model or peak_Ft this now a transformed parameter? I also assume I only have to do this for the first time since peak_Ft[t] is a function of the previous time point.

There's also a time-series example in the manual that suggests
replacing something like

> for (x in 1:Nx){
>        if(x>1)
>         {
>             //Create Random Walk
>             peak1_Ft[x] ~ normal(peak1_Ft[x-1], peak1_sig2Ft);

with

  peak1_Ft[2:Nx] ~ normal(peak1_Ft[1:(Nx - 1)], peak1_sig2Ft);

this should provide a substantial speedup.

Very helpful. I'll give this a try.

I didn't understand the rest of the model, but if there are
any shared terms in all those computations, saving them and
reusing is a big gain.

I'd strongly recommend writing that cut-and-paste formula in
the following as a function to avoid errors (assuming these really
are all the same formula---that's one of the reasons to use
a function---to highlight the repetition):

A good suggestion! Does stan does pass by value or reference? I was afraid that too many function calls might cause additional overhead.

Thanks so much!
Mike

 

[Bob Carpenter]

> On Sep 30, 2016, at 4:19 PM, ralphe...@gmail.com wrote:
>
> Thanks so much! This is very helpful!
>
> If you can't rescale the data, make these
>
> peak_Ft_raw[1] ~ normal(0, 1);
>
> peak_Ft[1] = 20 + 5 * peak_Ft_raw;
>
>
> The data unfortunately can't be rescaled in an obvious way. Is the idea to make every variable unit normal and then fit on scaled version of it so as to make the gradient equal in all dimensions for the MCMC?

Exactly. This helps in a lot of ways.

> Can this be done directly in the model or peak_Ft this now a transformed parameter? I also assume I only have to do this for the first time since peak_Ft[t] is a function of the previous time point.

That would work the same way, by breaking down the normal

foo[t] ~ normal(g(foo[t-1]), sigma);

as

foo_raw[t] ~ normal(0, 1);

foo[t] = g(foo[t - 1]) + sigma * foo_raw[t];

where g is any kind of function you want.

>
> There's also a time-series example in the manual that suggests
> replacing something like
>
> > for (x in 1:Nx){
> > if(x>1)
> > {
> > //Create Random Walk
> > peak1_Ft[x] ~ normal(peak1_Ft[x-1], peak1_sig2Ft);
>
> with
>
> peak1_Ft[2:Nx] ~ normal(peak1_Ft[1:(Nx - 1)], peak1_sig2Ft);
>
> this should provide a substantial speedup.
>
>
> Very helpful. I'll give this a try.
>
> I didn't understand the rest of the model, but if there are
> any shared terms in all those computations, saving them and
> reusing is a big gain.
>
> I'd strongly recommend writing that cut-and-paste formula in
> the following as a function to avoid errors (assuming these really
> are all the same formula---that's one of the reasons to use
> a function---to highlight the repetition):
>

> Does stan does pass by value or reference? I was afraid that too many function calls might cause additional overhead.

By reference. Specifically as (const T&) in C++ terms for
any class argument T like autodiff variables or matrices or
arrays and by value for primitives (double) or (int).

So function calls are basically free compared to everything
else going on.

- Bob

[ralphe...@gmail.com]

That would work the same way, by breaking down the normal

  foo[t] ~ normal(g(foo[t-1]), sigma);

as

  foo_raw[t] ~ normal(0, 1);

  foo[t] = g(foo[t - 1]) + sigma * foo_raw[t];

where g is any kind of function you want.

Since foo_raw is always ~normal(0,1) is it necessary to have T versions of it?  Could I do something like:

parameters {
       real foo_raw;
       real sigma;
}

      real foo[N]//Is foo a transformed parameter?

    //Does this go in the model of transformed parameters?
      foo_raw ~ normal(0,1);
      sigma ~ cauchy(0, .5);

      foo[1] = 20 + sigma*foo_raw;

      foo[2:N]=foo[1:N-1] + sigma*foo_raw;

Is foo now a transformed parameter? If so, to constrain the bounds of foo (upper=a, lower=b), must I now constrain foo_raw to constrain foo indirectly?
 

[Bob Carpenter]

That would be a different model and probably not a very
good one for a time series.

And

On Oct 1, 2016, at 3:51 PM, ralphe...@gmail.com wrote:
>
>
> That would work the same way, by breaking down the normal
>
> foo[t] ~ normal(g(foo[t-1]), sigma);
>
> as
>
> foo_raw[t] ~ normal(0, 1);
>
> foo[t] = g(foo[t - 1]) + sigma * foo_raw[t];
>
> where g is any kind of function you want.
>
>
> Since foo_raw is always ~normal(0,1) is it necessary to have T versions of it?

Yes---it's characterizing the movement of the elements
of the time series on the unit scale. It might be helpful
for you to try it and see the difference. It won't be a time
series model.

> Could I do something like:
>
> parameters {
> real foo_raw;
> real sigma;
> }
>
> real foo[N]//Is foo a transformed parameter?
>
> //Does this go in the model of transformed parameters?
> foo_raw ~ normal(0,1);
> sigma ~ cauchy(0, .5);
>
> foo[1] = 20 + sigma*foo_raw;
>
> foo[2:N]=foo[1:N-1] + sigma*foo_raw;

You could vectorize the calculation of foo[2:N].

>
> Is foo now a transformed parameter?

That or a local variable in the model block.

> If so, to constrain the bounds of foo (upper=a, lower=b), must I now constrain foo_raw to constrain foo indirectly?

Are you sure you really want to do that? Where do the bounds
come from in the generative process you're envisioning here?

If you do need those bounds, you'd need Jacobian adjustments.

- Bob

[ralphe...@gmail.com]

> Since foo_raw is always ~normal(0,1) is it necessary to have T versions of it?  

Yes---it's characterizing the movement of the elements
of the time series on the unit scale.  It might be helpful
for you to try it and see the difference.  It won't be a time
series model.

I guess my confusion was difference between

for t in (2:Nt){
  foo_raw[t] ~ normal(0,1);
  foo[t] = foo[t-1] + sigma*foo_raw[t];
}

and

foo_raw ~ normal(0,1);

for t in (2:Nt){
  foo[t] = foo[t-1] + sigma*foo_raw;
}

I assume the former necessary to do independent draws?

> If so, to constrain the bounds of foo (upper=a, lower=b), must I now constrain foo_raw to constrain foo indirectly?

Are you sure you really want to do that?  Where do the bounds
come from in the generative process you're envisioning here?

If you do need those bounds, you'd need Jacobian adjustments.

For this problem, I'm fitting multiple peaks with time-varying parameters to a surface. In many cases, those peak parameters (such as the center of the peak) have known ranges (from experiments). For example, let's say peaks 1,2, and 3 can have centers that fall between 0 and 10, 20 and 30, and 40 and 50, respectively. My thought was by imposing a hard restriction on the parameters I could keep peaks from crossing out of their ranges (as they time-varied) and to make sure that the parameters in each chain correspond to the same peaks. Even though the parameters evolve, these bounds would have to be imposed at each time point.

Thanks!!!

[Bob Carpenter]

> On Oct 2, 2016, at 12:21 AM, ralphe...@gmail.com wrote:
>
>
> > Since foo_raw is always ~normal(0,1) is it necessary to have T versions of it?
>
> Yes---it's characterizing the movement of the elements
> of the time series on the unit scale. It might be helpful
> for you to try it and see the difference. It won't be a time
> series model.
>
> I guess my confusion was difference between
>
> for t in (2:Nt){
> foo_raw[t] ~ normal(0,1);
> foo[t] = foo[t-1] + sigma*foo_raw[t];
> }

You want to code the above as

for_raw[2:Nt] ~ normal(0, 1);

foo[2:Nt] = foo[1:(Nt-1)] + sigma * foo_raw;

but note that this doesn't set foo[1] or give a distribution to foo_raw[1]
as written.

>
> and
>
> foo_raw ~ normal(0,1);

>
> I assume the former necessary to do independent draws?
>
> > If so, to constrain the bounds of foo (upper=a, lower=b), must I now constrain foo_raw to constrain foo indirectly?

If you must have constraints, that's what you'd want to do.

> Are you sure you really want to do that? Where do the bounds
> come from in the generative process you're envisioning here?
>
> If you do need those bounds, you'd need Jacobian adjustments.
>
> For this problem, I'm fitting multiple peaks with time-varying parameters to a surface. In many cases, those peak parameters (such as the center of the peak) have known ranges (from experiments). For example, let's say peaks 1,2, and 3 can have centers that fall between 0 and 10, 20 and 30, and 40 and 50, respectively. My thought was by imposing a hard restriction on the parameters I could keep peaks from crossing out of their ranges (as they time-varied) and to make sure that the parameters in each chain correspond to the same peaks. Even though the parameters evolve, these bounds would have to be imposed at each time point.

You're probably better off not constraining and letting the data
constrain the shape of the posterior. A data sample is often consistent
with future draws being outside the observed data bounds (e.g., collect ten
people's heights and it's not unlikely the eleventh person will be taller or
shorter than all of the first ten, but as you get a hundred or a thousand people,
it's still possible, but less likely the new one will be the most extreme).

- Bob

[ralphe...@gmail.com]

So I have done my best to implement all of your suggestions, however I keep getting the error:

Exception thrown at line 215: normal_log: Location parameter[2236] is -nan, but must be finite!
If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

without any convergence. Before, I'd get this a few times and then it would converge.

Here is the updated code:

functions {
    //Compute the k parameter
    real get_k(real Ft, real Bt){
        return (sqrt(Ft^2*(4*Bt+Ft^2))+2*Bt+Ft^2)/(2*Bt);
    }

    //Compute the theta parameter
    real get_theta(real Ft, real Bt){
        return (sqrt(Ft^2*(4*Bt+Ft^2))-Ft^2)/(2*Ft);
    }

    //Compute the fixed term in the gamma function
    real get_gamma_term(real At, real k, real theta){
        return At/(exp(1-k)*((k-1)*theta)^(k-1));
    }

    //Compute the gamma function
    real gamma(real gamma_term, real k, real theta, real freq){
        return gamma_term*freq^(k-1)*exp(-freq/theta);
    }
}

data {
    int<lower=0>         Nt; //number of times
    int<lower=0>         Nf; //number of frequencies

    vector[Nf]          yvals; //y values
    matrix[Nt, Nf]         data_matrix; //spectrogram
}

parameters{
    //Peak 1 Oscillation
    vector[Nt]                  peak1_Ft_raw;
    vector[Nt]                  peak1_At_raw;
    vector[Nt]                  peak1_Bt_raw;

    real<lower=0>               peak1_sig2Ft;

    real<lower=0>                 peak1_sig2At;
    real<lower=0>               peak1_sig2Bt;

    //Peak 1 Oscillation
    vector[Nt]                  peak2_Ft_raw;
    vector[Nt]                  peak2_At_raw;
    vector[Nt]                  peak2_Bt_raw;

    real<lower=0>               peak2_sig2Ft;

    real<lower=0>                 peak2_sig2At;
    real<lower=0>               peak2_sig2Bt;

    //Peak 3

    vector[Nt]                  peak3_Ft_raw;
    vector[Nt]                  peak3_At_raw;
    vector[Nt]                  peak3_Bt_raw;

    real<lower=0>               peak3_sig2Ft;

    real<lower=0>                 peak3_sig2At;
    real<lower=0>                 peak3_sig2Bt;

    //Background
    vector[Nt]                  baseline_At_raw;
    vector[Nt]                  baseline_Rt_raw;
    vector[Nt]                  baseline_Ot_raw;

    real<lower=0>               baseline_sig2At;
    real<lower=0>               baseline_sig2Rt;
    real<lower=0>               baseline_sig2Ot;

    //Observation Noise

    real<lower=0>               sig2obs;
}


transformed parameters{
    vector[Nt]                  peak1_Ft;
    vector[Nt]                  peak1_At;
    vector[Nt]                  peak1_Bt;

    vector[Nt]                  peak2_Ft;
    vector[Nt]                  peak2_At;
    vector[Nt]                  peak2_Bt;

    vector[Nt]                  peak3_Ft;
    vector[Nt]                  peak3_At;
    vector[Nt]                  peak3_Bt;

    vector[Nt]                  baseline_At;
    vector[Nt]                  baseline_Rt;
    vector[Nt]                  baseline_Ot;

    //Peak 1 Priors
    peak1_Ft[1] =.1 + peak1_Ft_raw[1]*peak1_sig2Ft;
    peak1_At[1] = 5 + peak1_At_raw[1]*peak1_sig2At;
    peak1_Bt[1] = 1 + peak1_Bt_raw[1]*peak1_sig2Bt;

    //Peak 1 Priors
    peak2_Ft[1] = 1 + peak2_Ft_raw[1]*peak2_sig2Ft;
    peak2_At[1] = 5 + peak2_At_raw[1]*peak2_sig2At;
    peak2_Bt[1] = 1 + peak2_Bt_raw[1]*peak2_sig2Bt;

    //Peak 3 Priors
    peak3_Ft[1] = 20 + peak3_Ft_raw[1]*peak3_sig2Ft;
    peak3_At[1] = 5 + peak3_At_raw[1]*peak3_sig2At;
    peak3_Bt[1] = 1 + peak3_Bt_raw[1]*peak3_sig2Bt;

    //Baseline Priors
    baseline_At[1] =  50 + baseline_At_raw[1]*baseline_sig2At;
    baseline_Rt[1] = .01 + baseline_Rt_raw[1]*baseline_sig2Rt;
    baseline_Ot[1] = -25 + baseline_Ot_raw[1]*baseline_sig2Ot;

    for (x in 2:Nt){
            //Transformed Random Walks
            peak1_Ft[x] =  peak1_Ft[x-1] + peak1_Ft_raw[x]*peak1_sig2Ft;
            peak1_At[x] =  peak1_At[x-1] + peak1_At_raw[x]*peak1_sig2At;
            peak1_Bt[x] =  peak1_Bt[x-1] + peak1_Bt_raw[x]*peak1_sig2Bt;

            peak2_Ft[x] =  peak2_Ft[x-1] + peak2_Ft_raw[x]*peak2_sig2Ft;
            peak2_At[x] =  peak2_At[x-1] + peak2_At_raw[x]*peak2_sig2At;
            peak2_Bt[x] =  peak2_Bt[x-1] + peak2_Bt_raw[x]*peak2_sig2Bt;

            peak3_Ft[x] =  peak3_Ft[x-1] + peak3_Ft_raw[x]*peak3_sig2Ft;
            peak3_At[x] =  peak3_At[x-1] + peak3_At_raw[x]*peak3_sig2At;
            peak3_Bt[x] =  peak3_Bt[x-1] + peak3_Bt_raw[x]*peak3_sig2Bt;

            baseline_At[x] = baseline_At[x-1] + baseline_At_raw[x]*baseline_sig2At;
            baseline_Rt[x] = baseline_Rt[x-1] + baseline_Rt_raw[x]*baseline_sig2Rt;
            baseline_Ot[x] = baseline_Ot[x-1] + baseline_Ot_raw[x]*baseline_sig2Ot;  
    }
}

model {
    //Model variables
    matrix[Nt, Nf] mu; //Holds the mean

    //Gamma parameters
    real peak1_k;
    real peak1_theta;
    real peak2_k;
    real peak2_theta;
    real peak3_k;
    real peak3_theta;

    real peak1_gamma;
    real peak2_gamma;
    real peak3_gamma;
   
    //Keep estimated parameters on the same scale
    peak1_Ft_raw ~ normal(0,1);
    peak1_At_raw ~ normal(0,1);
    peak1_Bt_raw ~ normal(0,1);

    peak2_Ft_raw ~ normal(0,1);
    peak2_At_raw ~ normal(0,1);
    peak2_Bt_raw ~ normal(0,1);

    peak3_Ft_raw ~ normal(0,1);
    peak3_At_raw ~ normal(0,1);
    peak3_Bt_raw ~ normal(0,1);

    baseline_Ot_raw ~ normal(0,1);
    baseline_Rt_raw ~ normal(0,1);
    baseline_At_raw ~ normal(0,1);

    for (x in 1:Nt){
        //Compute Gamma Terms
        peak1_k=get_k(peak1_Ft[x], peak1_Bt[x]);
        peak1_theta=get_theta(peak1_Ft[x], peak1_Bt[x]);
        peak1_gamma=get_gamma_term(peak1_At[x], peak1_k, peak1_theta);

        peak2_k=get_k(peak2_Ft[x], peak2_Bt[x]);
        peak2_theta=get_theta(peak2_Ft[x], peak2_Bt[x]);
        peak2_gamma=get_gamma_term(peak2_At[x], peak2_k, peak2_theta);

        peak3_k=get_k(peak3_Ft[x], peak3_Bt[x]);
        peak3_theta=get_theta(peak3_Ft[x], peak3_Bt[x]);
        peak3_gamma=get_gamma_term(peak3_At[x], peak3_k, peak3_theta);

        for (y in 1:Nf){
            //Add Background
            mu[x,y] = (baseline_At[x])*(1-baseline_Rt[x])^yvals[y]+baseline_Ot[x];
   
            //Add Peak 3
            mu[x,y] = mu[x,y] + gamma(peak1_gamma, peak1_k, peak1_theta, yvals[y]);

            //Add peak1 Oscillation
            mu[x,y] = mu[x,y] + gamma(peak2_gamma, peak2_k, peak2_theta, yvals[y]);
   
            //Add peak2 Oscillation
            mu[x,y] = mu[x,y] + gamma(peak3_gamma, peak3_k, peak3_theta, yvals[y]);
        }
    }
 
      to_vector(data_matrix) ~ normal(to_vector(mu), sig2obs);
}

[Bob Carpenter]

You should see what line 215 is and what might be happening.
One cause would be a local variable that's used before it's defined.

- Bob

> for (t in 2:Nt){
> //Transformed Random Walks
> peak1_Ft[t] = peak1_Ft[t-1] + peak1_Ft_raw[t]*peak1_sig2Ft;
> peak1_At[t] = peak1_At[t-1] + peak1_At_raw[t]*peak1_sig2At;
> peak1_Bt[t] = peak1_Bt[t-1] + peak1_Bt_raw[t]*peak1_sig2Bt;
>
> peak2_Ft[t] = peak2_Ft[t-1] + peak2_Ft_raw[t]*peak2_sig2Ft;
> peak2_At[t] = peak2_At[t-1] + peak2_At_raw[t]*peak2_sig2At;
> peak2_Bt[t] = peak2_Bt[t-1] + peak2_Bt_raw[t]*peak2_sig2Bt;
>
> peak3_Ft[t] = peak3_Ft[t-1] + peak3_Ft_raw[t]*peak3_sig2Ft;
> peak3_At[t] = peak3_At[t-1] + peak3_At_raw[t]*peak3_sig2At;
> peak3_Bt[t] = peak3_Bt[t-1] + peak3_Bt_raw[t]*peak3_sig2Bt;
>
> baseline_At[t] = baseline_At[t-1] + baseline_At_raw[t]*baseline_sig2At;
> baseline_Rt[t] = baseline_Rt[t-1] + baseline_Rt_raw[t]*baseline_sig2Rt;
> baseline_Ot[t] = baseline_Ot[t-1] + baseline_Ot_raw[t]*baseline_sig2Ot;

> //Set up priors for the sigmas
> /*peak1_sig2Ft ~ cauchy(0,.5);
> peak1_sig2At ~ cauchy(0, 1);
> peak1_sig2Bt ~ cauchy(0,.5);
>
> peak2_sig2Ft ~ cauchy(0,.5);
> peak2_sig2At ~ cauchy(0, 1);
> peak2_sig2Bt ~ cauchy(0,.5);
>
> peak3_sig2Ft ~ cauchy(0, 10);
> peak3_sig2At ~ cauchy(0, 1);
> peak3_sig2Bt ~ cauchy(0,.5);
>
> baseline_sig2At ~ cauchy(0, .1);
> baseline_sig2Rt ~ cauchy(0, .01);
> baseline_sig2Ot ~ cauchy(0, .1);
>
> sig2obs ~ cauchy(0,.5); */
>
> for (t in 1:Nt){
> //Compute Gamma Terms
> peak1_k=get_k(peak1_Ft[t], peak1_Bt[t]);
> peak1_theta=get_theta(peak1_Ft[t], peak1_Bt[t]);
> peak1_gamma=get_gamma_term(peak1_At[t], peak1_k, peak1_theta);
>
> peak2_k=get_k(peak2_Ft[t], peak2_Bt[t]);
> peak2_theta=get_theta(peak2_Ft[t], peak2_Bt[t]);
> peak2_gamma=get_gamma_term(peak2_At[t], peak2_k, peak2_theta);
>
> peak3_k=get_k(peak3_Ft[t], peak3_Bt[t]);
> peak3_theta=get_theta(peak3_Ft[t], peak3_Bt[t]);
> peak3_gamma=get_gamma_term(peak3_At[t], peak3_k, peak3_theta);
>
> for (f in 1:Nf){
> //Add Background
> mu[t,f] = (baseline_At[t])*(1-baseline_Rt[t])^freqs[f]+baseline_Ot[t];
>
> //Add Peak 3
> mu[t,f] = mu[t,f] + gamma(peak1_gamma, peak1_k, peak1_theta, freqs[f]);
>
> //Add peak1 Oscillation
> mu[t,f] = mu[t,f] + gamma(peak2_gamma, peak2_k, peak2_theta, freqs[f]);
>
> //Add peak2 Oscillation
> mu[t,f] = mu[t,f] + gamma(peak3_gamma, peak3_k, peak3_theta, freqs[f]);
> }
> }
>
> to_vector(spect) ~ normal(to_vector(mu), sig2obs);

[ralphe...@gmail.com]

215 is the sampling statement line

to_vector(data_matrix) ~ normal(to_vector(mu), sig2obs);

After playing around, this syntax seems to work (simplified to baseline only below)

functions {
  ...
}

data {
...
}

parameters{
    //Background
    vector[Nx]                  baseline_At_raw;
    vector[Nx]                  baseline_Rt_raw;
    vector[Nx]                  baseline_Ot_raw;

    real<lower=0>               baseline_sig2At;
    real<lower=0>               baseline_sig2Rt;
    real<lower=0>               baseline_sig2Ot;

    //Observation Noise
    real<lower=0>               sig2obs;
}


transformed parameters{

    vector[Nx]                  baseline_At;
    vector[Nx]                  baseline_Rt;
    vector[Nx]                  baseline_Ot;

    //Baseline Priors

    baseline_At[1] =  50 + baseline_At_raw[1]*baseline_sig2At;
    baseline_Rt[1] = .01 + baseline_Rt_raw[1]*baseline_sig2Rt;
    baseline_Ot[1] = -25 + baseline_Ot_raw[1]*baseline_sig2Ot;

    for (x in 2:Nx){

            baseline_At[x] = baseline_At[x-1] + baseline_At_raw[x]*baseline_sig2At;
            baseline_Rt[x] = baseline_Rt[x-1] + baseline_Rt_raw[x]*baseline_sig2Rt;
            baseline_Ot[x] = baseline_Ot[x-1] + baseline_Ot_raw[x]*baseline_sig2Ot;  
    }
}

model {
    //Model variables

    matrix[Nx, Ny] mu; //Holds the mean

    //Set up priors for the sigmas
    sig2obs ~ cauchy(0,.5);
   
    bg_sig2At ~ cauchy(0, 1);  
    bg_sig2Rt ~ cauchy(0,.01);  
    bg_sig2Ot ~ cauchy(0, 1);

    //Keep estimated parameters on the same scale
    bg_Ot_raw[1:Nx] ~ normal(0,1);
    bg_Rt_raw[1:Nx] ~ normal(0,1);
    bg_At_raw[1:NX] ~ normal(0,1);

    for (x in 1:Nx){

        for (y in 1:Ny){

            //Add Background
            mu[x,y] = (baseline_At[x])*(1-baseline_Rt[x])^yvals[y]+baseline_Ot[x];
        }
    }

 
      to_vector(data_matrix) ~ normal(to_vector(mu), sig2obs);
}

but fails when I convert the loop with lines such as

baseline_At[x] = baseline_At[x-1] + baseline_At_raw[x]*baseline_sig2At;

into the vectorized form

baseline_At[2:Nx] = baseline_At[1:(Nx-1)] + baseline_At_raw[2:Nx]*baseline_sig2At;

or if add more peak parameters. Could this have to do with the order of the statements, in that the distributions for the sigmas and the raws are declared after the variable transform? It will not let me place those statements in the main body of the model.

Thanks so much for all of your help!! I deeply appreciate it.

[Bob Carpenter]

You didn't say what the error was. It helps to get the full
output. Then it'll tell you if there's an undefined variable,
which would indicate an ordering problem or if some of the types
aren't right (usually need to convert elementwise operations
to regular ones or vice-versa or declare arrays as vectors/matrices).

- Bob

[Ralph Emaurice]

Here is the error, which I get for all 3 parameters:

output-8.csv:   Error evaluating the log probability at the initial value.
output-8.csv: Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
output-8.csv: Exception thrown at line 70: normal_log: Location parameter[3] is nan, but must be finite!
output-8.csv: If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
output-8.csv: but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
output-8.csv: Rejecting initial value:
output-8.csv:   Error evaluating the log probability at the initial value.
output-8.csv: Initialization between (-2, 2) failed after 100 attempts.
output-8.csv:  Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.

Line 70 is the sampling statement:
to_vector(data_matrix) ~ normal(to_vector(mu), sig2obs);
 

 

> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/5jGpRYU0q1A/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

[Daniel Lee]

Sorry -- what do you mean by "I get for all 3 parameters"? Can you just show the whole error message?

This error message is pretty straightforward. The third element of mu is nan. Did you define it? The next step is to just print mu and find out which ones are nan and why.

Daniel

[ralphe...@gmail.com]

During the run, it says the same message over and over but it varying "Location parameter[1]","Location parameter[2]", and "Location parameter[3]".

output-8.csv:   Error evaluating the log probability at the initial value.
output-8.csv: Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

output-8.csv: Exception thrown at line 70: normal_log: Location parameter[1] is -nan, but must be finite!

output-8.csv: If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
output-8.csv: but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
output-8.csv: Rejecting initial value:
output-8.csv:   Error evaluating the log probability at the initial value.

output-8.csv: Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

output-8.csv: Exception thrown at line 70: normal_log: Location parameter[2] is -nan, but must be finite!

output-8.csv: If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
output-8.csv: but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
output-8.csv: Rejecting initial value:
output-8.csv:   Error evaluating the log probability at the initial value.
output-8.csv: Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
output-8.csv: Exception thrown at line 70: normal_log: Location parameter[3] is nan, but must be finite!
output-8.csv: If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
output-8.csv: but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

 Eventually it fails.

All I have done is changed the code (posted earlier in the thread)

//Baseline Priors
    bg_At[1] =  50 + bg_At_raw[1]*bg_sig2At;
    bg_Rt[1] = .01 + bg_Rt_raw[1]*bg_sig2Rt;
    bg_Ot[1] = -25 + bg_Ot_raw[1]*bg_sig2Ot;

    //Transformed Random Walk

   for (t in 2:Nt)
    {

        bg_At[t] = bg_At[t-1] + bg_At_raw[t]*bg_sig2At;
        bg_Rt[t] = bg_Rt[t-1] + bg_Rt_raw[t]*bg_sig2Rt;
        bg_Ot[t] = bg_Ot[t-1] + bg_Ot_raw[t]*bg_sig2Ot;  
    }

which works, but is very slow, to

   //Baseline Priors
    bg_At[1] =  50 + bg_At_raw[1]*bg_sig2At;
    bg_Rt[1] = .01 + bg_Rt_raw[1]*bg_sig2Rt;
    bg_Ot[1] = -25 + bg_Ot_raw[1]*bg_sig2Ot;

    //Transformed Random Walk
    bg_At[2:Nt] = bg_At[1:(Nt-1)] + bg_At_raw[2:Nt]*bg_sig2At;
    bg_Rt[2:Nt] = bg_Rt[1:(Nt-1)] + bg_Rt_raw[2:Nt]*bg_sig2At;
    bg_Ot[2:Nt] = bg_Ot[1:(Nt-1)] + bg_Ot_raw[2:Nt]*bg_sig2At;

And all three variables bg_At, bg_Rt, and bg_Ot quickly go to NaN. Are these programmatically identical?

Also, it seems that putting in a print statement kills my processes. Maybe overloading the IO?

> 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.

--

You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/5jGpRYU0q1A/unsubscribe.

To unsubscribe from this group and all its topics, 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.

--
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.

[Daniel Lee]

Sorry for begin so pedantic, but can you actually just attach both models and the output from both models (run with the same data, the same seed, etc.)? I don't believe that's the only bit you changed. If you look, the index variable has changed from "x" to "t". I'm not sure if that's really the only change or if there's something else that is different. It's really, really hard to help without actually knowing what you've run. (This one isn't a conceptual question so it needs the program in order to help.)

To answer your question, those two should be identical. If they're not, there's a bug in the language -- let's see if you've triggered a bug here.

Regarding the IO... not sure. What interface are you using?

Daniel

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

[Ralph Emaurice]

This is truly the same file (I've been changing some variable names around between posts), and I've just been commenting in and out the two blocks of code on this very simplified version of my model. See attached (model_baseline_only.stan).

I am using MatlabStan for the interface.

To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

[Daniel Lee]

Additional questions:

- what version of CmdStan do you have?

- can you post the output the error messages for both runs with the same seed? You didn't post your data, so it's a program that I really can't run. My guess is that the behavior is seed dependent, which would indicate problems with how you're computing mu. For example, what's freqs?

My guess: your problem is with line 76. (But this is just a guess.)

I don't know what the limitations of MatlabStan are with IO.

Daniel

[ralphe...@gmail.com]

- what version of CmdStan do you have?

I am using 2.12.0

- can you post the output the error messages for both runs with the same seed? You didn't post your data, so it's a program that I really can't run. My guess is that the behavior is seed dependent, which would indicate problems with how you're computing mu. For example, what's freqs?

Ok. Seed = 1234. See vectorized_output.txt (which gives errors then exits) and forloop_output.txt (which takes a gives errors and then quickly converges) attached.

My guess: your problem is with line 76. (But this is just a guess.)      

   
           a*b^X+c, where X>=0  can only be nan or inf if a,b, or c are nan or inf. So I think it could do with the initialization of the variables.

Just to keep the big picture in mind, my goal is to optimize my working (but very slow) model (see surface_model.stan) by changing the variables, which are on different scales (which can't be changed) to transformations on unit normal distributions (see surface_model_opt.stan, which fails).

[Daniel Lee]

Hi -- I'm sorry, but I don't think you're understanding the importance of having the exact files that generate the output. The program you sent earlier doesn't line up with both these outputs. The reason I can tell is that in the forloop_output.txt, the exception happens in line 81 and in the vectorized output it happens in line 80. As far as I know, there's at least one difference that you haven't specified. I really can't tell if there are others and there's an infinite number of ways things could be different.

Space is cheap, so please just send both programs and both outputs exactly. If you attach the data, I could try to run it too.

Regarding your slow model, "optimize" isn't the right term. Vectorizing code and making computation more efficient helps, but only by small percentages. Reparameterizing models and writing models that are better fits to your model will make large jumps in time.

So... if you want more help with this model:

1) send the Stan program with the for loop in one file (.stan)

2) send the Stan program with the slicing in another file (.stan)

3) send the output with the for loop as a text file (like you've done)

4) send the output with the slicing in another file (like you've done)

5) send the data that goes along with the runs

I'm actually trying to help, but there isn't more than can be done when there are lurking issues in the background. If your model is showing a bug, that's really important and I want to be able to hunt it down, but it's really hard when I can't replicate and can't track through enough of your changes.

Daniel

--

[ralphe...@gmail.com]

Sorry for the confusion. I assure you it's the same file. I have both blocks in the same code and commented one in and one out for the different runs.The difference in line number only because of an adjustment in comments. Attached is the exact code with the blocks to comment in or out.

What format would you like the data in? I am doing everything in matlab so I can send a .mat file or convert it if it is easier.

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.

[Daniel Lee]

Please convert to .data.R. Since you're in MatlabStan, it should dump the data out for you. (you might have to search your temp directory or something else, but it should be there.)

And please, please, please, just send two Stan programs. The space is cheap. I'll need to run two programs anyway.

Daniel

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

[ralphe...@gmail.com]

Absolutely. I kept it as one file just to prove to myself that only that one block of code was changing.

Attached:
for loop model: model_baseline_only_forloop.stan
for loop output: forloop_output.txt

vectorized model: model_baseline_only_vectorized.stan
vectorized output: vectorized_output.txt

data:  baseline.data.R (this was temp.data.R, generated by the run)

Many thanks!!!

[Daniel Lee]

Thanks for sending all the files over. It took a direct comparison to actually find the problem.

You have to use the for loop. It's because of this:

bg_At[2:Nx] = bg_At[1:(Nx-1)] + bg_At_raw[2:Nx]*bg_sig2At;

Stan is imperative. Before the loop, you've only defined bg_At[1]. bg_At[2:Nx] is NaN, so when you use bg_At[1:(Nx-1)] elements bg_At[2:(Nx-1)] are NaN and that'll result in the assignment being NaN.

For things like this, Stan isn't slow. Using loops are just as fast as vectorizing it. It's not like Matlab, R, or Python in that aspect. In this case, you have to loop because it's the right thing to do. 

And for future reference, it really helps to send Stan programs with data. Just in case it's useful to anyone else, I debugged this by:

1. Running both programs without changing anything. Verifying the behavior. 

2. Inspecting both programs for consistency. 

3. Reducing the number of iterations. I wanted to make sure the problem occurred with the Stan program.

4. Copying over the working code into the non working code and comparing results. This verifies that everything was the same. 

5. Changing one of the three assignments. It's always easier for me to reason about one change over multiple changes. 

6. Verifying failure with one change. 

7. Digging through the logic for the one change. If I didn't spot it, I would have modified the program to print. 

I couldn't have done any of that with the trickle of changing programs strewn across multiple emails. I was only able to help once you actually sent programs and data over. As an exercise, try going back through your this thread and only use the info you provided to debug. Hopefully you'll see why I kept asking you to send full programs. If everything was attached at the very beginning, we would have been able to find this much sooner. 

And try using print to debug. It's useful. 

Daniel

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

[Avraham Adler]

On Thursday, October 6, 2016 at 1:16:20 AM UTC-4, Daniel Lee wrote:

Thanks for sending all the files over. It took a direct comparison to actually find the problem.

And for future reference, it really helps to send Stan programs with data. Just in case it's useful to anyone else, I debugged this by:

1. Running both programs without changing anything. Verifying the behavior. 

2. Inspecting both programs for consistency. 

3. Reducing the number of iterations. I wanted to make sure the problem occurred with the Stan program.

4. Copying over the working code into the non working code and comparing results. This verifies that everything was the same. 

5. Changing one of the three assignments. It's always easier for me to reason about one change over multiple changes. 

6. Verifying failure with one change. 

7. Digging through the logic for the one change. If I didn't spot it, I would have modified the program to print. 

I couldn't have done any of that with the trickle of changing programs strewn across multiple emails. I was only able to help once you actually sent programs and data over. As an exercise, try going back through your this thread and only use the info you provided to debug. Hopefully you'll see why I kept asking you to send full programs. If everything was attached at the very beginning, we would have been able to find this much sooner. 

And try using print to debug. It's useful. 

Another post that should, in some form or another, be added to the FAQ/blog/Wiki whatever!

Avi

 

[ralphe...@gmail.com]

First off, thank you so much for your help. I really appreciate the deep effort that everyone on this group has shown. It truly means a lot to have this level of support when adopting a new tool, and speaks a great deal to the strength of the community.

Regarding your finding: In the case we just examined, the lack of direct equivalence between the iterative and vectorized forms suggests that the vectorized form of the statement is doing the equivalent of the pass-by-value of the entire vector block (or nonsequental pass-by-reference) at the moment of first execution. Is that correct? That makes it the same as most languages, and makes a great deal of sense.

The reason I tried to vectorize these statements is because it was originally suggested to me (for my non-unit-norm model) that I could change my iterative random walk statements into vectorized statements for efficiency:

There's also a time-series example in the manual that suggests
replacing something like

> for (x in 1:Nx){

>        if(x>1)
>         {
>             //Create Random Walk
>             peak1_Ft[x] ~ normal(peak1_Ft[x-1], peak1_sig2Ft);

with

  peak1_Ft[2:Nx] ~ normal(peak1_Ft[1:(Nx - 1)], peak1_sig2Ft);

this should provide a substantial speedup.

Find attached:

•  Iterative full model (surface_model_iterativeRW.stan)
• The iterative output (iterative_output.txt)
• Vectorized full model (surface_model_vectorizedRW.stan)
• The vectorized output (vectorized_output.txt)
• The data (fullmodel.data.R)

In this case, the vectorization works. While not leading to identical terminal outputs (with the same seed), the vectorized model converged and had almost identical results as the iterative model.

I wonder now, why does the vectorization work in the case of the sampling statement, but not so in the case of the transformed parameters? Do the sampling statements on parameters behave differently then the assignments on the transformed variables?

Thanks so much!

On Thursday, October 6, 2016 at 1:16:20 AM UTC-4, Daniel Lee wrote:

[Bob Carpenter]

Thanks --- we should publicize the Daniel Lee debugging process.

And indeed, the takeaway message is: please send us *everything*
for as small a reproducible example as you can create. And that
means the programs, the driver scripts and all the output that you
got when you tried to run it. Space is cheap.

- Bob

[Bob Carpenter]

C++ is pass by value, but the values can be references (don't
blame me, I just wrestle with C++).

In the C++ translation of a Stan function call, primitives are
passed by pure value, as in (double x) or (int x), whereas any structures
(arrays or matrices) are passed by const reference, as in
(const Eigen::Matrix<T, L, R>& x).

If you do

X = foo(X)

then X is only set after foo(X) is evaluated.

- Bob

> <fullmodel.data.R><surface_model_vectorizedRW.stan><surface_model_iterativeRW.stan><vectorized_output.txt><iterative_output.txt>