Help with pymc3 theano calculation/optimization

Hello all,

I have gotten decent at pymc3 lately but have a heck of a time still with doing calculations within the with loop.

Here goes what I am trying to do, I suspect it is not that hard but I do need some help. Or maybe I am way off base.

I am working to reproduce the calculations in Meeks 2008 that show that:
\frac{\langle (t-d)^n \rangle}{\langle (t-d) \rangle^n} = n!
where t is the time between counts and d is the instrument dead time, and n is a power 2,3,4,5,6…
What I am trying to do is find d where the above inequality holds for all the data t and various powers, n (say 2,3,4).

t is a bunch of points distributed exponential like this:

OK away from the application and on to the question. What I am not sure how to accomplish is building this model to optimize and give a credible interval on d where the equality holds to some fraction (or optimized).

So what I could use examples for are:

  1. Computing the moments ratios within the pymc3 with block
  2. Optimizing the equality to find d.

This scipy code does the job but I don’t know how to do the error bars with is leading me toward pymc3 credible intervals.

def calc_value(inval, D=0.0, M=2):
    return ((np.mean((inval.values-D)**M))/(np.mean((inval.values-D))**M))

def do_full_test(t_bet, minM, maxM, D=0.0):
    ans = pd.DataFrame(columns=['Value', 'Factorial'], dtype=float)
    for M in np.arange(minM, maxM+1):
        ans.loc[M] = (calc_value(t_bet, M=M, D=D), scipy.math.factorial(M))
    ans['Diff'] = ans['Value'] - ans['Factorial']
    return ans.dropna()

ans = do_full_test(t_bet, 1, 6, D=0.1)

t_bet is a pandas dataframe of time between following the above plot.

Any help or hints would be awesome, thanks!