I’m attempting to perform inference on a Weibull mixture using observed data and censored data. While the computation of the mixture’s probability density function (pdf) appears to be functioning correctly, I encounter an error when trying to compile the Weibull lccdfmixt (for censored data) on the survival function (SF).
Below, you’ll find my code along with the specific error message. Your assistance in resolving this matter would be greatly appreciated.
import pymc as pm
import numpy as np
import pytensor
import pytensor.tensor as pt
import arviz as az
import matplotlib.pyplot as plt
import math
# Set values for the prior parameters
gamma_mode = 1
gamma_variance = 0.5
gamma_2_mean = 40.0
gamma_2_variance = 100.0
z = np.array([0.05, 1, 0.18, 0.19])
t = np.array([9.045859001, 9.188807134, 9.188807134 ,9.188807134, 8.74804738, 8.74804738, 8.74804738, 5.525762415, 5.525762415, 5.525762415, 4.455276067, 4.455276067, 4.455276067, 4.455276067])
def weibull_lccdfmixt(x, gammamixt1, gammamixt2, gammamixt3, gammamixt4, p):
return pt.log(p) + ((-x / gammamixt2) ** gammamixt1) + pt.log(1 - p) + ((-x / gammamixt4) ** gammamixt3)
# Definition of prior for gammamixt1
aux1 = 2.0 + (gamma_mode**2.0) / gamma_variance
a = 0.5 * (aux1 + np.sqrt(aux1**2.0 - 4.0))
b = gamma_mode / (a - 1.0)
# Definition of prior for gammamixt2
ap = (gamma_2_mean)**2 / gamma_2_variance
bp = gamma_2_mean / gamma_2_variance
# Definition of prior for gammamixt3
aux2 = 2.0 + (gamma_mode**2.0) / gamma_variance
aa = 0.5 * (aux1 + np.sqrt(aux1**2.0 - 4.0))
bb = gamma_mode / (a - 1.0)
# Definition of prior for gammamixt4
app = (gamma_2_mean)**2 / gamma_2_variance
bpp = gamma_2_mean / gamma_2_variance
# Définition des paramètres du modèle
with pm.Model() as model:
# Paramètre p
p = 0.5
# Paramètres de la distribution Weibull pour PDF1
gammamixt1 = pm.Gamma('gammamixt1', alpha=a, beta=1/b)
gammamixt2 = pm.Gamma('gammamixt2', alpha=ap, beta= bp)
gammamixt3 = pm.Gamma('gammamixt3', alpha=aa, beta=1/bb)
gammamixt4 = pm.Gamma('gammamixt4', alpha=app, beta=bpp)
# Définition des PDF Weibull
components = [
pm.Weibull.dist(alpha=gammamixt1, beta=gammamixt2),
pm.Weibull.dist(alpha=gammamixt3, beta=gammamixt4)
]
# PDF mixte
pdf_mixture = pm.Mixture('pdf_mixture', w=[p, 1 - p], comp_dists=components, observed=z)
y_cens = pm.Potential("y_cens", weibull_lccdfmixt(t, gammamixt1, gammamixt2, gammamixt3, gammamixt4, p))
# Échantillonnage
with model:
trace = pm.sample(10000, tune=1000)
# Résultats
pm.summary(trace)
SamplingError Traceback (most recent call last)
<ipython-input-14-3888d295bcb5> in <module>
63 # Échantillonnage
64 with model:
---> 65 trace = pm.sample(10000, tune=1000)
66
67 # Résultats
~\Anaconda3\lib\site-packages\pymc\sampling\mcmc.py in sample(draws, tune, chains, cores, random_seed, progressbar, step, nuts_sampler, initvals, init, jitter_max_retries, n_init, trace, discard_tuned_samples, compute_convergence_checks, keep_warning_stat, return_inferencedata, idata_kwargs, nuts_sampler_kwargs, callback, mp_ctx, model, **kwargs)
706 ip: Dict[str, np.ndarray]
707 for ip in initial_points:
--> 708 model.check_start_vals(ip)
709 _check_start_shape(model, ip)
710
~\Anaconda3\lib\site-packages\pymc\model.py in check_start_vals(self, start)
1698
1699 if not all(np.isfinite(v) for v in initial_eval.values()):
-> 1700 raise SamplingError(
1701 "Initial evaluation of model at starting point failed!\n"
1702 f"Starting values:\n{elem}\n\n"
SamplingError: Initial evaluation of model at starting point failed!
Starting values:
{'gammamixt1_log__': array(0.21755129), 'gammamixt2_log__': array(3.4161265), 'gammamixt3_log__': array(-0.3050275), 'gammamixt4_log__': array(4.5230101)}
Logp initial evaluation results:
{'gammamixt1': -0.3, 'gammamixt2': -0.08, 'gammamixt3': -0.87, 'gammamixt4': -7.04, 'pdf_mixture': -14.54, 'y_cens': nan}
You can call `model.debug()` for more details.
model.debug()
point={'gammamixt1_log__': array(0.31190536), 'gammamixt2_log__': array(3.68887945), 'gammamixt3_log__': array(0.31190536), 'gammamixt4_log__': array(3.68887945)}
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
<ipython-input-15-57c965fcd523> in <module>
----> 1 model.debug()
~\Anaconda3\lib\site-packages\pymc\model.py in debug(self, point, fn, verbose)
1807 rvs_to_check = list(self.basic_RVs)
1808 if fn in ("logp", "dlogp"):
-> 1809 rvs_to_check += [self.replace_rvs_by_values(p) for p in self.potentials]
1810
1811 found_problem = False
~\Anaconda3\lib\site-packages\pymc\model.py in <listcomp>(.0)
1807 rvs_to_check = list(self.basic_RVs)
1808 if fn in ("logp", "dlogp"):
-> 1809 rvs_to_check += [self.replace_rvs_by_values(p) for p in self.potentials]
1810
1811 found_problem = False
~\Anaconda3\lib\site-packages\pymc\model.py in replace_rvs_by_values(self, graphs, **kwargs)
1537 The graphs in which to perform the replacements.
1538 """
-> 1539 return replace_rvs_by_values(
1540 graphs,
1541 rvs_to_values=self.rvs_to_values,
~\Anaconda3\lib\site-packages\pymc\pytensorf.py in replace_rvs_by_values(graphs, rvs_to_values, rvs_to_transforms, **kwargs)
356 # Clone original graphs so that we don't modify variables in place
357 inputs = [i for i in graph_inputs(graphs) if not isinstance(i, Constant)]
--> 358 equiv = clone_get_equiv(inputs, graphs, False, False, {})
359 graphs = [equiv[n] for n in graphs]
360
~\Anaconda3\lib\site-packages\pytensor\graph\basic.py in clone_get_equiv(inputs, outputs, copy_inputs, copy_orphans, memo, clone_inner_graphs, **kwargs)
1282
1283 # go through the inputs -> outputs graph cloning as we go
-> 1284 for apply in io_toposort(inputs, outputs):
1285 for input in apply.inputs:
1286 if input not in memo:
~\Anaconda3\lib\site-packages\pytensor\graph\basic.py in io_toposort(inputs, outputs, orderings, clients)
1437 # This will change the order returned.
1438 computed = set(inputs)
-> 1439 todo = [o.owner for o in reversed(outputs) if o.owner]
1440 order = []
1441 while todo:
TypeError: object of type 'TensorVariable' has no len()