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   = 100
D = norm(mu, sd)
droplet_sizes = D.rvs(n_drops,n_columns)

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

# Use elementwise multiplication to zero
# out inactive drops, rather than indexing
# them individually
kill = np.sum(droplet_sizes*active, axis=0)

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('μ', pm.math.sum(droplet_sizes*active_))
    kill_pred = pm.Normal('kill_pred', mu=μ, sd=ϵ, observed=kill)

    trace_b = pm.sample()
                         
pm.traceplot(trace_b);