# Weibull mixed distribution

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()``````