conc-model

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
from theano import shared
import scipy.stats as stats
from scipy.stats import gamma, norm
import pymc3 as pm
import theano.tensor as tt
import arviz as az

#Make up fake data
mu = 3
sd = 0.5
n_columns = 5
n_drops   = 1000
D = norm(mu, sd)
droplet_sizes = D.rvs(n_drops,n_columns)

## duplicate the same droplet sizes to create second set of data
droplet_sizes = np.tile(droplet_sizes,(2,1))

n = 2. # what I am trying to solve for
active = droplet_sizes > mu + n * sd

kill = np.multiply(droplet_sizes, active)
kill = np.sum(kill, axis=1)
kill = kill.reshape(2,1)

with pm.Model() as model_b:
    #Priors
    tau = pm.Normal('tau', mu = 2, sd=2)
    ϵ = pm.HalfCauchy('ϵ', 5)

    #Observed
    active_ = droplet_sizes > mu + tau * sd
    
    # Likelihood
    μ = pm.Deterministic('μ', (droplet_sizes*active_).sum(axis=1))
   
    kill_pred = pm.Normal('kill_pred', mu=μ, sd=ϵ, observed=kill)

    trace_b = pm.sample(init='advi',cores=5)
                         
pm.traceplot(trace_b);