Second derivative of measured function is PDF of hidden RV

I’m currently struggling with the following model, as I don’t have any good idea how to implement it in PyMC3. In a nutshell the problem looks like this:

There is a hidden time series S_t which probability distribution shall be inferred. This time series cannot be measured directly, but for several values x measurements of a related function H_{x,t} are available. From physical considerations it is known that the second derivative of H_{x,t} is exactly the probability density function of S_t. For the sake of simplicity it can be assumed that the measurements \hat{H}_{x,t} include some normal distributed noise.

S_t \sim \frac{\partial^2 H_{x,t}}{\partial x^2} \Bigg\rvert_{x=S_t} \\ \hat{H}_{x,t} \sim N(H_{x,t},\sigma)

My big problem for the moment is to reformulate the second derivative in a way that it could be implemented in PyMC3 - without making too many assumptions about the underlying process.

Had anyone already some experience with similar models? Or does anyone have some ideas how this problem could be handled? (I actually don’t know all hidden features of PyMC3 & Theano yet… :wink: )

second derivative of H_{x,t} is exactly the probability density function of S_t

I’m not sure if I understood this part correctly, but I’d suggest to think of your problem from a generative perspective:

  1. You assume that there’s some timeseries S_t - how would you model just this part? Gaussian Process, Random Walk…
  2. from this timeseries, you can calculate another (plus unknown constant C_1)
  3. from this timeseries, you can then calculate H_{x,t} (again plus unknown constant C_2)
  4. finally create a likelihood function that compares H_{x,t} with data
1 Like

Thanks for your advise! I guess I had yesterday a similar idea, so I tried to reformulate the model with an integral condition, instead of a differential one. Assuming (at a first glance) a gaussian process as prior for S_t, I obtain the following model:

S_t \sim GP(m(t),k(t,t^\prime)) \\ H_{x,t} = \int_{-\infty}^x(x-s) \, p_{S_t}(s) \, ds \\ \hat{H}_{x,t} \sim N(H_{x,t},\sigma)

For my opinion this is already quite close to the generative structure you have described. However, to implement the integral condition, I have to do two step which are not too intuitive for me:

  1. I need to evaluate directly the probability density function of S_t - Would this be possible via tt.exp(s.distribution.logp(s_)) or are there any caveats? (I actually don’t know if this is correctly evaluated during sampling)
  2. In a naive implementation I would discretize the integral at those values s = \hat{x} where I got measurement of \hat{H}_{x,t} and use the trapezoidal rule to approximate H_{x,t} - Or are there any functionalities in Theano that allow to me to treat this integral symbolically?

If s is a GP, I’m not sure if there’s a logp available that you could use. With any standard distribution it should be fine. But I’m also not sure if I’m following your steps correctly, so I guess just try it out.

For the integral the discretization and a corresponding (cum)sum are the best candidates I can think of at this point.

So, the last days I’ve played around with some experimental implementations and, so far, I guess I found something which is (roughly) doing the job. For the moment I’ve dropped the time-dependency of the problem as it just adds complexity - I think it shouldn’t be a problem to add this dimension later once I’ve found a good implementation.

This is the current code of my model:

import pandas as pd
import matplotlib.pyplot as plt
import pymc3 as pm
import theano.tensor as tt
import arviz as az
import numpy as np

