I am trying to model data which comes in the form of “two cones”. You can imagine this as two cones of different negative slopes and heights positioned on a x,y grid and at every x,y the observed value is the maximum between two cones + some noise (so kind of like maximum of two mountains over a grid). Everything works good when it is just a single cone but when I tried to do double cones it initially did not work very well. Looking at the posteriors revealed the problem. Since there was no ordering constraint, in some chains the first apex coordinate parameter would sample from the first real apex coordinate and in some other chains from the second real apex.

So next what I tried to do was to impose an order on the x variables (so basically the x coordinate of the second cone would be first + positive value). In this case now some of the chains sample as expected (ac0 samples the apex coordinate with smaller x etc) and some of the chains are stuck in some possibly local optima where both apex x coordinates are very close to each other. Is there are more efficient way to re-parametrize this model?

In fact the most efficient way seems to be to run the unordered model but then go through each chain and shuffle the parameters so that they agree with each other. However this seems a bit impractical for large number of cones. So I am wondering if there is way to do it more efficiently by better priors?

```
import pymc as pm
import arviz as az
import numpy as np
import pytensor.tensor as ptt
from scipy.spatial.distance import cdist
def fit_cones(ncones, observed_vals, measurement_coordinates):
with pm.Model():
sd = pm.InverseGamma("sd", 2, 6)
x0 = pm.Uniform("x0", -7, 7, size=(1,))
offsets = pm.Uniform("offsets", 0, 12, size=(ncones-1,))
xcoordinates = ptt.concatenate([x0, x0+offsets])
ycoordinates = pm.Normal("ycoordinates", np.zeros((ncones,)),
10, shape=(ncones, ))
apex_coordinates =\
pm.Deterministic("apex_coordinates",
ptt.concatenate([xcoordinates[:,None],
ycoordinates[:,None]], axis=1))
heights = pm.Uniform("heights", np.zeros((ncones,)),
10*np.ones((ncones,)), shape=(ncones))
slopes = 0.1 + pm.HalfNormal("slopes", np.ones((ncones,)), shape=(ncones))
#difference between apex_coordinates and all the measurement_coordinates
#apex coordinates is ncones x None x ndims
#measurement coordinates is None x npoints x ndims
dif_c = apex_coordinates[:,None,:] - measurement_coordinates[None, :, :]
#distance between apex corodinates and all the measurement_coordinates
#norm_dif_c is ncones x npoints
norm_dif_c = pm.math.sqrt(pm.math.sqr(dif_c[:,:,0])
+ pm.math.sqr(dif_c[:,:,1]))
fitted_y = heights[:, None] - slopes[:, None]*norm_dif_c
#fitted_y is npoints
fitted_y = fitted_y.max(axis=0)
#observed_vals is nrepeats x npoints
pm.Normal("likelihood", mu=fitted_y[None,:], sigma=sd, observed=observed_vals)
trace = pm.sample(draws=2000, tune=1000, chains=6)
return trace
def make_cone(s, h, acs, coordinates):
dists = cdist(acs, coordinates)
cones = h[:,None] - dists*s[:,None]
return np.max(np.array(cones), axis=0)
np.random.seed(139111)
apex = np.array([[0,0],[5,-2]]).T
ncones = 2
npoints = 100
nrepeats = 10
real_sd = 1
measurement_coordinates = np.random.rand(npoints,2)*10 - 5
s = np.random.rand(ncones)*0.25 + 0.75 #slopes
h = np.random.uniform(3,8, (ncones,)) #heights
acs = measurement_coordinates[np.random.randint(0, npoints+1, ncones),:]
real_vals = make_cone(s, h, acs, measurement_coordinates)
observed_vals = real_vals[None,:] + np.random.normal(0, real_sd, size=(nrepeats,npoints))
trace=\
fit_cones(ncones, observed_vals, measurement_coordinates)
```