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:
-
ps = pm.Deterministic('ps',tt.exp(s.distribution.logp(data_s))): For my understanding I’m evaluating here the PDF of the prior distribution ofs, not of the actual distribution ofsduring sampling. What if I have also some measurements forsin my model - would this line still allow correct inference? -
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?