# Excerpt of original dataset for a single point in time
data_x = np.array([2000.,  2000.,  3000.,  3000.,  3000.,  3000.,
        3500.,  3500.,  3500.,  3500.,  3600.,  3600.,  3600.,  3600.,
        3700.,  3700.,  3700.,  3700.,  3800.,  3800.,  3800.,  3800.,
        3900.,  3900.,  3900.,  3900.,  4000.,  4000.,  4000.,  4000.,
        4100.,  4100.,  4100.,  4100.,  4200.,  4200.,  4200.,  4200.,
        4300.,  4300.,  4300.,  4300.,  4400.,  4400.,  4400.,  4400.,
        4500.,  4500.,  4500.,  4500.,  4600.,  4600.,  4600.,  4600.,
        4700.,  4700.,  4700.,  4700.,  4800.,  4800.,  4800.,  4800.,
        4900.,  4900.,  4900.,  4900.,  5000.,  5000.,  5000.,  5000.,
        5100.,  5100.,  5100.,  5100.,  5200.,  5200.,  5200.,  5200.,
        5300.,  5300.,  5300.,  5300.,  5400.,  5400.,  5400.,  5400.,
        5500.,  5500.,  5500.,  5500.,  5600.,  5600.,  5600.,  5600.,
        5700.,  5700.,  5700.,  5700.,  5800.,  5800.,  5800.,  5800.,
        5900.,  5900.,  5900.,  5900.,  6000.,  6000.,  6000.,  6000.,
        6100.,  6100.,  6100.,  6100.,  6200.,  6200.,  6200.,  6200.,
        6300.,  6300.,  6300.,  6300.,  6400.,  6400.,  6400.,  6400.,
        6500.,  6500.,  6500.,  6500.,  6600.,  6600.,  6600.,  6600.,
        6700.,  6700.,  6700.,  6700.,  6800.,  6800.,  6800.,  6800.,
        6900.,  6900.,  6900.,  6900.,  7000.,  7000.,  7000.,  7000.,
        7100.,  7100.,  7100.,  7100.,  7200.,  7200.,  7200.,  7200.,
        7300.,  7300.,  7300.,  7300.,  7400.,  7400.,  7400.,  7400.,
        7500.,  7500.,  7500.,  7500.,  7600.,  7600.,  7600.,  7600.,
        7700.,  7700.,  7700.,  7700.,  7800.,  7800.,  7800.,  7800.,
        7900.,  7900.,  7900.,  7900.,  8000.,  8000.,  8000.,  8000.,
        8100.,  8100.,  8100.,  8100.,  8200.,  8200.,  8200.,  8200.,
        8300.,  8300.,  8300.,  8300.,  8400.,  8400.,  8400.,  8400.,
        8500.,  8500.,  8500.,  8500.,  8600.,  8600.,  8600.,  8600.,
        8700.,  8700.,  8700.,  8700.,  8800.,  8800.,  8800.,  8800.,
        8900.,  8900.,  8900.,  8900.,  9000.,  9000.,  9000.,  9000.,
        9050.,  9050.,  9050.,  9050.,  9100.,  9100.,  9100.,  9100.,
        9150.,  9150.,  9150.,  9150.,  9200.,  9200.,  9200.,  9200.,
        9250.,  9250.,  9250.,  9250.,  9300.,  9300.,  9300.,  9300.,
        9350.,  9350.,  9350.,  9350.,  9400.,  9400.,  9400.,  9400.,
        9450.,  9450.,  9450.,  9450.,  9500.,  9500.,  9500.,  9500.,
        9550.,  9550.,  9550.,  9550.,  9600.,  9600.,  9600.,  9600.,
        9650.,  9650.,  9650.,  9650.,  9700.,  9700.,  9700.,  9700.,
        9750.,  9750.,  9750.,  9750.,  9800.,  9800.,  9800.,  9800.,
        9850.,  9850.,  9850.,  9850.,  9900.,  9900.,  9900.,  9900.,
        9950.,  9950.,  9950.,  9950., 10000., 10000., 10000., 10000.,
       10050., 10050., 10050., 10050., 10100., 10100., 10100., 10100.,
       10150., 10150., 10150., 10150., 10200., 10200., 10200., 10200.,
       10250., 10250., 10250., 10250., 10300., 10300., 10300., 10300.,
       10350., 10350., 10350., 10350., 10400., 10400., 10400., 10400.,
       10450., 10450., 10450., 10450., 10500., 10500., 10500., 10500.,
       10550., 10550., 10550., 10550., 10600., 10600., 10600., 10600.,
       10650., 10650., 10650., 10650., 10700., 10700., 10700., 10700.,
       10750., 10750., 10750., 10750., 10800., 10800., 10800., 10800.,
       10850., 10850., 10850., 10850., 10900., 10900., 10900., 10900.,
       10950., 10950., 10950., 10950., 11000., 11000., 11000., 11000.,
       11100., 11100., 11100., 11100., 11200., 11200., 11200., 11200.,
       11250., 11250., 11250., 11250., 11300., 11300., 11300., 11300.,
       11350., 11350., 11350., 11350., 11400., 11400., 11400., 11400.,
       11450., 11450., 11450., 11450., 11500., 11500., 11500., 11500.,
       11550., 11550., 11550., 11550., 11600., 11600., 11600., 11600.,
       11650., 11650., 11650., 11650., 11700., 11700., 11700., 11700.,
       11750., 11750., 11750., 11750., 11800., 11800., 11800., 11800.,
       11850., 11850., 11850., 11850., 11900., 11900., 11900., 11900.,
       11950., 11950., 11950., 11950., 12000., 12000., 12000., 12000.,
       12050., 12050., 12050., 12050., 12100., 12100., 12100., 12100.,
       12150., 12150., 12150., 12150., 12200., 12200., 12200., 12200.,
       12250., 12250., 12250., 12250., 12300., 12300., 12300., 12300.,
       12350., 12350., 12350., 12350., 12400., 12400., 12400., 12400.,
       12450., 12450., 12450., 12450., 12500., 12500., 12500., 12500.,
       12550., 12550., 12550., 12550., 12600., 12600., 12600., 12600.,
       12650., 12650., 12650., 12650., 12700., 12700., 12700., 12700.,
       12750., 12750., 12750., 12750., 12800., 12800., 12800., 12800.,
       12850., 12850., 12850., 12850., 12900., 12900., 12900., 12900.,
       12950., 12950., 12950., 12950., 13000., 13000., 13000., 13000.,
       13050., 13050., 13050., 13050., 13100., 13100., 13100., 13100.,
       13150., 13150., 13150., 13150., 13200., 13200., 13200., 13200.,
       13250., 13250., 13250., 13250., 13300., 13300., 13300., 13300.,
       13350., 13350., 13350., 13350., 13400., 13400., 13400., 13400.,
       13450., 13450., 13450., 13450., 13500., 13500., 13500., 13500.,
       13550., 13550., 13550., 13550., 13600., 13600., 13600., 13600.,
       13650., 13650., 13650., 13650., 13700., 13700., 13700., 13700.,
       13750., 13750., 13750., 13750., 13800., 13800., 13800., 13800.,
       13850., 13850., 13850., 13850., 13900., 13900., 13900., 13900.,
       13950., 13950., 13950., 13950., 14000., 14000., 14000., 14000.,
       14050., 14050., 14050., 14050., 14100., 14100., 14100., 14100.,
       14150., 14150., 14150., 14150., 14200., 14200., 14200., 14200.,
       14250., 14250., 14250., 14250., 14300., 14300., 14300., 14300.,
       14350., 14350., 14350., 14350., 14400., 14400., 14400., 14400.,
       14450., 14450., 14450., 14450., 14500., 14500., 14500., 14500.,
       14600., 14600., 14600., 14600., 14700., 14700., 14700., 14700.,
       14800., 14800., 14800., 14800., 15000., 15000., 15000., 15000.,
       15500., 15500., 15500., 15500., 16000., 16000., 16000., 16000.,
       16500., 16500., 16500., 16500., 17000., 17000., 17000., 17000.,
       17500., 17500., 17500., 17500., 18000., 18000., 18000., 18000.,
       18500., 18500., 18500., 18500., 19000., 19000., 19000., 19000.,
       19500., 19500., 19500., 19500., 20000., 20000., 20000., 20000.,
       21000., 21000., 21000., 21000., 21500., 21500., 21500., 21500.,
       25000., 25000., 25000., 25000.])
