speed up hierarchical model

[inti pedroso]

Hi,

I have implemented the method presented on this, Skelly et al 2011, article using STAN. the authors of the article provided R code implementing a Gibbs sampler but it is a bit slow. It takes 25 min on some small simulated example data. The STAN implementation, in here https://gist.github.com/inti/f8b91ee7bcb98288e314, is faster and it takes a few mins to produce the same results. I compared the results an they are correct. 

I need to run the model on much larger datasets, orders of magnitud larger. Running on a real sized dataset took a day and it did not finish 100 samples, so it does not look very promising.

- My question is, how could I speed up the model?

is it correct to re-write 

for (i in 1:M){
	increment_log_prob(
	log_mix(pi_0,
	beta_log(p[i],a_hat,a_hat) +  beta_log(e[i],1,d_hat),
	beta_log(p[i],f    ,g)     +  beta_log(e[i],1,h)
	)
	);
	}

as following?

	increment_log_prob(
	log_mix(pi_0,
	beta_log(p,a_hat,a_hat) +  beta_log(e,1,d_hat),
	beta_log(p,f    ,g)     +  beta_log(e,1,h)
	)
	);

if correct could it be done similarly with the second increment_log_prob statement

I tried the ADVI and it is certainly much faster, it did run on the full dataset in under one hour. 

- How could I implement stochastic variational approximation on this model? I checked the example on the ADVI arXiv paper but it is not very clear to me yet how the subsampling of the data actually happens.

many thanks in advance!

[Jonah Gabry]

Unfortunately I don't believe log_mix is vectorized. But you can vectorize the beta_binomial that you use later in your model. Instead of

for (i in 1:J) increment_log_prob( beta_binomial_log(Y[i],R[i],alpha[ gene_index[i] ],beta[ gene_index[i] ]) );

you can do something like

model {
  real alphas[J];
  real betas[J];

// ... other stuff in model block

  for (j in 1:J) {
    alphas[j] <- alpha[gene_index[j]];
    betas[j] <- beta[gene_index[j]];
  }
  increment_log_prob(beta_binomial_log(Y, R, alphas, betas)); // or Y ~ beta_binomial(R, alpha, beta)
}

which should be faster than putting the beta_binomial_log within the loop. 

Jonah

[Jonah Gabry]

I made a few typos that I just edited in my reply, but in case you're viewing this over email I meant to write

increment_log_prob(beta_binomial_log(Y, R, alphas, betas));

[Inti pedroso]

Thanks! I will do some timing and I will let you know if it improves.

Inti

On Mar 3, 2016, at 22:20, Jonah Gabry <jga...@gmail.com> wrote:

I made a few typos that I just edited in my reply, but in case you're viewing this over email I meant to write

increment_log_prob(beta_binomial_log(Y, R, alphas, betas));

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

[Inti pedroso]

Hi, you are right, the log_mix cannot be vectorized … it is on the manual. sorry about that. 

having the the second likelihood vectorized makes the code faster, significantly so since J is the largest dimension of the data.

Regards,

Inti

On Mar 3, 2016, at 22:18, Jonah Gabry <jga...@gmail.com> wrote:

Unfortunately I don't believe log_mix is vectorized. But you can vectorize the beta_binomial that you use later in your model. Instead of

for (i in 1:J) increment_log_prob( beta_binomial_log(Y[i],R[i],alpha[ gene_index[i] ],beta[ gene_index[i] ]) );

you can do something like

model {
  real alphas[J];
  real betas[J];

// ... other stuff in model block

  for (j in 1:J) {
    alphas[j] <- alpha[gene_index[j]];
    betas[j] <- beta[gene_index[j]];
  }

  increment_log_prob(beta_binomial_log(Y, R, alpha, beta); // or Y ~ beta_binomial(R, alpha, beta)
}

which should be faster than putting the beta_binomial_log within the loop. 

Jonah

[Jonah Gabry]

On Friday, March 4, 2016 at 12:58:11 PM UTC-5, inti pedroso wrote:

having the the second likelihood vectorized makes the code faster, significantly so since J is the largest dimension of the data.

Great, glad that helped!