I have been trying to incorporate a multinomial sampling within theano.scan. I have managed to get it work as the theano function, and I think it is possible to write it completely without implementing pymc3. However, I wonder if I can do as a part of pymc3. Then I will be able to sample parameters from different distributions in nicer way.
For example, I have the following program:
import theano import theano.tensor as tt K_sym = tt.iscalar("K_sym") n_sym = tt.lmatrix("n_sym") p_sym = tt.dscalar("p_sym") results_step, updates_step = theano.scan(lambda i, n, p: 2*theano_rng.multinomial(n=n[i,:],pvals=[p,1-p])[:,0], sequences=[tt.arange(K_sym)], n_steps=K_sym, outputs_info=None, non_sequences=[n_sym,p_sym]) compute_step = theano.function(inputs=[K_sym,n_sym,p_sym], outputs=results_step, updates=updates_step) res = compute_step(2,np.transpose([[100,50]]*100),.6) res
and I would like to re-write it as something like
import pymc3 as pm import theano import theano.tensor as tt with pm.Model() as model: # some definition of n, p, and K n = pm.RV(shape=2) #or tt.as_tensor_variable() p = pm.RV() K = 2 results_step, updates_step = theano.scan(lambda i, n, p: 2*theano_rng.multinomial(n=n[i,:],pvals=[p,1-p])[:,0], sequences=[tt.arange(K)], n_steps=K, outputs_info=None, non_sequences=[n,p]) compute_step = theano.function(inputs=[K,n,p], outputs=results_step, updates=updates_step) pm.Deterministic('result',compute_step) sampling = pm.sample_prior_predictive(100)
which is wrong, because
theano.function should not be a part of the model. But I need to have it since there is sampling from multinomial inside
Maybe some of you have had something similar, and would not mind to give some hints or an advice. Thank you in advance!
I have seen Sample within theano.scan?, but that one looked a bit different to me.