Parameterizations

[thodson-usgs]

This issue tracks different parameterizations

[thodson-usgs]

First parameterization using a switch

h = pm.MutableData("h", self.h_obs)

# fixed parameters
# taking the log of h0_offset produces the clipping boundaries in Fig 1, from Reitan et al. 2019
ho = np.ones((self.segments, 1))
ho[0] = 0

# priors
w_mu = np.zeros(self.segments)
# see Le Coz 2014 for default values, but typically between 1.5 and 2.5
w_mu[0] = 1.6
w = pm.Normal("w", mu=w_mu, sigma=0.5, dims="splines")
a = pm.Normal("a", mu=0, sigma=2)

# likelihood
b = pm.Deterministic('b', at.switch(at.le(h, hs), at.log(ho), at.log(h-hs+ho)))
sigma = pm.HalfCauchy("sigma", beta=0.1)
mu = pm.Normal("mu", a + at.dot(w, b), sigma + self.q_sigma, observed=self.y)

[thodson-usgs]

Clip (max) parameterization. The following seem to be equivalent.

b =  pm.Deterministic('b', at.switch(at.gt(h,hs), at.log(h - hs + self.ho), 0
#b = pm.Deterministic('b', at.log( at.clip(h - hs, 0, np.inf) + self.ho)) # b
#b = pm.Deterministic('b', at.log( pm.math.maximum(h - hs, 0) + self.ho)) # b
mu = pm.Normal("mu", a + at.dot(w, b), sigma + self.q_sigma, observed=self.y)

[thodson-usgs]

Now ADVI seems to fail to converge for large sample n>40.

The version in commit a4534369ae294f56a1a7a094b54d60290c83ee31 might give slightly better convergence, but I don't see a difference in model structure.

 b = pm.Deterministic('b', at.log( at.clip(h - hs + self.ho, self.ho, np.inf)))

[thodson-usgs]

The current choice is

b = pm.Deterministic('b', at.log( at.clip(h - hs, 0, np.inf) + self.ho))

which seems equivalent to

      b = pm.Deterministic('b', at.switch(h > hs, at.log(h - hs + self.ho), 0))