data_h = np.array([2.3000e+00, 2.6000e+00, 6.8000e+00,
       6.8000e+00, 7.0500e+00, 6.4500e+00, 9.4500e+00, 9.0000e+00,
       1.0500e+01, 9.3500e+00, 9.7500e+00, 1.1700e+01, 1.0750e+01,
       1.1200e+01, 1.1700e+01, 1.0650e+01, 1.1600e+01, 1.2550e+01,
       1.1750e+01, 1.2800e+01, 1.3500e+01, 1.3850e+01, 1.4150e+01,
       1.5150e+01, 1.4800e+01, 1.2850e+01, 1.4150e+01, 1.6650e+01,
       1.6000e+01, 1.5450e+01, 1.7600e+01, 1.8050e+01, 1.5300e+01,
       1.6900e+01, 1.9050e+01, 1.8400e+01, 1.9750e+01, 1.6750e+01,
       1.9950e+01, 2.0600e+01, 1.8250e+01, 2.1350e+01, 2.1600e+01,
       2.3200e+01, 2.2450e+01, 1.9900e+01, 2.1800e+01, 2.5150e+01,
       2.4100e+01, 2.3300e+01, 2.7150e+01, 2.6350e+01, 2.5200e+01,
       2.3450e+01, 2.6850e+01, 2.8250e+01, 2.5500e+01, 2.9500e+01,
       3.1750e+01, 2.9500e+01, 2.7850e+01, 3.0750e+01, 3.0000e+01,
       3.3000e+01, 3.4250e+01, 3.1250e+01, 3.2750e+01, 3.4000e+01,
       3.7250e+01, 3.5500e+01, 3.5250e+01, 3.6250e+01, 4.0250e+01,
       3.8500e+01, 3.9500e+01, 4.3250e+01, 4.1500e+01, 3.8000e+01,
       4.1250e+01, 4.6750e+01, 4.2500e+01, 4.4750e+01, 4.8250e+01,
       5.0250e+01, 4.6000e+01, 4.4750e+01, 4.9500e+01, 5.1750e+01,
       4.8750e+01, 5.3750e+01, 5.7750e+01, 5.1750e+01, 5.5750e+01,
       5.3250e+01, 5.9750e+01, 6.1750e+01, 5.7250e+01, 5.5500e+01,
       6.0500e+01, 6.1250e+01, 6.3750e+01, 6.6000e+01, 7.0500e+01,
       6.5750e+01, 6.8250e+01, 6.4500e+01, 6.9250e+01, 7.0750e+01,
       7.5250e+01, 7.2750e+01, 7.2750e+01, 7.7250e+01, 7.5250e+01,
       8.0000e+01, 8.5500e+01, 8.0000e+01, 8.3000e+01, 7.8250e+01,
       8.3500e+01, 9.1000e+01, 8.7750e+01, 8.5750e+01, 8.9000e+01,
       9.4000e+01, 9.6500e+01, 9.1250e+01, 9.7250e+01, 1.0275e+02,
       9.5000e+01, 9.9500e+01, 1.0375e+02, 1.0600e+02, 1.0900e+02,
       1.0125e+02, 1.1025e+02, 1.1225e+02, 1.1575e+02, 1.0775e+02,
       1.1475e+02, 1.2275e+02, 1.1750e+02, 1.1925e+02, 1.2625e+02,
       1.2200e+02, 1.2425e+02, 1.3025e+02, 1.2950e+02, 1.3800e+02,
       1.3200e+02, 1.3400e+02, 1.3800e+02, 1.4625e+02, 1.4225e+02,
       1.4000e+02, 1.5500e+02, 1.4625e+02, 1.4825e+02, 1.5075e+02,
       1.5950e+02, 1.5475e+02, 1.5750e+02, 1.6400e+02, 1.7325e+02,
       1.6850e+02, 1.6400e+02, 1.6650e+02, 1.7850e+02, 1.7575e+02,
       1.7375e+02, 1.8325e+02, 1.8375e+02, 1.9375e+02, 1.8875e+02,
       1.8550e+02, 1.9425e+02, 1.9950e+02, 2.0450e+02, 1.9575e+02,
       2.1575e+02, 2.1050e+02, 2.0700e+02, 2.0525e+02, 2.1675e+02,
       2.1850e+02, 2.2200e+02, 2.2750e+02, 2.3550e+02, 2.2875e+02,
       2.4000e+02, 2.2700e+02, 2.4650e+02, 2.4250e+02, 2.5275e+02,
       2.4100e+02, 2.5450e+02, 2.6000e+02, 2.6600e+02, 2.5550e+02,
       2.8000e+02, 2.6750e+02, 2.6900e+02, 2.7350e+02, 2.8150e+02,
       2.8750e+02, 2.8300e+02, 2.9400e+02, 2.9650e+02, 3.0250e+02,
       2.9800e+02, 3.0950e+02, 3.1800e+02, 3.2550e+02, 3.1200e+02,
       3.1350e+02, 3.3400e+02, 3.2900e+02, 3.4150e+02, 3.2800e+02,
       3.5900e+02, 3.5050e+02, 3.4500e+02, 3.4600e+02, 3.6800e+02,
       3.6350e+02, 3.7650e+02, 3.6250e+02, 3.8600e+02, 3.9550e+02,
       3.8200e+02, 3.8100e+02, 3.9600e+02, 3.9100e+02, 3.9000e+02,
       4.0500e+02, 4.1450e+02, 4.0050e+02, 4.0500e+02, 4.0000e+02,
       4.1050e+02, 4.1500e+02, 4.1000e+02, 4.2500e+02, 4.2050e+02,
       4.3500e+02, 4.2500e+02, 4.2000e+02, 4.3050e+02, 4.4550e+02,
       4.3550e+02, 4.3100e+02, 4.4100e+02, 4.4550e+02, 4.5600e+02,
       4.4100e+02, 4.5200e+02, 4.5650e+02, 4.5200e+02, 4.6700e+02,
       4.6750e+02, 4.6300e+02, 4.7800e+02, 4.6300e+02, 4.7850e+02,
       4.7450e+02, 4.7400e+02, 4.8950e+02, 4.8950e+02, 5.0150e+02,
       4.8600e+02, 4.8650e+02, 4.9750e+02, 5.0150e+02, 5.1350e+02,
       4.9750e+02, 5.0900e+02, 5.2550e+02, 5.1350e+02, 5.0950e+02,
       5.2150e+02, 5.3750e+02, 5.2150e+02, 5.2550e+02, 5.3450e+02,
       5.3500e+02, 5.5050e+02, 5.3800e+02, 5.4750e+02, 5.6350e+02,
       5.5100e+02, 5.4650e+02, 5.6050e+02, 5.6400e+02, 5.7700e+02,
       5.6000e+02, 5.7300e+02, 5.9050e+02, 5.7700e+02, 5.7400e+02,
       5.8800e+02, 5.8650e+02, 6.0450e+02, 5.9100e+02, 6.1850e+02,
       6.0150e+02, 6.0450e+02, 6.0200e+02, 6.3300e+02, 6.1600e+02,
       6.1900e+02, 6.1600e+02, 6.4800e+02, 6.3050e+02, 6.3100e+02,
       6.3350e+02, 6.6300e+02, 6.4600e+02, 6.4600e+02, 6.4850e+02,
       6.6050e+02, 6.6100e+02, 6.6350e+02, 6.7850e+02, 6.7650e+02,
       6.9400e+02, 6.7900e+02, 6.7750e+02, 6.9350e+02, 6.9250e+02,
       6.9450e+02, 7.1000e+02, 7.2700e+02, 7.0900e+02, 7.0850e+02,
       7.1100e+02, 7.2750e+02, 7.2600e+02, 7.4400e+02, 7.2550e+02,
       7.4450e+02, 7.4250e+02, 7.4350e+02, 7.6100e+02, 7.6000e+02,
       7.6200e+02, 7.7850e+02, 7.6100e+02, 7.7900e+02, 7.9650e+02,
       7.7950e+02, 7.7800e+02, 7.9600e+02, 7.9750e+02, 7.9800e+02,
       8.1450e+02, 8.1550e+02, 8.3350e+02, 8.1600e+02, 8.1500e+02,
       8.5300e+02, 8.3400e+02, 8.3500e+02, 8.3500e+02, 8.5500e+02,
       8.5450e+02, 8.7250e+02, 8.5300e+02, 8.9300e+02, 8.7400e+02,
       8.7300e+02, 8.7500e+02, 8.9350e+02, 8.9500e+02, 8.9400e+02,
       9.1300e+02, 9.1450e+02, 9.3400e+02, 9.1450e+02, 9.1600e+02,
       9.3550e+02, 9.5550e+02, 9.3600e+02, 9.3750e+02, 9.5950e+02,
       9.5700e+02, 9.7750e+02, 9.5750e+02, 9.7900e+02, 9.9950e+02,
       9.8150e+02, 9.7950e+02, 1.0280e+03, 1.0255e+03, 1.0240e+03,
       1.0460e+03, 1.0710e+03, 1.0940e+03, 1.0740e+03, 1.0760e+03,
       1.0955e+03, 1.0995e+03, 1.1190e+03, 1.0995e+03, 1.1210e+03,
       1.1270e+03, 1.1245e+03, 1.1445e+03, 1.1505e+03, 1.1455e+03,
       1.1505e+03, 1.1705e+03, 1.1765e+03, 1.1970e+03, 1.1715e+03,
       1.1790e+03, 1.2060e+03, 1.1975e+03, 1.2035e+03, 1.2240e+03,
       1.2315e+03, 1.2510e+03, 1.2330e+03, 1.2250e+03, 1.2610e+03,
       1.2610e+03, 1.2790e+03, 1.2525e+03, 1.2895e+03, 1.2805e+03,
       1.2900e+03, 1.3075e+03, 1.3190e+03, 1.3085e+03, 1.3365e+03,
       1.3180e+03, 1.3670e+03, 1.3480e+03, 1.3375e+03, 1.3480e+03,
       1.3670e+03, 1.3790e+03, 1.3960e+03, 1.3780e+03, 1.4270e+03,
       1.4095e+03, 1.4090e+03, 1.3975e+03, 1.4410e+03, 1.4295e+03,
       1.4580e+03, 1.4410e+03, 1.4730e+03, 1.4910e+03, 1.4595e+03,
       1.4730e+03, 1.5050e+03, 1.5055e+03, 1.4915e+03, 1.5220e+03,
       1.5255e+03, 1.5380e+03, 1.5385e+03, 1.5560e+03, 1.5890e+03,
       1.5585e+03, 1.5720e+03, 1.5720e+03, 1.6230e+03, 1.6060e+03,
       1.5915e+03, 1.6060e+03, 1.6260e+03, 1.6410e+03, 1.6570e+03,
       1.6410e+03, 1.6760e+03, 1.6920e+03, 1.6605e+03, 1.6760e+03,
       1.7280e+03, 1.7120e+03, 1.7120e+03, 1.6955e+03, 1.7640e+03,
       1.7480e+03, 1.7490e+03, 1.7320e+03, 1.8010e+03, 1.7675e+03,
       1.7860e+03, 1.7850e+03, 1.8230e+03, 1.8390e+03, 1.8045e+03,
       1.8230e+03, 1.8760e+03, 1.8430e+03, 1.8630e+03, 1.8610e+03,
       1.9150e+03, 1.9000e+03, 1.9000e+03, 1.8810e+03, 1.9390e+03,
       1.9195e+03, 1.9390e+03, 1.9540e+03, 1.9930e+03, 1.9790e+03,
       1.9585e+03, 1.9805e+03, 2.0205e+03, 2.0330e+03, 2.0180e+03,
       1.9980e+03, 2.0600e+03, 2.0380e+03, 2.0730e+03, 2.0610e+03,
       2.1140e+03, 2.1020e+03, 2.0785e+03, 2.1010e+03, 2.1420e+03,
       2.1550e+03, 2.1420e+03, 2.1195e+03, 2.1610e+03, 2.1840e+03,
       2.1860e+03, 2.1970e+03, 2.2030e+03, 2.2270e+03, 2.2390e+03,
       2.2275e+03, 2.2810e+03, 2.2450e+03, 2.2690e+03, 2.2705e+03,
       2.2880e+03, 2.3140e+03, 2.3130e+03, 2.3240e+03, 2.3310e+03,
       2.3560e+03, 2.3580e+03, 2.3670e+03, 2.4020e+03, 2.4110e+03,
       2.3740e+03, 2.4015e+03, 2.4440e+03, 2.4445e+03, 2.4180e+03,
       2.4550e+03, 2.4890e+03, 2.4990e+03, 2.4895e+03, 2.4615e+03,
       2.5340e+03, 2.5065e+03, 2.5440e+03, 2.5345e+03, 2.5810e+03,
       2.5515e+03, 2.5890e+03, 2.5795e+03, 2.5965e+03, 2.6345e+03,
       2.6250e+03, 2.6245e+03, 2.6715e+03, 2.6810e+03, 2.6425e+03,
       2.6710e+03, 2.6875e+03, 2.7250e+03, 2.7170e+03, 2.7150e+03,
       2.7630e+03, 2.7315e+03, 2.7620e+03, 2.7710e+03, 2.7780e+03,
       2.8100e+03, 2.8080e+03, 2.8175e+03, 2.8645e+03, 2.8550e+03,
       2.8235e+03, 2.8570e+03, 2.8705e+03, 2.9040e+03, 2.9020e+03,
       2.9120e+03, 2.9175e+03, 2.9490e+03, 2.9595e+03, 2.9510e+03,
       3.0080e+03, 3.0000e+03, 2.9640e+03, 2.9970e+03, 3.0440e+03,
       3.0530e+03, 3.0460e+03, 3.0115e+03, 3.0940e+03, 3.0920e+03,
       3.0605e+03, 3.1010e+03, 3.1065e+03, 3.1420e+03, 3.1485e+03,
       3.1400e+03, 3.1545e+03, 3.1900e+03, 3.1955e+03, 3.1880e+03,
       3.2025e+03, 3.2380e+03, 3.2360e+03, 3.2450e+03, 3.2505e+03,
       3.2860e+03, 3.2930e+03, 3.2850e+03, 3.3410e+03, 3.3350e+03,
       3.3330e+03, 3.2985e+03, 3.3475e+03, 3.3830e+03, 3.3815e+03,
       3.3910e+03, 3.4380e+03, 3.4305e+03, 3.3955e+03, 3.4320e+03,
       3.4790e+03, 3.4880e+03, 3.4810e+03, 3.4445e+03, 3.5350e+03,
       3.5280e+03, 3.5290e+03, 3.4925e+03, 3.5840e+03, 3.5440e+03,
       3.5850e+03, 3.5805e+03, 3.6350e+03, 3.6295e+03, 3.5935e+03,
       3.6335e+03, 3.6755e+03, 3.6825e+03, 3.6425e+03, 3.6820e+03,
       3.7350e+03, 3.6915e+03, 3.7330e+03, 3.7245e+03, 3.7905e+03,
       3.8250e+03, 3.8290e+03, 3.8290e+03, 3.9285e+03, 3.8880e+03,
       3.9215e+03, 3.9290e+03, 4.0200e+03, 4.0265e+03, 3.9860e+03,
       4.0260e+03, 4.2250e+03, 4.2265e+03, 4.1845e+03, 4.2195e+03,
       4.7215e+03, 4.7160e+03, 4.6805e+03, 4.7235e+03, 5.2185e+03,
       5.2210e+03, 5.1800e+03, 5.2220e+03, 5.7230e+03, 5.7215e+03,
       5.6800e+03, 5.7180e+03, 6.1815e+03, 6.2155e+03, 6.2225e+03,
       6.2245e+03, 6.7250e+03, 6.7250e+03, 6.7220e+03, 6.6870e+03,
       7.2275e+03, 7.2210e+03, 7.1885e+03, 7.2260e+03, 7.7260e+03,
       7.7205e+03, 7.6890e+03, 7.7290e+03, 8.2295e+03, 8.2205e+03,
       8.2270e+03, 8.1870e+03, 8.7280e+03, 8.7210e+03, 8.6900e+03,
       8.7340e+03, 9.2320e+03, 9.1910e+03, 9.2290e+03, 9.2250e+03,
       1.0232e+04, 1.0193e+04, 1.0224e+04, 1.0238e+04, 1.0693e+04,
       1.0736e+04, 1.0726e+04, 1.0740e+04, 1.4203e+04, 1.4240e+04,
       1.4245e+04, 1.4245e+04])
