I would not call this a mixture model, but you can model it with PyMC. Getting the MCMC to converge might be a pain, though, because the parameter space is pretty high dimensional.
Let me refactor your data generating procedure to put the data in a DataFrame:
import pymc3 as pm, numpy as np, matplotlib.pyplot as plt, pandas as pd
from scipy import stats
N = 100
mu_1_true = 5
mu_2_true = 10
alpha_true = 5
beta_true = 2
X1 = stats.poisson.rvs(mu_1_true, size=N)
X2 = stats.poisson.rvs(mu_2_true, size=N)
t = stats.beta.rvs(alpha_true, beta_true, size=N)
df = pd.DataFrame()
df['Z1'] = X1*t
df['Z2'] = X2*t
Then one way to try to recover the mu_1 and mu_2 values from the data in PyMC is this:
with pm.Model() as model:
mu1 = pm.Uniform('mu1', lower=0, upper=20)
mu2 = pm.Uniform('mu2', lower=0, upper=20)
alpha = alpha_true
beta = beta_true
t = pm.Beta('t', alpha=alpha, beta=beta, shape=N)
Z1 = pm.Poisson('Z1', mu=mu1*t, observed=df.Z1)
Z2 = pm.Poisson('Z2', mu=mu2*t, observed=df.Z2)
trace = pm.sample(500, cores=2)
The estimates come in a little low when I try it, so maybe there is still room for improvement:
Here is a notebook with all the details in context: https://gist.github.com/aflaxman/b98d3212f59ce24f667225de1fbfafbd