Sorry for the late reply, I have been working on the data fusion for some time using a Bayesian filter. I am not sure if I am using Bayesian updating right, but, I got a function from the following link:
Updating model on PyMC3 with new observed data
However, when it reaches a stage of extracting posterior information from initial data, I keep on getting error talking about an unknown variable. I have modified my code to the following:
def from_posterior(param, samples):
# Assigning linspace variable:
smin, smax = np.min(samples), np.max(samples)
width = smax - smin
x = np.linspace(smin, smax, 100)
# Creating Gaussain KDE:
y = st.gaussian_kde(samples)(x)
# what was never sampled should have a small probability but not 0,
# so we'll extend the domain and use linear approximation of density on it
x = np.concatenate([[x[0] - 3 * width], x, [x[-1] + 3 * width]])
y = np.concatenate([[0], y, [0]])
return Interpolated(param, x, y)
def S2_0_Bayesian_Interface(data):
########################################################################################################################
with pm.Model() as model_0:
# Prior Distributions for unknown model parameters:
nu_0 = pm.HalfNormal('nu_0', sigma=10)
sigma_0 = pm.HalfNormal('sigma_0', sigma=10)
mu_0 = pm.Normal('mu_0', mu=0, sigma=10)
# Observed data is from a Likelihood distributions (Likelihood (sampling distribution) of observations):
observed_data_0 = pm.StudentT('observed_data_0', nu=nu_0, mu=mu_0, sigma=sigma_0, observed=data)
# draw 5000 posterior samples
trace_0 = pm.sample(draws=5000, tune=3000, chains=2, cores=1)
# Obtaining Posterior Predictive Sampling:
post_pred_0 = pm.sample_posterior_predictive(trace_0, samples=3000)
print(post_pred_0['observed_data_0'].shape)
print('\nSummary: ')
print(pm.stats.summary(data=trace_0))
print(pm.stats.summary(data=post_pred_0))
########################################################################################################################
return trace_0, post_pred_0
def S2_P_Bayesian_Interface(nu, mu, sigma, data):
########################################################################################################################
with pm.Model() as model_P:
# Prior Distributions for unknown model parameters:
nu_P = from_posterior('nu_0', trace['nu_0'])
sigma_P = from_posterior('sigma_0', trace['sigma_0'])
mu_P = from_posterior('mu_0', trace['mu_0'])
# Observed data is from a Likelihood distributions (Likelihood (sampling distribution) of observations):
observed_data_P = pm.StudentT('observed_data_P', nu=nu_P, mu=mu_P, sigma=sigma_P, observed=data)
# draw 5000 posterior samples
trace_P = pm.sample(draws=5000, tune=3000, chains=2, cores=1)
# Obtaining Posterior Predictive Sampling:
post_pred_P = pm.sample_posterior_predictive(trace_P, samples=3000)
print(post_pred_P['observed_data_P'].shape)
print('\nSummary: ')
print(pm.stats.summary(data=trace_P))
print(pm.stats.summary(data=post_pred_P))
########################################################################################################################
return trace_P, post_pred_P
def S1_C_Bayesian_Interface(nu, mu, sigma, data):
########################################################################################################################
with pm.Model() as model_C:
# Prior Distributions for unknown model parameters:
nu_C = from_posterior('nu_P', trace['nu_P'])
sigma_C = from_posterior('sigma_P', trace['sigma_P'])
mu_C = from_posterior('mu_P', trace['mu_P'])
# Observed data is from a Likelihood distributions (Likelihood (sampling distribution) of observations):
observed_data_C = pm.StudentT('observed_data_C', nu=nu_C, mu=mu_C, sigma=sigma_C, observed=data)
# draw 5000 posterior samples
trace_C = pm.sample(draws=5000, tune=3000, chains=2, cores=1)
# Obtaining Posterior Predictive Sampling:
post_pred_C = pm.sample_posterior_predictive(trace_C, samples=3000)
print(post_pred_C['observed_data_C'].shape)
print('\nSummary: ')
print(pm.stats.summary(data=trace_C))
print(pm.stats.summary(data=post_pred_C))
########################################################################################################################
return trace_C, post_pred_C
for m in range(len(S2_data_list)):
print('\n')
print('Acessing list number: ', m)
# print(S2_data_list)
# print(len(S2_data_list))
# Initializing the PDF Distribution:
print('#######################################################################################################')
print('Designing the Bayesian PDF for initial Sensor 2:')
trace_S2_0, post_pred_S2_0 = S2_0_Bayesian_Interface(data=S2_data_list[0])
best_dist = getattr(st, 't')
param_S2_0 = [np.mean(trace_S2_0['nu_0']), np.mean(trace_S2_0['mu_0']), np.mean(trace_S2_0['sigma_0'])]
print('Parameters values = ', param_S2_0)
pdf_S2_0, B_x_axis_pdf_S2_0, y_axis_pdf_S2_0 = make_pdf(data=S2_data_list[0], dist=best_dist, params=param_S2_0, size=len(trace_S2_0))
print('#######################################################################################################')
plt.plot(B_x_axis_pdf_S2_0, y_axis_pdf_S2_0, label='%i. step (Initial[S2])' % (m + 1))
# Predict:
print('#######################################################################################################')
print('Designing the Bayesian PDF for predictive Sensor 2:')
trace_S2_P, post_pred_S2_P = S2_P_Bayesian_Interface(nu=np.mean(trace_S2_0['nu_0']), mu=np.mean(trace_S2_0['mu_0']), sigma=np.mean(trace_S2_0['sigma_0']), data=S2_data_list[m])
best_dist = getattr(st, 't')
param_S2_P = [np.mean(trace_S2_P['nu_P']), np.mean(trace_S2_P['mu_P']), np.mean(trace_S2_P['sigma_P'])]
print('Parameters values = ', param_S2_P)
pdf_S2_P, B_x_axis_pdf_S2_P, y_axis_pdf_S2_P = make_pdf(data=S2_data_list[m], dist=best_dist, params=param_S2_P, size=len(trace_S2_P))
print('#######################################################################################################')
plt.plot(B_x_axis_pdf_S2_P, y_axis_pdf_S2_P, label='%i. step (Prediction[S2])' % (m + 1))
# Correct:
print('#######################################################################################################')
print('Designing the Bayesian PDF for correction Sensor 1:')
trace_S1_C, post_pred_S1_C = S1_C_Bayesian_Interface(nu=np.mean(trace_S2_P['nu_P']), mu=np.mean(trace_S2_P['mu_P']), sigma=np.mean(trace_S2_P['sigma_P']), data=S1_data_list[m])
best_dist = getattr(st, 't')
param_S1_C = [np.mean(trace_S1_C['nu_C']), np.mean(trace_S1_C['mu_C']), np.mean(trace_S1_C['sigma_C'])]
print('Parameters values = ', param_S1_C)
pdf_S1_C, B_x_axis_pdf_S1_C, y_axis_pdf_S1_C = make_pdf(data=S1_data_list[m], dist=best_dist, params=param_S1_C, size=len(trace_S1_C))
print('#######################################################################################################')
plt.plot(B_x_axis_pdf_S1_C, y_axis_pdf_S1_C, label='%i. step (Correction[S1])' % (m + 1))
corr_pdf_S1 = y_axis_pdf_S1_C
corr_pdf_S1_list.append(corr_pdf_S1)
print('Corrected PDF: ', corr_pdf_S1)
print('Length of Corrected PDF list = ', len(corr_pdf_S1_list))
x_corr_pdf_list.append(B_x_axis_pdf_S1_C)
When it gets to the predicting stage, I keep on getting:
Unknown variable nu_0
Also, side question:
How can I speed up the sampling for posterior? Because when I run the code sometimes it takes around 20sec (If I remember it correctly 
), but, now when I run the code it takes a long time to obtain posterior distribution like around 300 seconds.