data_s = np.unique(data_x)

def diff(x):
    dx = np.zeros_like(x,dtype='float')
    dx[:-1] += 0.5*np.diff(x)
    dx[1:] += 0.5*np.diff(x)
    dx[0] *=2
    dx[-1] *=2
    return dx

with pm.Model() as model:
    # Hyperpriors
    mu_s = pm.Cauchy('mu_s',alpha=12000,beta=1000)
    sigma_s = pm.HalfCauchy('sigma_s',beta=10000)
    sigma = pm.HalfNormal('sigma',sigma=1)
    # Hidden RV
    s = pm.Normal('s',mu=mu_s,sigma=sigma_s)
    ps = pm.Deterministic('ps',tt.exp(s.distribution.logp(data_s)))
    h = pm.Deterministic('h',tt.sum(np.maximum(data_x[None,:] - data_s[:,None],0) * diff(data_s)[:,None] * ps[:,None],axis=0))
    #Likelihood
    y = pm.Normal('y',mu=h,sigma=sigma,observed=data_h)
    trace = pm.sample(init='advi+adapt_diag')
    posterior = pm.sample_posterior_predictive(trace)
    idata = az.from_pymc3(trace=trace, posterior_predictive=posterior)

az.plot_trace(idata,var_names=['~ps','~h'])
plt.show()

