Censored pymc for power-law distribution of column density

Questions v5
, ,
7

I am not sure if I understood your suggestion correctly, but here is how I modified my code. But the fit still not looks good. Also, what if I want to include the errors also in my y values then how to do that? I understand the errors on the lower and upper censored would be well defined in that case I will not assign -10 or -999 values to them.

import numpy as np
import matplotlib.pyplot as plt
import pymc as pm
from pymc.distributions.dist_math import normal_lccdf, normal_lcdf
import arviz as az
from copy import copy
import xarray as xr
from scipy.optimize import curve_fit
from scipy.optimize import minimize
import scipy.stats as stats

Define your data

x = np.array([6.7, 14.9, 7.8, 11.3, 12.1, 20.3, 6.7, 23.4, 16.9, 28.8, 30.0, 29.8, 28.6, 8.0, 5.5, 12.6, 15.1, 29.0, 26.6, 15.2, 13.3, 22.6, 3.2, 9.7, 30.9, 33.9, 30.9, 32.5])
y = np.array([14.392, 13.393, 14.493, 14.898, 13.187, 13.512, 14.396, 13.489, 12.652, 13.838, 13.506, 12.796, 12.794, 14.586, 14.204, 13.107, 13.702, 12.483, 13.125, 13.392, 14.314, 14.458, 14.944, 14.777, 11.76328143, 12.30545265, 12.02984059, 11.66659327])
y_err = np.array([0.175, 0.059, 0.011, 0.417, 0.01, 0.024, 0.021, 0.013, 0.044, 0.029, 0.083, 0.049, 0.115, 0.011, 0.023, 0.078, 0.028, 0.098, 0.026, 0.034, 0.038, -1, -1, -1, -999, -999, -999, -999])
#######

id_low = (y_err < -10)
id_up = (y_err > -10) & (y_err < -0)
id_uncensored = (y_err > 0)

defining lower and upper arrays =========

lower = np.zeros(len(x)) - np.inf
lower[id_low] = y[id_low]

upper = np.zeros(len(x)) + np.inf
upper[id_up] = y[id_up]

def censored_regression_new(x, y, lower, upper):

with pm.Model() as model1:
    
    A = pm.Normal("A", mu=0, sigma=1)
    B = pm.Normal("B", mu=0, sigma=1)
    #σ = pm.HalfNormal("σ", sigma=1)

    normal_dist = pm.Normal.dist(mu= A +  x * B)   ###, sigma=σ)
    censored_normal = pm.Censored("censored_normal", normal_dist, lower=lower, upper=upper, observed=y)


with model1:
    trace1 = pm.sample()

return trace1

####==============
censored_fit_new = censored_regression_new(x, y, lower, upper)
az.plot_posterior(censored_fit_new)
####============

def pp_plot(x, y, fit, ax, clr):

##### plot data
ax.scatter(x, y, color="b", marker='*', s=100)




xi = np.arange(0.0,40,0.1)  

post = az.extract(fit)

A_alls = post["A"].values
B_alls = post["B"].values

for ii in range(len(A_alls)):
    y_ppc = A_alls[ii] + xi * B_alls[ii]
    ax.plot(xi, y_ppc, c=clr, alpha=0.05, rasterized=True)

ax.legend()
ax.set(xlabel="x", ylabel="y")

fig, ax = plt.subplots(1, 1, figsize=(10, 5), sharex=True, sharey=True)

pp_plot(x, y, censored_fit_new , ax, ‘gray’)
plt.show()

fit
fit1920×991 64.6 KB