fig, ax = plt.subplots(3,1,sharex=True,figsize=(16,16))
pm.gp.util.plot_gp_dist(ax[0],trace['ps'],data_s,palette='Blues',plot_samples=True)
ax[0].set_ylabel('ps')
ax[0].grid()
pm.gp.util.plot_gp_dist(ax[1],posterior['y'],data_x,palette='Reds')
ax[1].plot(data_x,data_h,'xk')
ax[1].set_ylabel('y')
ax[1].grid()
pm.gp.util.plot_gp_dist(ax[2],posterior['y'] - np.quantile(posterior['y'],0.5,0),data_x,palette='Reds',plot_samples=False)
ax[2].plot(data_x,(data_h - np.quantile(posterior['y'],0.5,0)),'xk')
ax[2].set_ylabel('y - y_median')
ax[2].set_xlabel('s')
ax[2].grid()
plt.show()

The results are looking reasonable for the moment - However, I’m still quite unsure about these two points:

  1. ps = pm.Deterministic('ps',tt.exp(s.distribution.logp(data_s))): For my understanding I’m evaluating here the PDF of the prior distribution of s, not of the actual distribution of s during sampling. What if I have also some measurements for s in my model - would this line still allow correct inference?
  2. h = pm.Deterministic('h',tt.sum(np.maximum(data_x[None,:]-data_s[:,None],0)*diff(data_s)[:,None]*ps[:,None],axis=0)): This is some quick & dirty approximation of the integral constraint. However, I also found this custom integral operator which should perform better - but so far I was not able to get it running in my model… Or are there any build in functions of theano which allow for symbolic integration?