Repository: pymc-devs/pymc

Issue:
Latex representation broken for models with SymbolicDistributions
```python
with pm.Model() as model:
    x = pm.ZeroInflatedPoisson("x", 0.5, 5)
print(model)
```

```
  File "/home/ricardo/miniconda3/envs/roche_preclinical_hd/lib/python3.9/site-packages/IPython/core/formatters.py", line 224, in catch_format_error
    r = method(self, *args, **kwargs)
  File "/home/ricardo/miniconda3/envs/roche_preclinical_hd/lib/python3.9/site-packages/IPython/core/formatters.py", line 702, in __call__
    printer.pretty(obj)
  File "/home/ricardo/miniconda3/envs/roche_preclinical_hd/lib/python3.9/site-packages/IPython/lib/pretty.py", line 377, in pretty
    return self.type_pprinters[cls](obj, self, cycle)
  File "/home/ricardo/miniconda3/envs/roche_preclinical_hd/lib/python3.9/site-packages/pymc/printing.py", line 185, in _default_repr_pretty
    output = obj.str_repr()
  File "/home/ricardo/miniconda3/envs/roche_preclinical_hd/lib/python3.9/site-packages/pymc/printing.py", line 83, in str_for_model
    names = [s[: s.index("~") - 1] for s in rv_reprs]
  File "/home/ricardo/miniconda3/envs/roche_preclinical_hd/lib/python3.9/site-packages/pymc/printing.py", line 83, in <listcomp>
    names = [s[: s.index("~") - 1] for s in rv_reprs]
ValueError: substring not found
```

Discussion:

Patch:
diff --git a/RELEASE-NOTES.md b/RELEASE-NOTES.md
index 8c896ae66..4aeb873d2 100644
--- a/RELEASE-NOTES.md
+++ b/RELEASE-NOTES.md
@@ -1,152 +1,146 @@
 # Release Notes
 
-<!--
-⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠
-Do not create individual sections for the beta releases.
-Instead update the vNext section until 4.0.0 is out.
-⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠⚠
--->
+## PyMC 4.0.0
 
-## PyMC vNext (4.0.0b1 → ... → 4.0.0b6 → 4.0.0)
-⚠ The changes below are the delta between the releases `v3.11.5` →...→ `v4.0.0`.
+**If you want a description of the highlights of this release, check out the [release announcement](https://www.pymc.io/blog/v4_announcement.html) on our [new website](https://www.pymc.io)**.
+Feel free to read it, print it out, and give it to people on the street -- because _everybody_ has to know PyMC 4.0 is officially out 🍾
 
-### Not-yet working features
-We plan to get these working again, but at this point their inner workings have not been refactored.
-- Timeseries distributions (see [#4642](https://github.com/pymc-devs/pymc/issues/4642))
-- Nested Mixture distributions (see [#5533](https://github.com/pymc-devs/pymc/issues/5533))
-- `pm.sample_posterior_predictive_w` (see [#4807](https://github.com/pymc-devs/pymc/issues/4807))
-- Partially observed Multivariate distributions (see [#5260](https://github.com/pymc-devs/pymc/issues/5260))
+### Do not miss 🚨
 
-Also check out the [milestones](https://github.com/pymc-devs/pymc/milestones) for a potentially more complete list.
-
-### Unexpected breaking changes (action needed)
-+ New API is not available in `v3.11.5`.
-+ Old API does not work in `v4.0.0`.
+  - ⚠️ The project was renamed to "PyMC". Now the library is installed as "pip install pymc" and imported like `import pymc as pm`. See this [migration guide](https://www.pymc-labs.io/blog-posts/the-quickest-migration-guide-ever-from-pymc3-to-pymc-v40/) for more details.
+  - ⚠️ Theano-PyMC has been replaced with Aesara, so all external references to `theano` and `tt` need to be replaced with `aesara` and `at`, respectively (see [4471](https://github.com/pymc-devs/pymc/pull/4471)).
+  - ⚠️ Support for JAX and JAX samplers, also allows sampling on GPUs. [This benchmark](https://www.pymc-labs.io/blog-posts/pymc-stan-benchmark/) shows speed-ups of up to 11x.
+  - ⚠️ Random seeding behavior changed (see [#5787](https://github.com/pymc-devs/pymc/pull/5787))!
+    - Sampling results will differ from those of v3 when passing the same `random_seed` as before. They will be consistent across subsequent v4 releases unless mentioned otherwise.
+    - Sampling functions no longer respect user-specified global seeding! Always pass `random_seed` to ensure reproducible behavior.
+    - `random_seed` now accepts RandomState and Generators besides integers.
+  - ⚠️ The GLM submodule was removed, please use [Bambi](https://bambinos.github.io/bambi/) instead.
+  - ⚠️ PyMC now requires SciPy version `>= 1.4.1` (see [#4857](https://github.com/pymc-devs/pymc/pull/4857)).
+
+#### v3 features not yet working in v4 ⏳
+
+⚠️  We plan to get these working again, but at this point their inner workings have not been refactored.
+  - MvNormalRandomWalk, MvStudentTRandomWalk, GARCH11 and EulerMaruyama distributions (see [#4642](https://github.com/pymc-devs/pymc/issues/4642))
+  - Nested Mixture distributions (see [#5533](https://github.com/pymc-devs/pymc/issues/5533))
+  - `pm.sample_posterior_predictive_w` (see [#4807](https://github.com/pymc-devs/pymc/issues/4807))
+  - Partially observed Multivariate distributions (see [#5260](https://github.com/pymc-devs/pymc/issues/5260))
+
+### New features 🥳
+
+- **Distributions**:
+  - Univariate censored distributions are now available via `pm.Censored`. [#5169](https://github.com/pymc-devs/pymc/pull/5169)
+  - The `CAR` distribution has been added to allow for use of conditional autoregressions which often are used in spatial and network models.
+  - Added a `logcdf` implementation for the Kumaraswamy distribution (see [#4706](https://github.com/pymc-devs/pymc/pull/4706)).
+  - The `OrderedMultinomial` distribution has been added for use on ordinal data which are _aggregated_ by trial, like multinomial observations, whereas `OrderedLogistic` only accepts ordinal data in a _disaggregated_ format, like categorical observations (see [#4773](https://github.com/pymc-devs/pymc/pull/4773)).
+  - The `Polya-Gamma` distribution has been added (see [#4531](https://github.com/pymc-devs/pymc/pull/4531)). To make use of this distribution, the [`polyagamma>=1.3.1`](https://pypi.org/project/polyagamma/) library must be installed and available in the user's environment.
+  - `pm.DensityDist` can now accept an optional `logcdf` keyword argument to pass in a function to compute the cummulative density function of the distribution (see [5026](https://github.com/pymc-devs/pymc/pull/5026)).
+  - `pm.DensityDist` can now accept an optional `moment` keyword argument to pass in a function to compute the moment of the distribution (see [5026](https://github.com/pymc-devs/pymc/pull/5026)).
+
+  - Added an alternative parametrization, `logit_p` to `pm.Binomial` and `pm.Categorical` distributions (see [5637](https://github.com/pymc-devs/pymc/pull/5637)).
+
+- **Model dimensions**:
+
+  - The dimensionality of model variables can now be parametrized through either of `shape` or `dims` (see [#4696](https://github.com/pymc-devs/pymc/pull/4696)):
+    - With `shape` the length of dimensions must be given numerically or as scalar Aesara `Variables`. Numeric entries in `shape` restrict the model variable to the exact length and re-sizing is no longer possible.
+    - `dims` keeps model variables re-sizeable (for example through `pm.Data`) and leads to well defined coordinates in `InferenceData` objects.
+    - An `Ellipsis` (`...`) in the last position of `shape` or `dims` can be used as short-hand notation for implied dimensions.
+  - New features for `pm.Data` containers:
+    - With `pm.Data(..., mutable=False)`, or by using `pm.ConstantData()` one can now create `TensorConstant` data variables. These can be more performant and compatible in situations where a variable doesn't need to be changed via `pm.set_data()`. See [#5295](https://github.com/pymc-devs/pymc/pull/5295). If you do need to change the variable, use `pm.Data(..., mutable=True)`, or `pm.MutableData()`.
+    - New named dimensions can be introduced to the model via `pm.Data(..., dims=...)`. For mutable data variables (see above) the lengths of these dimensions are symbolic, so they can be re-sized via `pm.set_data()`.
+    - `pm.Data` now passes additional kwargs to `aesara.shared`/`at.as_tensor`. [#5098](https://github.com/pymc-devs/pymc/pull/5098).
+  - The length of `dims` in the model is now tracked symbolically through `Model.dim_lengths` (see [#4625](https://github.com/pymc-devs/pymc/pull/4625)).
+
+- **Sampling**:
+  - A small change to the mass matrix tuning methods jitter+adapt_diag (the default) and adapt_diag improves performance early on during tuning for some models. [#5004](https://github.com/pymc-devs/pymc/pull/5004)
+  - New experimental mass matrix tuning method jitter+adapt_diag_grad. [#5004](https://github.com/pymc-devs/pymc/pull/5004)
+  - Support for samplers written in JAX:
+    - Adding support for [numpyro](https://github.com/pyro-ppl/numpyro)'s NUTS sampler via `pymc.sampling_jax.sample_numpyro_nuts()`
+    - Adding support for [blackjax](https://github.com/blackjax-devs/blackjax)'s NUTS sampler via `pymc.sampling_jax.sample_blackjax_nuts()` (see [#5477](ihttps://github.com/pymc-devs/pymc/pull/5477))
+    - `pymc.sampling_jax` samplers support `log_likelihood`, `observed_data`, and `sample_stats` in returned `InferenceData` object (see [#5189](https://github.com/pymc-devs/pymc/pull/5189))
+    - Adding support for `pm.Deterministic` in `pymc.sampling_jax` (see [#5182](https://github.com/pymc-devs/pymc/pull/5182))
+
+- **Miscellaneous**:
+  - The new `pm.find_constrained_prior` function can be used to find optimized prior parameters of a distribution under some
+  constraints (e.g lower and upper bound). See [#5231](https://github.com/pymc-devs/pymc/pull/5231).
+  - Nested models now inherit the parent model's coordinates. [#5344](https://github.com/pymc-devs/pymc/pull/5344)
+  - `softmax` and `log_softmax` functions added to `math` module (see [#5279](https://github.com/pymc-devs/pymc/pull/5279)).
+  - Added the low level `compile_forward_sampling_function` method to compile the aesara function responsible for generating forward samples (see [#5759](https://github.com/pymc-devs/pymc/pull/5759)).
 
-All of the above apply to:
+### Expected breaking changes 💔
 
-- ⚠ Random seeding behavior changed (see [#5787](https://github.com/pymc-devs/pymc/pull/5787))!
-    - Sampling results will differ from those of V3 when passing the same `random_seed` as before. They will be consistent across subsequent V4 releases unless mentioned otherwise.
-    - Sampling functions no longer respect user-specified global seeding! Always pass `random_seed` to ensure reproducible behavior.
-- Signature and default parameters changed for several distributions:
+- `pm.sample(return_inferencedata=True)` is now the default (see [#4744](https://github.com/pymc-devs/pymc/pull/4744)).
+- ArviZ `plots` and `stats` *wrappers* were removed. The functions are now just available by their original names (see [#4549](https://github.com/pymc-devs/pymc/pull/4471) and `3.11.2` release notes).
+- `pm.sample_posterior_predictive(vars=...)` kwarg was removed in favor of `var_names` (see [#4343](https://github.com/pymc-devs/pymc/pull/4343)).
+- `ElemwiseCategorical` step method was removed (see [#4701](https://github.com/pymc-devs/pymc/pull/4701))
+- `LKJCholeskyCov`'s `compute_corr` keyword argument is now set to `True` by default (see[#5382](https://github.com/pymc-devs/pymc/pull/5382))
+-  Alternative `sd` keyword argument has been removed from all distributions. `sigma` should be used instead (see [#5583](https://github.com/pymc-devs/pymc/pull/5583)).
+
+_Read on if you're a developer. Or curious. Or both._
+
+### Unexpected breaking changes (action needed) 😲
+
+#### Very important ⚠️
+  - `pm.Bound` interface no longer accepts a callable class as argument, instead it requires an instantiated distribution (created via the `.dist()` API) to be passed as an argument. In addition, Bound no longer returns a class instance but works as a normal PyMC distribution. Finally, it is no longer possible to do predictive random sampling from Bounded variables. Please, consult the new documentation for details on how to use Bounded variables (see [4815](https://github.com/pymc-devs/pymc/pull/4815)).
+  - BART has received various updates ([5091](https://github.com/pymc-devs/pymc3/pull/5091), [5177](https://github.com/pymc-devs/pymc3/pull/5177), [5229](https://github.com/pymc-devs/pymc3/pull/5229), [4914](https://github.com/pymc-devs/pymc3/pull/4914)) but was removed from the main package in [#5566](https://github.com/pymc-devs/pymc/pull/5566). It is now available from [pymc-experimental](https://github.com/pymc-devs/pymc-experimental).
+  - Removed `AR1`. `AR` of order 1 should be used instead. (see [5734](https://github.com/pymc-devs/pymc/pull/5734)).
+  - The `pm.EllipticalSlice` sampler was removed (see [#5756](https://github.com/pymc-devs/pymc/issues/5756)).
+  - `BaseStochasticGradient` was removed (see [#5630](https://github.com/pymc-devs/pymc/pull/5630))
+  - `pm.Distribution(...).logp(x)` is now `pm.logp(pm.Distribution(...), x)`.
+  - `pm.Distribution(...).logcdf(x)` is now `pm.logcdf(pm.Distribution(...), x)`.
+  - `pm.Distribution(...).random(size=x)` is now `pm.draw(pm.Distribution(...), draws=x)`.
+  - `pm.draw_values(...)` and `pm.generate_samples(...)` were removed.
+  - `pm.fast_sample_posterior_predictive` was removed.
+  - `pm.sample_prior_predictive`, `pm.sample_posterior_predictive` and `pm.sample_posterior_predictive_w` now return an `InferenceData` object by default, instead of a dictionary (see [#5073](https://github.com/pymc-devs/pymc/pull/5073)).
+  - `pm.sample_prior_predictive` no longer returns transformed variable values by default. Pass them by name in `var_names` if you want to obtain these draws (see [4769](https://github.com/pymc-devs/pymc/pull/4769)).
+  - `pm.sample(trace=...)` no longer accepts `MultiTrace` or `len(.) > 0` traces ([see 5019#](https://github.com/pymc-devs/pymc/pull/5019)).
+  - Setting of initial values:
+    - Setting initial values through `pm.Distribution(testval=...)` is now `pm.Distribution(initval=...)`.
+    - `Model.update_start_values(...)` was removed. Initial values can be set in the `Model.initial_values` dictionary directly.
+    - Test values can no longer be set through `pm.Distribution(testval=...)` and must be assigned manually.
+  - `transforms` module is no longer accessible at the root level. It is accessible at `pymc.distributions.transforms` (see[#5347](https://github.com/pymc-devs/pymc/pull/5347)).
+  - `logp`, `dlogp`, and `d2logp` and `nojac` variations were removed. Use `Model.compile_logp`, `compile_dlgop` and `compile_d2logp` with `jacobian` keyword instead.
+  - `pm.DensityDist` no longer accepts the `logp` as its first position argument. It is now an optional keyword argument. If you pass a callable as the first positional argument, a `TypeError` will be raised (see [5026](https://github.com/pymc-devs/pymc/pull/5026)).
+  - `pm.DensityDist` now accepts distribution parameters as positional arguments. Passing them as a dictionary in the `observed` keyword argument is no longer supported and will raise an error (see [5026](https://github.com/pymc-devs/pymc/pull/5026)).
+  - The signature of the `logp` and `random` functions that can be passed into a `pm.DensityDist` has been changed (see [5026](https://github.com/pymc-devs/pymc/pull/5026)).
+
+#### Important:
+  - Signature and default parameters changed for several distributions:
     - `pm.StudentT` now requires either `sigma` or `lam` as kwarg (see [#5628](https://github.com/pymc-devs/pymc/pull/5628))
     - `pm.StudentT` now requires `nu` to be specified (no longer defaults to 1) (see [#5628](https://github.com/pymc-devs/pymc/pull/5628))
     - `pm.AsymmetricLaplace` positional arguments re-ordered (see [#5628](https://github.com/pymc-devs/pymc/pull/5628))
     - `pm.AsymmetricLaplace` now requires `mu` to be specified (no longer defaults to 0) (see [#5628](https://github.com/pymc-devs/pymc/pull/5628))
     - `ZeroInflatedPoisson` `theta` parameter was renamed to `mu` (see [#5584](https://github.com/pymc-devs/pymc/pull/5584)).
-    - `pm.GaussianRandomWalk` initial distribution defaults to unit normal instead of flat (see[#5779](https://github.com/pymc-devs/pymc/pull/5779))
+    - `pm.GaussianRandomWalk` initial distribution defaults to zero-centered normal with sigma=100 instead of flat (see[#5779](https://github.com/pymc-devs/pymc/pull/5779))
     - `pm.AR` initial distribution defaults to unit normal instead of flat (see[#5779](https://github.com/pymc-devs/pymc/pull/5779))
-- BART was removed [#5566](https://github.com/pymc-devs/pymc/pull/5566). It is now available from [pymc-experimental](https://github.com/pymc-devs/pymc-experimental)
-- The `pm.EllipticalSlice` sampler was removed (see [#5756](https://github.com/pymc-devs/pymc/issues/5756)).
-- `BaseStochasticGradient` was removed (see [#5630](https://github.com/pymc-devs/pymc/pull/5630))
-- ⚠ The library is now named, installed and imported as "pymc". For example: `pip install pymc`.
-- ⚠ Theano-PyMC has been replaced with Aesara, so all external references to `theano`, `tt`, and `pymc3.theanof` need to be replaced with `aesara`, `at`, and `pymc.aesaraf` (see [4471](https://github.com/pymc-devs/pymc/pull/4471)).
-- `pm.Distribution(...).logp(x)` is now `pm.logp(pm.Distribution(...), x)`.
-- `pm.Distribution(...).logcdf(x)` is now `pm.logcdf(pm.Distribution(...), x)`.
-- `pm.Distribution(...).random(size=x)` is now `pm.draw(pm.Distribution(...), draws=x)`.
-- `pm.draw_values(...)` and `pm.generate_samples(...)` were removed.
-- `pm.fast_sample_posterior_predictive` was removed.
-- `pm.sample_prior_predictive`, `pm.sample_posterior_predictive` and `pm.sample_posterior_predictive_w` now return an `InferenceData` object by default, instead of a dictionary (see [#5073](https://github.com/pymc-devs/pymc/pull/5073)).
-- `pm.sample_prior_predictive` no longer returns transformed variable values by default. Pass them by name in `var_names` if you want to obtain these draws (see [4769](https://github.com/pymc-devs/pymc/pull/4769)).
-- `pm.sample(trace=...)` no longer accepts `MultiTrace` or `len(.) > 0` traces ([see 5019#](https://github.com/pymc-devs/pymc/pull/5019)).
-- `transforms` module is no longer accessible ta the root level. It is accessible at `pymc.distributions.transforms` (see[#5347](https://github.com/pymc-devs/pymc/pull/5347)).
-- `logpt`, `logpt_sum`, `logp_elemwiset` and `nojac` variations were removed. Use `Model.logpt(jacobian=True/False, sum=True/False)` instead.
-- `dlogp_nojact` and `d2logp_nojact` were removed. Use `Model.dlogpt` and `d2logpt` with `jacobian=False` instead.
-- `logp`, `dlogp`, and `d2logp` and `nojac` variations were removed. Use `Model.compile_logp`, `compile_dlgop` and `compile_d2logp` with `jacobian` keyword instead.
-- `model.makefn` is now called `Model.compile_fn`, and `model.fn` was removed.
-- Methods starting with `fast_*`, such as `Model.fast_logp`, were removed. Same applies to `PointFunc` classes
-- The GLM submodule was removed, please use [Bambi](https://bambinos.github.io/bambi/) instead.
-- `pm.Bound` interface no longer accepts a callable class as argument, instead it requires an instantiated distribution (created via the `.dist()` API) to be passed as an argument. In addition, Bound no longer returns a class instance but works as a normal PyMC distribution. Finally, it is no longer possible to do predictive random sampling from Bounded variables. Please, consult the new documentation for details on how to use Bounded variables (see [4815](https://github.com/pymc-devs/pymc/pull/4815)).
-- `Model(model=...)` kwarg was removed
-- `Model(theano_config=...)` kwarg was removed
-- `Model.size` property was removed (use `Model.ndim` instead).
-- `dims` and `coords` handling:
+
+  - `logpt`, `logpt_sum`, `logp_elemwiset` and `nojac` variations were removed. Use `Model.logpt(jacobian=True/False, sum=True/False)` instead.
+  - `dlogp_nojact` and `d2logp_nojact` were removed. Use `Model.dlogpt` and `d2logpt` with `jacobian=False` instead.
+  - `model.makefn` is now called `Model.compile_fn`, and `model.fn` was removed.
+  - Methods starting with `fast_*`, such as `Model.fast_logp`, were removed. Same applies to `PointFunc` classes
+  - `Model(model=...)` kwarg was removed
+  - `Model(theano_config=...)` kwarg was removed
+  - `Model.size` property was removed (use `Model.ndim` instead).
+
+  - `dims` and `coords` handling:
     - `Model.RV_dims` and `Model.coords` are now read-only properties. To modify the `coords` dictionary use `Model.add_coord`.
     - `dims` or coordinate values that are `None` will be auto-completed (see [#4625](https://github.com/pymc-devs/pymc/pull/4625)).
     - Coordinate values passed to `Model.add_coord` are always converted to tuples (see [#5061](https://github.com/pymc-devs/pymc/pull/5061)).
-- `Model.update_start_values(...)` was removed. Initial values can be set in the `Model.initial_values` dictionary directly.
-- Test values can no longer be set through `pm.Distribution(testval=...)` and must be assigned manually.
-- `Transform.forward` and `Transform.backward` signatures changed.
-- `pm.DensityDist` no longer accepts the `logp` as its first position argument. It is now an optional keyword argument. If you pass a callable as the first positional argument, a `TypeError` will be raised (see [5026](https://github.com/pymc-devs/pymc/pull/5026)).
-- `pm.DensityDist` now accepts distribution parameters as positional arguments. Passing them as a dictionary in the `observed` keyword argument is no longer supported and will raise an error (see [5026](https://github.com/pymc-devs/pymc/pull/5026)).
-- The signature of the `logp` and `random` functions that can be passed into a `pm.DensityDist` has been changed (see [5026](https://github.com/pymc-devs/pymc/pull/5026)).
-- Changes to the Gaussian process (`gp`) submodule:
-  - The `gp.prior(..., shape=...)` kwarg was renamed to `size`.
-  - Multiple methods including `gp.prior` now require explicit kwargs.
-- Changes to the BART implementation:
-  - A BART variable can be combined with other random variables. The `inv_link` argument has been removed (see [4914](https://github.com/pymc-devs/pymc3/pull/4914)).
-  - Moved BART to its own module (see [5058](https://github.com/pymc-devs/pymc3/pull/5058)).
-- Changes to the Gaussian Process (GP) submodule (see [5055](https://github.com/pymc-devs/pymc/pull/5055)):
-  - For all implementations, `gp.Latent`, `gp.Marginal` etc., `cov_func` and `mean_func` are required kwargs.
-  - In Windows test conda environment the `mkl` version is fixed to verison 2020.4, and `mkl-service` is fixed to `2.3.0`.  This was required for `gp.MarginalKron` to function properly.
-  - `gp.MvStudentT` uses rotated samples from `StudentT` directly now, instead of sampling from `pm.Chi2` and then from `pm.Normal`.
-  - The "jitter" parameter, or the diagonal noise term added to Gram matrices such that the Cholesky is numerically stable, is now exposed to the user instead of hard-coded.  See the function `gp.util.stabilize`.
-  - The `is_observed` arguement for `gp.Marginal*` implementations has been deprecated.
-  - In the gp.utils file, the `kmeans_inducing_points` function now passes through `kmeans_kwargs` to scipy's k-means function.
-  - The function `replace_with_values` function has been added to `gp.utils`.
-  - `MarginalSparse` has been renamed `MarginalApprox`.
-  - Removed `MixtureSameFamily`. `Mixture` is now capable of handling batched multivariate components (see [#5438](https://github.com/pymc-devs/pymc/pull/5438)).
-  - Removed `AR1`, `AR` of order 1 should be used instead. (see [5734](https://github.com/pymc-devs/pymc/pull/5734)).
-- ...
-
-### Expected breaks
-+ New API was already available in `v3`.
-+ Old API had deprecation warnings since at least `3.11.0` (2021-01).
-+ Old API stops working in `v4` (preferably with informative errors).
 
-All of the above apply to:
+  - `Transform.forward` and `Transform.backward` signatures changed.
 
-- `pm.sample(return_inferencedata=True)` is now the default (see [#4744](https://github.com/pymc-devs/pymc/pull/4744)).
-- ArviZ `plots` and `stats` *wrappers* were removed. The functions are now just available by their original names (see [#4549](https://github.com/pymc-devs/pymc/pull/4471) and `3.11.2` release notes).
-- `pm.sample_posterior_predictive(vars=...)` kwarg was removed in favor of `var_names` (see [#4343](https://github.com/pymc-devs/pymc/pull/4343)).
-- `ElemwiseCategorical` step method was removed (see [#4701](https://github.com/pymc-devs/pymc/pull/4701))
-- `LKJCholeskyCov` `compute_corr` keyword argument is now set to `True` by default (see[#5382](https://github.com/pymc-devs/pymc/pull/5382))
-
-### Ongoing deprecations
-- Old API still works in `v4` and has a deprecation warning.
-- Prefereably the new API should be available in `v3` already.
-
-This includes API changes we did not warn about since at least `3.11.0` (2021-01).
+  - Changes to the Gaussian Process (GP) submodule (see [5055](https://github.com/pymc-devs/pymc/pull/5055)):
+    - The `gp.prior(..., shape=...)` kwarg was renamed to `size`.
+    - Multiple methods including `gp.prior` now require explicit kwargs.
+    - For all implementations, `gp.Latent`, `gp.Marginal` etc., `cov_func` and `mean_func` are required kwargs.
+    - In Windows test conda environment the `mkl` version is fixed to verison 2020.4, and `mkl-service` is fixed to `2.3.0`.  This was required for `gp.MarginalKron` to function properly.
+    - `gp.MvStudentT` uses rotated samples from `StudentT` directly now, instead of sampling from `pm.Chi2` and then from `pm.Normal`.
+    - The "jitter" parameter, or the diagonal noise term added to Gram matrices such that the Cholesky is numerically stable, is now exposed to the user instead of hard-coded.  See the function `gp.util.stabilize`.
+    - The `is_observed` arguement for `gp.Marginal*` implementations has been deprecated.
+    - In the gp.utils file, the `kmeans_inducing_points` function now passes through `kmeans_kwargs` to scipy's k-means function.
+    - The function `replace_with_values` function has been added to `gp.utils`.
+    - `MarginalSparse` has been renamed `MarginalApprox`.
 
-- Setting initial values through `pm.Distribution(testval=...)` is now `pm.Distribution(initval=...)`.
-- Alternative `sd` keyword argument has been removed from all distributions. `sigma` should be used instead (see [#5583](https://github.com/pymc-devs/pymc/pull/5583)).
-
-
-### New features
-- The length of `dims` in the model is now tracked symbolically through `Model.dim_lengths` (see [#4625](https://github.com/pymc-devs/pymc/pull/4625)).
-- The `CAR` distribution has been added to allow for use of conditional autoregressions which often are used in spatial and network models.
-- The dimensionality of model variables can now be parametrized through either of `shape`, `dims` or `size` (see [#4696](https://github.com/pymc-devs/pymc/pull/4696)):
-  - With `shape` the length of dimensions must be given numerically or as scalar Aesara `Variables`. Numeric entries in `shape` restrict the model variable to the exact length and re-sizing is no longer possible.
-  - `dims` keeps model variables re-sizeable (for example through `pm.Data`) and leads to well defined coordinates in `InferenceData` objects.
-  - The `size` kwarg behaves like it does in Aesara/NumPy. For univariate RVs it is the same as `shape`, but for multivariate RVs it depends on how the RV implements broadcasting to dimensionality greater than `RVOp.ndim_supp`.
-  - An `Ellipsis` (`...`) in the last position of `shape` or `dims` can be used as short-hand notation for implied dimensions.
-- Added a `logcdf` implementation for the Kumaraswamy distribution (see [#4706](https://github.com/pymc-devs/pymc/pull/4706)).
-- The `OrderedMultinomial` distribution has been added for use on ordinal data which are _aggregated_ by trial, like multinomial observations, whereas `OrderedLogistic` only accepts ordinal data in a _disaggregated_ format, like categorical
-  observations (see [#4773](https://github.com/pymc-devs/pymc/pull/4773)).
-- The `Polya-Gamma` distribution has been added (see [#4531](https://github.com/pymc-devs/pymc/pull/4531)). To make use of this distribution, the [`polyagamma>=1.3.1`](https://pypi.org/project/polyagamma/) library must be installed and available in the user's environment.
-- A small change to the mass matrix tuning methods jitter+adapt_diag (the default) and adapt_diag improves performance early on during tuning for some models. [#5004](https://github.com/pymc-devs/pymc/pull/5004)
-- New experimental mass matrix tuning method jitter+adapt_diag_grad. [#5004](https://github.com/pymc-devs/pymc/pull/5004)
-- `pm.DensityDist` can now accept an optional `logcdf` keyword argument to pass in a function to compute the cummulative density function of the distribution (see [5026](https://github.com/pymc-devs/pymc/pull/5026)).
-- `pm.DensityDist` can now accept an optional `get_moment` keyword argument to pass in a function to compute the moment of the distribution (see [5026](https://github.com/pymc-devs/pymc/pull/5026)).
-- New features for BART:
-  - Added partial dependence plots and individual conditional expectation plots [5091](https://github.com/pymc-devs/pymc3/pull/5091).
-  - Modify how particle weights are computed. This improves accuracy of the modeled function (see [5177](https://github.com/pymc-devs/pymc3/pull/5177)).
-  - Improve sampling, increase default number of particles [5229](https://github.com/pymc-devs/pymc3/pull/5229).
-- The new `pm.find_constrained_prior` function can be used to find optimized prior parameters of a distribution under some
-  constraints (e.g lower and upper bound). See [#5231](https://github.com/pymc-devs/pymc/pull/5231).
-- New features for `pm.Data` containers:
-  - With `pm.Data(..., mutable=True/False)`, or by using `pm.MutableData` vs. `pm.ConstantData` one can now create `TensorConstant` data variables. They can be more performant and compatible in situations where a variable doesn't need to be changed via `pm.set_data()`. See [#5295](https://github.com/pymc-devs/pymc/pull/5295).
-  - New named dimensions can be introduced to the model via `pm.Data(..., dims=...)`. For mutable data variables (see above) the lengths of these dimensions are symbolic, so they can be re-sized via `pm.set_data()`.
-  - `pm.Data` now passes additional kwargs to `aesara.shared`/`at.as_tensor`. [#5098](https://github.com/pymc-devs/pymc/pull/5098).
-- Univariate censored distributions are now available via `pm.Censored`. [#5169](https://github.com/pymc-devs/pymc/pull/5169)
-- Nested models now inherit the parent model's coordinates. [#5344](https://github.com/pymc-devs/pymc/pull/5344)
-- `softmax` and `log_softmax` functions added to `math` module (see [#5279](https://github.com/pymc-devs/pymc/pull/5279)).
-- Adding support for blackjax's NUTS sampler `pymc.sampling_jax` (see [#5477](ihttps://github.com/pymc-devs/pymc/pull/5477))
-- `pymc.sampling_jax` samplers support `log_likelihood`, `observed_data`, and `sample_stats` in returned InferenceData object (see [#5189](https://github.com/pymc-devs/pymc/pull/5189))
-- Adding support for `pm.Deterministic` in `pymc.sampling_jax` (see [#5182](https://github.com/pymc-devs/pymc/pull/5182))
-- Added an alternative parametrization, `logit_p` to `pm.Binomial` and `pm.Categorical` distributions (see [5637](https://github.com/pymc-devs/pymc/pull/5637)).
-- Added the low level `compile_forward_sampling_function` method to compile the aesara function responsible for generating forward samples (see [#5759](https://github.com/pymc-devs/pymc/pull/5759)).
-- ...
+  - Removed `MixtureSameFamily`. `Mixture` is now capable of handling batched multivariate components (see [#5438](https://github.com/pymc-devs/pymc/pull/5438)).
 
 ### Documentation
 - Switched to the [pydata-sphinx-theme](https://pydata-sphinx-theme.readthedocs.io/en/latest/)
@@ -157,8 +151,7 @@ This includes API changes we did not warn about since at least `3.11.0` (2021-01
 - Updated API docs to document objects at the path users should use to import them
 
 ### Maintenance
-- ⚠ Fixed old-time bug in Slice sampler that resulted in biased samples (see [#5816](https://github.com/pymc-devs/pymc/pull/5816)).
-- ⚠ PyMC now requires Scipy version `>= 1.4.1` (see [4857](https://github.com/pymc-devs/pymc/pull/4857)).
+- ⚠️ Fixed old-time bug in Slice sampler that resulted in biased samples (see [#5816](https://github.com/pymc-devs/pymc/pull/5816)).
 - Removed float128 dtype support (see [#4514](https://github.com/pymc-devs/pymc/pull/4514)).
 - Logp method of `Uniform` and `DiscreteUniform` no longer depends on `pymc.distributions.dist_math.bound` for proper evaluation (see [#4541](https://github.com/pymc-devs/pymc/pull/4541)).
 - We now include `cloudpickle` as a required dependency, and no longer depend on `dill` (see [#4858](https://github.com/pymc-devs/pymc/pull/4858)).
@@ -167,8 +160,6 @@ This includes API changes we did not warn about since at least `3.11.0` (2021-01
 - Changed name of `Lognormal` distribution to `LogNormal` to harmonize CamelCase usage for distribution names.
 - Attempt to iterate over MultiTrace will raise NotImplementedError.
 - Removed silent normalisation of `p` parameters in Categorical and Multinomial distributions (see [#5370](https://github.com/pymc-devs/pymc/pull/5370)).
-- ...
-
 
 ## PyMC 3.11.2 (14 March 2021)
 
diff --git a/docs/source/learn/books.md b/docs/source/learn/books.md
index 8a0e03599..8d7f92fb8 100644
--- a/docs/source/learn/books.md
+++ b/docs/source/learn/books.md
@@ -4,6 +4,17 @@
 ::::{grid} 1 2 2 3
 :gutter: 3
 
+:::{grid-item-card} Bayesian Modeling and Computation in Python
+:img-top: https://bayesiancomputationbook.com/_images/Cover.jpg
+
+By Osvaldo Martin, Ravin Kumar and Junpeng Lao
+
+Hands on approach with PyMC and ArviZ focusing on the practice of applied statistics.
+
+[Website + code](https://bayesiancomputationbook.com/welcome.html)
+
+:::
+
 :::{grid-item-card} Bayesian Methods for Hackers
 :img-top: https://camo.githubusercontent.com/4a0aca82ca82efab71747d00db30f3a68de98e82/687474703a2f2f692e696d6775722e636f6d2f36444b596250622e706e673f31
 By Cameron Davidson-Pilon
@@ -56,17 +67,6 @@ A Bayesian Course with Examples in R and Stan.
 
 :::
 
-:::{grid-item-card} Bayesian Modeling and Computation in Python
-:img-top: https://bayesiancomputationbook.com/_images/Cover.jpg
-
-By Osvaldo Martin, Ravin Kumar and Junpeng Lao
-
-Hands on approach with PyMC and ArviZ focusing on the practice of applied statistics.
-
-[Website + code](https://bayesiancomputationbook.com/welcome.html)
-
-:::
-
 :::{grid-item-card} Bayesian Cognitive Modeling: A Practical Course
 :img-top: https://images-na.ssl-images-amazon.com/images/I/51K33XI2I8L._SX330_BO1,204,203,200_.jpg
 
diff --git a/pymc/distributions/distribution.py b/pymc/distributions/distribution.py
index 22d766257..dd5968df7 100644
--- a/pymc/distributions/distribution.py
+++ b/pymc/distributions/distribution.py
@@ -47,7 +47,7 @@ from pymc.distributions.shape_utils import (
     resize_from_dims,
     resize_from_observed,
 )
-from pymc.printing import str_for_dist
+from pymc.printing import str_for_dist, str_for_symbolic_dist
 from pymc.util import UNSET
 from pymc.vartypes import string_types
 
@@ -483,15 +483,11 @@ class SymbolicDistribution:
             transform=transform,
             initval=initval,
         )
-
-        # TODO: Refactor this
         # add in pretty-printing support
-        rv_out.str_repr = lambda *args, **kwargs: name
-        rv_out._repr_latex_ = f"\\text{name}"
-        # rv_out.str_repr = types.MethodType(str_for_dist, rv_out)
-        # rv_out._repr_latex_ = types.MethodType(
-        #     functools.partial(str_for_dist, formatting="latex"), rv_out
-        # )
+        rv_out.str_repr = types.MethodType(str_for_symbolic_dist, rv_out)
+        rv_out._repr_latex_ = types.MethodType(
+            functools.partial(str_for_symbolic_dist, formatting="latex"), rv_out
+        )
 
         return rv_out
 
diff --git a/pymc/model.py b/pymc/model.py
index 185e30644..4f089e573 100644
--- a/pymc/model.py
+++ b/pymc/model.py
@@ -60,7 +60,7 @@ from pymc.data import GenTensorVariable, Minibatch
 from pymc.distributions import joint_logpt
 from pymc.distributions.logprob import _get_scaling
 from pymc.distributions.transforms import _default_transform
-from pymc.exceptions import ImputationWarning, SamplingError, ShapeError
+from pymc.exceptions import ImputationWarning, SamplingError, ShapeError, ShapeWarning
 from pymc.initial_point import make_initial_point_fn
 from pymc.math import flatten_list
 from pymc.util import (
@@ -1067,12 +1067,6 @@ class Model(WithMemoization, metaclass=ContextMeta):
             raise ValueError(
                 f"Either `values` or `length` must be specified for the '{name}' dimension."
             )
-        if isinstance(length, int):
-            length = at.constant(length)
-        elif length is not None and not isinstance(length, Variable):
-            raise ValueError(
-                f"The `length` passed for the '{name}' coord must be an Aesara Variable or None."
-            )
         if values is not None:
             # Conversion to a tuple ensures that the coordinate values are immutable.
             # Also unlike numpy arrays the's tuple.index(...) which is handy to work with.
@@ -1080,12 +1074,19 @@ class Model(WithMemoization, metaclass=ContextMeta):
         if name in self.coords:
             if not np.array_equal(values, self.coords[name]):
                 raise ValueError(f"Duplicate and incompatible coordinate: {name}.")
-        else:
+        if length is not None and not isinstance(length, (int, Variable)):
+            raise ValueError(
+                f"The `length` passed for the '{name}' coord must be an int, Aesara Variable or None."
+            )
+        if length is None:
+            length = len(values)
+        if not isinstance(length, Variable):
             if mutable:
-                self._dim_lengths[name] = length or aesara.shared(len(values))
+                length = aesara.shared(length)
             else:
-                self._dim_lengths[name] = length or aesara.tensor.constant(len(values))
-            self._coords[name] = values
+                length = aesara.tensor.constant(length)
+        self._dim_lengths[name] = length
+        self._coords[name] = values
 
     def add_coords(
         self,
@@ -1101,6 +1102,36 @@ class Model(WithMemoization, metaclass=ContextMeta):
         for name, values in coords.items():
             self.add_coord(name, values, length=lengths.get(name, None))
 
+    def set_dim(self, name: str, new_length: int, coord_values: Optional[Sequence] = None):
+        """Update a mutable dimension.
+
+        Parameters
+        ----------
+        name
+            Name of the dimension.
+        new_length
+            New length of the dimension.
+        coord_values
+            Optional sequence of coordinate values.
+        """
+        if not isinstance(self.dim_lengths[name], ScalarSharedVariable):
+            raise ValueError(f"The dimension '{name}' is immutable.")
+        if coord_values is None and self.coords.get(name, None) is not None:
+            raise ValueError(
+                f"'{name}' has coord values. Pass `set_dim(..., coord_values=...)` to update them."
+            )
+        if coord_values is not None:
+            len_cvals = len(coord_values)
+            if len_cvals != new_length:
+                raise ShapeError(
+                    f"Length of new coordinate values does not match the new dimension length.",
+                    actual=len_cvals,
+                    expected=new_length,
+                )
+            self._coords[name] = tuple(coord_values)
+        self.dim_lengths[name].set_value(new_length)
+        return
+
     def set_data(
         self,
         name: str,
@@ -1163,16 +1194,40 @@ class Model(WithMemoization, metaclass=ContextMeta):
                             f"{new_length}, so new coord values for the {dname} dimension are required."
                         )
                 if isinstance(length_tensor, TensorConstant):
+                    # The dimension was fixed in length.
+                    # Resizing a data variable in this dimension would
+                    # definitely lead to shape problems.
                     raise ShapeError(
                         f"Resizing dimension '{dname}' is impossible, because "
-                        "a 'TensorConstant' stores its length. To be able "
+                        "a `TensorConstant` stores its length. To be able "
                         "to change the dimension length, pass `mutable=True` when "
                         "registering the dimension via `model.add_coord`, "
                         "or define it via a `pm.MutableData` variable."
                     )
-                else:
+                elif length_tensor.owner is not None:
+                    # The dimension was created from a model variable.
+                    # Get a handle on the tensor from which this dimension length was
+                    # obtained by doing subindexing on the shape as in `.shape[i]`.
+                    # Needed to check if it was another shared variable.
+                    # TODO: Consider checking the graph is what we are assuming it is
+                    # isinstance(length_tensor.owner.op, Subtensor)
+                    # isinstance(length_tensor.owner.inputs[0].owner.op, Shape)
                     length_belongs_to = length_tensor.owner.inputs[0].owner.inputs[0]
-                    if not isinstance(length_belongs_to, SharedVariable):
+                    if length_belongs_to is shared_object:
+                        # No surprise it's changing.
+                        pass
+                    elif isinstance(length_belongs_to, SharedVariable):
+                        # The dimension is mutable through a SharedVariable other than the one being modified.
+                        # But the other variable was not yet re-sized! Warn the user to do that!
+                        warnings.warn(
+                            f"You are resizing a variable with dimension '{dname}' which was initialized "
+                            f"as a mutable dimension by another variable ('{length_belongs_to}')."
+                            " Remember to update that variable with the correct shape to avoid shape issues.",
+                            ShapeWarning,
+                            stacklevel=2,
+                        )
+                    else:
+                        # The dimension is immutable.
                         raise ShapeError(
                             f"Resizing dimension '{dname}' with values of length {new_length} would lead to incompatibilities, "
                             f"because the dimension was initialized from '{length_belongs_to}' which is not a shared variable. "
@@ -1182,8 +1237,9 @@ class Model(WithMemoization, metaclass=ContextMeta):
                             expected=old_length,
                         )
                 if isinstance(length_tensor, ScalarSharedVariable):
-                    # Updating the shared variable resizes dependent nodes that use this dimension for their `size`.
-                    length_tensor.set_value(new_length)
+                    # The dimension is mutable, but was defined without being linked
+                    # to a shared variable. This is allowed, but a little less robust.
+                    self.set_dim(dname, new_length, coord_values=new_coords)
 
             if new_coords is not None:
                 # Update the registered coord values (also if they were None)
diff --git a/pymc/printing.py b/pymc/printing.py
index 828d07cee..87c43924e 100644
--- a/pymc/printing.py
+++ b/pymc/printing.py
@@ -55,6 +55,18 @@ def str_for_dist(rv: TensorVariable, formatting: str = "plain", include_params:
             return rf"{print_name} ~ {dist_name}"
 
 
+def str_for_symbolic_dist(
+    rv: TensorVariable, formatting: str = "plain", include_params: bool = True
+) -> str:
+    """Make a human-readable string representation of a SymbolicDistribution in a model,
+    either LaTeX or plain, optionally with distribution parameter values included."""
+
+    if "latex" in formatting:
+        return rf"$\text{{{rv.name}}} \sim \text{{{rv.owner.op}}}$"
+    else:
+        return rf"{rv.name} ~ {rv.owner.op}"
+
+
 def str_for_model(model: Model, formatting: str = "plain", include_params: bool = True) -> str:
     """Make a human-readable string representation of Model, listing all random variables
     and their distributions, optionally including parameter values."""


Test patch:
diff --git a/pymc/tests/test_distributions.py b/pymc/tests/test_distributions.py
index ba1611404..06ba3428d 100644
--- a/pymc/tests/test_distributions.py
+++ b/pymc/tests/test_distributions.py
@@ -127,7 +127,7 @@ from pymc.distributions import (
 )
 from pymc.distributions.shape_utils import to_tuple
 from pymc.math import kronecker
-from pymc.model import Deterministic, Model, Point, Potential
+from pymc.model import Model, Point
 from pymc.tests.helpers import select_by_precision
 from pymc.vartypes import continuous_types, discrete_types
 
@@ -2886,139 +2886,6 @@ class TestBoundedContinuous:
         assert upper_interval is None
 
 
-class TestStrAndLatexRepr:
-    def setup_class(self):
-        # True parameter values
-        alpha, sigma = 1, 1
-        beta = [1, 2.5]
-
-        # Size of dataset
-        size = 100
-
-        # Predictor variable
-        X = np.random.normal(size=(size, 2)).dot(np.array([[1, 0], [0, 0.2]]))
-
-        # Simulate outcome variable
-        Y = alpha + X.dot(beta) + np.random.randn(size) * sigma
-        with Model() as self.model:
-            # TODO: some variables commented out here as they're not working properly
-            # in v4 yet (9-jul-2021), so doesn't make sense to test str/latex for them
-
-            # Priors for unknown model parameters
-            alpha = Normal("alpha", mu=0, sigma=10)
-            b = Normal("beta", mu=0, sigma=10, size=(2,), observed=beta)
-            sigma = HalfNormal("sigma", sigma=1)
-
-            # Test Cholesky parameterization
-            Z = MvNormal("Z", mu=np.zeros(2), chol=np.eye(2), size=(2,))
-
-            # NegativeBinomial representations to test issue 4186
-            # nb1 = pm.NegativeBinomial(
-            #     "nb_with_mu_alpha", mu=pm.Normal("nbmu"), alpha=pm.Gamma("nbalpha", mu=6, sigma=1)
-            # )
-            nb2 = pm.NegativeBinomial("nb_with_p_n", p=pm.Uniform("nbp"), n=10)
-
-            # Expected value of outcome
-            mu = Deterministic("mu", floatX(alpha + at.dot(X, b)))
-
-            # add a bounded variable as well
-            # bound_var = Bound(Normal, lower=1.0)("bound_var", mu=0, sigma=10)
-
-            # KroneckerNormal
-            n, m = 3, 4
-            covs = [np.eye(n), np.eye(m)]
-            kron_normal = KroneckerNormal("kron_normal", mu=np.zeros(n * m), covs=covs, size=n * m)
-
-            # MatrixNormal
-            # matrix_normal = MatrixNormal(
-            #     "mat_normal",
-            #     mu=np.random.normal(size=n),
-            #     rowcov=np.eye(n),
-            #     colchol=np.linalg.cholesky(np.eye(n)),
-            #     size=(n, n),
-            # )
-
-            # DirichletMultinomial
-            dm = DirichletMultinomial("dm", n=5, a=[1, 1, 1], size=(2, 3))
-
-            # Likelihood (sampling distribution) of observations
-            Y_obs = Normal("Y_obs", mu=mu, sigma=sigma, observed=Y)
-
-            # add a potential as well
-            pot = Potential("pot", mu**2)
-
-        self.distributions = [alpha, sigma, mu, b, Z, nb2, Y_obs, pot]
-        self.deterministics_or_potentials = [mu, pot]
-        # tuples of (formatting, include_params
-        self.formats = [("plain", True), ("plain", False), ("latex", True), ("latex", False)]
-        self.expected = {
-            ("plain", True): [
-                r"alpha ~ N(0, 10)",
-                r"sigma ~ N**+(0, 1)",
-                r"mu ~ Deterministic(f(beta, alpha))",
-                r"beta ~ N(0, 10)",
-                r"Z ~ N(f(), f())",
-                r"nb_with_p_n ~ NB(10, nbp)",
-                r"Y_obs ~ N(mu, sigma)",
-                r"pot ~ Potential(f(beta, alpha))",
-            ],
-            ("plain", False): [
-                r"alpha ~ N",
-                r"sigma ~ N**+",
-                r"mu ~ Deterministic",
-                r"beta ~ N",
-                r"Z ~ N",
-                r"nb_with_p_n ~ NB",
-                r"Y_obs ~ N",
-                r"pot ~ Potential",
-            ],
-            ("latex", True): [
-                r"$\text{alpha} \sim \operatorname{N}(0,~10)$",
-                r"$\text{sigma} \sim \operatorname{N^{+}}(0,~1)$",
-                r"$\text{mu} \sim \operatorname{Deterministic}(f(\text{beta},~\text{alpha}))$",
-                r"$\text{beta} \sim \operatorname{N}(0,~10)$",
-                r"$\text{Z} \sim \operatorname{N}(f(),~f())$",
-                r"$\text{nb_with_p_n} \sim \operatorname{NB}(10,~\text{nbp})$",
-                r"$\text{Y_obs} \sim \operatorname{N}(\text{mu},~\text{sigma})$",
-                r"$\text{pot} \sim \operatorname{Potential}(f(\text{beta},~\text{alpha}))$",
-            ],
-            ("latex", False): [
-                r"$\text{alpha} \sim \operatorname{N}$",
-                r"$\text{sigma} \sim \operatorname{N^{+}}$",
-                r"$\text{mu} \sim \operatorname{Deterministic}$",
-                r"$\text{beta} \sim \operatorname{N}$",
-                r"$\text{Z} \sim \operatorname{N}$",
-                r"$\text{nb_with_p_n} \sim \operatorname{NB}$",
-                r"$\text{Y_obs} \sim \operatorname{N}$",
-                r"$\text{pot} \sim \operatorname{Potential}$",
-            ],
-        }
-
-    def test__repr_latex_(self):
-        for distribution, tex in zip(self.distributions, self.expected[("latex", True)]):
-            assert distribution._repr_latex_() == tex
-
-        model_tex = self.model._repr_latex_()
-
-        # make sure each variable is in the model
-        for tex in self.expected[("latex", True)]:
-            for segment in tex.strip("$").split(r"\sim"):
-                assert segment in model_tex
-
-    def test_str_repr(self):
-        for str_format in self.formats:
-            for dist, text in zip(self.distributions, self.expected[str_format]):
-                assert dist.str_repr(*str_format) == text
-
-            model_text = self.model.str_repr(*str_format)
-            for text in self.expected[str_format]:
-                if str_format[0] == "latex":
-                    for segment in text.strip("$").split(r"\sim"):
-                        assert segment in model_text
-                else:
-                    assert text in model_text
-
-
 def test_discrete_trafo():
     with Model():
         with pytest.raises(ValueError) as err:
diff --git a/pymc/tests/test_model.py b/pymc/tests/test_model.py
index a4c218190..2f8683e9e 100644
--- a/pymc/tests/test_model.py
+++ b/pymc/tests/test_model.py
@@ -12,6 +12,7 @@
 #   See the License for the specific language governing permissions and
 #   limitations under the License.
 import unittest
+import warnings
 
 from functools import reduce
 
@@ -37,7 +38,7 @@ import pymc as pm
 from pymc import Deterministic, Potential
 from pymc.blocking import DictToArrayBijection, RaveledVars
 from pymc.distributions import Normal, transforms
-from pymc.exceptions import ShapeError
+from pymc.exceptions import ShapeError, ShapeWarning
 from pymc.model import Point, ValueGradFunction
 from pymc.tests.helpers import SeededTest
 
@@ -704,25 +705,45 @@ def test_nested_model_coords():
     assert set(m2.RV_dims) < set(m1.RV_dims)
 
 
-def test_shapeerror_from_resize_immutable_dims():
+def test_shapeerror_from_set_data_dimensionality():
+    with pm.Model() as pmodel:
+        pm.MutableData("m", np.ones((3,)), dims="one")
+        with pytest.raises(ValueError, match="must have 1 dimensions"):
+            pmodel.set_data("m", np.ones((3, 4)))
+
+
+def test_shapeerror_from_resize_immutable_dim_from_RV():
     """
     Trying to resize an immutable dimension should raise a ShapeError.
     Even if the variable being updated is a SharedVariable and has other
     dimensions that are mutable.
     """
     with pm.Model() as pmodel:
-        a = pm.Normal("a", mu=[1, 2, 3], dims="fixed")
+        pm.Normal("a", mu=[1, 2, 3], dims="fixed")
+        assert isinstance(pmodel.dim_lengths["fixed"], TensorVariable)
 
-        m = pm.MutableData("m", [[1, 2, 3]], dims=("one", "fixed"))
+        pm.MutableData("m", [[1, 2, 3]], dims=("one", "fixed"))
 
         # This is fine because the "fixed" dim is not resized
-        pm.set_data({"m": [[1, 2, 3], [3, 4, 5]]})
+        pmodel.set_data("m", [[1, 2, 3], [3, 4, 5]])
 
     with pytest.raises(ShapeError, match="was initialized from 'a'"):
-        # Can't work because the "fixed" dimension is linked to a constant shape:
+        # Can't work because the "fixed" dimension is linked to a
+        # TensorVariable with constant shape.
         # Note that the new data tries to change both dimensions
-        with pmodel:
-            pm.set_data({"m": [[1, 2], [3, 4]]})
+        pmodel.set_data("m", [[1, 2], [3, 4]])
+
+
+def test_shapeerror_from_resize_immutable_dim_from_coords():
+    with pm.Model(coords={"immutable": [1, 2]}) as pmodel:
+        assert isinstance(pmodel.dim_lengths["immutable"], TensorConstant)
+        pm.MutableData("m", [1, 2], dims="immutable")
+        # Data can be changed
+        pmodel.set_data("m", [3, 4])
+
+    with pytest.raises(ShapeError, match="`TensorConstant` stores its length"):
+        # But the length is linked to a TensorConstant
+        pmodel.set_data("m", [1, 2, 3], coords=dict(immutable=[1, 2, 3]))
 
 
 def test_valueerror_from_resize_without_coords_update():
@@ -765,6 +786,100 @@ def test_add_coord_mutable_kwarg():
         assert isinstance(m._dim_lengths["mutable2"], TensorVariable)
 
 
+def test_set_dim():
+    """Test the concious re-sizing of dims created through add_coord()."""
+    with pm.Model() as pmodel:
+        pmodel.add_coord("fdim", mutable=False, length=1)
+        pmodel.add_coord("mdim", mutable=True, length=2)
+        a = pm.Normal("a", dims="mdim")
+    assert a.eval().shape == (2,)
+
+    with pytest.raises(ValueError, match="is immutable"):
+        pmodel.set_dim("fdim", 3)
+
+    pmodel.set_dim("mdim", 3)
+    assert a.eval().shape == (3,)
+
+
+def test_set_dim_with_coords():
+    """Test the concious re-sizing of dims created through add_coord() with coord value."""
+    with pm.Model() as pmodel:
+        pmodel.add_coord("mdim", mutable=True, length=2, values=["A", "B"])
+        a = pm.Normal("a", dims="mdim")
+    assert len(pmodel.coords["mdim"]) == 2
+
+    with pytest.raises(ValueError, match="has coord values"):
+        pmodel.set_dim("mdim", new_length=3)
+
+    with pytest.raises(ShapeError, match="does not match"):
+        pmodel.set_dim("mdim", new_length=3, coord_values=["A", "B"])
+
+    pmodel.set_dim("mdim", 3, ["A", "B", "C"])
+    assert a.eval().shape == (3,)
+    assert pmodel.coords["mdim"] == ("A", "B", "C")
+
+
+def test_set_data_indirect_resize():
+    with pm.Model() as pmodel:
+        pmodel.add_coord("mdim", mutable=True, length=2)
+        pm.MutableData("mdata", [1, 2], dims="mdim")
+
+    # First resize the dimension.
+    pmodel.dim_lengths["mdim"].set_value(3)
+    # Then change the data.
+    with warnings.catch_warnings():
+        warnings.simplefilter("error")
+        pmodel.set_data("mdata", [1, 2, 3])
+
+    # Now the other way around.
+    with warnings.catch_warnings():
+        warnings.simplefilter("error")
+        pmodel.set_data("mdata", [1, 2, 3, 4])
+
+
+def test_set_data_warns_on_resize_of_dims_defined_by_other_mutabledata():
+    with pm.Model() as pmodel:
+        pm.MutableData("m1", [1, 2], dims="mutable")
+        pm.MutableData("m2", [3, 4], dims="mutable")
+
+        # Resizing the non-defining variable first gives a warning
+        with pytest.warns(ShapeWarning, match="by another variable"):
+            pmodel.set_data("m2", [4, 5, 6])
+            pmodel.set_data("m1", [1, 2, 3])
+
+        # Resizing the definint variable first is silent
+        with warnings.catch_warnings():
+            warnings.simplefilter("error")
+            pmodel.set_data("m1", [1, 2])
+            pmodel.set_data("m2", [3, 4])
+
+
+def test_set_data_indirect_resize_with_coords():
+    with pm.Model() as pmodel:
+        pmodel.add_coord("mdim", ["A", "B"], mutable=True, length=2)
+        pm.MutableData("mdata", [1, 2], dims="mdim")
+
+    assert pmodel.coords["mdim"] == ("A", "B")
+
+    # First resize the dimension.
+    pmodel.set_dim("mdim", 3, ["A", "B", "C"])
+    assert pmodel.coords["mdim"] == ("A", "B", "C")
+    # Then change the data.
+    with warnings.catch_warnings():
+        warnings.simplefilter("error")
+        pmodel.set_data("mdata", [1, 2, 3])
+
+    # Now the other way around.
+    with warnings.catch_warnings():
+        warnings.simplefilter("error")
+        pmodel.set_data("mdata", [1, 2, 3, 4], coords=dict(mdim=["A", "B", "C", "D"]))
+    assert pmodel.coords["mdim"] == ("A", "B", "C", "D")
+
+    # This time with incorrectly sized coord values
+    with pytest.raises(ShapeError, match="new coordinate values"):
+        pmodel.set_data("mdata", [1, 2], coords=dict(mdim=[1, 2, 3]))
+
+
 @pytest.mark.parametrize("jacobian", [True, False])
 def test_model_logp(jacobian):
     with pm.Model() as m:
diff --git a/pymc/tests/test_printing.py b/pymc/tests/test_printing.py
new file mode 100644
index 000000000..2d082caa3
--- /dev/null
+++ b/pymc/tests/test_printing.py
@@ -0,0 +1,156 @@
+import numpy as np
+
+from pymc.aesaraf import floatX
+from pymc.distributions import (
+    DirichletMultinomial,
+    HalfNormal,
+    KroneckerNormal,
+    MvNormal,
+    NegativeBinomial,
+    Normal,
+    Uniform,
+    ZeroInflatedPoisson,
+)
+from pymc.math import dot
+from pymc.model import Deterministic, Model, Potential
+
+
+# TODO: This test is a bit too monolithic
+class TestStrAndLatexRepr:
+    def setup_class(self):
+        # True parameter values
+        alpha, sigma = 1, 1
+        beta = [1, 2.5]
+
+        # Size of dataset
+        size = 100
+
+        # Predictor variable
+        X = np.random.normal(size=(size, 2)).dot(np.array([[1, 0], [0, 0.2]]))
+
+        # Simulate outcome variable
+        Y = alpha + X.dot(beta) + np.random.randn(size) * sigma
+        with Model() as self.model:
+            # TODO: some variables commented out here as they're not working properly
+            # in v4 yet (9-jul-2021), so doesn't make sense to test str/latex for them
+
+            # Priors for unknown model parameters
+            alpha = Normal("alpha", mu=0, sigma=10)
+            b = Normal("beta", mu=0, sigma=10, size=(2,), observed=beta)
+            sigma = HalfNormal("sigma", sigma=1)
+
+            # Test Cholesky parameterization
+            Z = MvNormal("Z", mu=np.zeros(2), chol=np.eye(2), size=(2,))
+
+            # NegativeBinomial representations to test issue 4186
+            # nb1 = pm.NegativeBinomial(
+            #     "nb_with_mu_alpha", mu=pm.Normal("nbmu"), alpha=pm.Gamma("nbalpha", mu=6, sigma=1)
+            # )
+            nb2 = NegativeBinomial("nb_with_p_n", p=Uniform("nbp"), n=10)
+
+            # Symbolic distribution
+            zip = ZeroInflatedPoisson("zip", 0.5, 5)
+
+            # Expected value of outcome
+            mu = Deterministic("mu", floatX(alpha + dot(X, b)))
+
+            # add a bounded variable as well
+            # bound_var = Bound(Normal, lower=1.0)("bound_var", mu=0, sigma=10)
+
+            # KroneckerNormal
+            n, m = 3, 4
+            covs = [np.eye(n), np.eye(m)]
+            kron_normal = KroneckerNormal("kron_normal", mu=np.zeros(n * m), covs=covs, size=n * m)
+
+            # MatrixNormal
+            # matrix_normal = MatrixNormal(
+            #     "mat_normal",
+            #     mu=np.random.normal(size=n),
+            #     rowcov=np.eye(n),
+            #     colchol=np.linalg.cholesky(np.eye(n)),
+            #     size=(n, n),
+            # )
+
+            # DirichletMultinomial
+            dm = DirichletMultinomial("dm", n=5, a=[1, 1, 1], size=(2, 3))
+
+            # Likelihood (sampling distribution) of observations
+            Y_obs = Normal("Y_obs", mu=mu, sigma=sigma, observed=Y)
+
+            # add a potential as well
+            pot = Potential("pot", mu**2)
+
+        self.distributions = [alpha, sigma, mu, b, Z, nb2, zip, Y_obs, pot]
+        self.deterministics_or_potentials = [mu, pot]
+        # tuples of (formatting, include_params
+        self.formats = [("plain", True), ("plain", False), ("latex", True), ("latex", False)]
+        self.expected = {
+            ("plain", True): [
+                r"alpha ~ N(0, 10)",
+                r"sigma ~ N**+(0, 1)",
+                r"mu ~ Deterministic(f(beta, alpha))",
+                r"beta ~ N(0, 10)",
+                r"Z ~ N(f(), f())",
+                r"nb_with_p_n ~ NB(10, nbp)",
+                r"zip ~ MarginalMixtureRV{inline=False}",
+                r"Y_obs ~ N(mu, sigma)",
+                r"pot ~ Potential(f(beta, alpha))",
+            ],
+            ("plain", False): [
+                r"alpha ~ N",
+                r"sigma ~ N**+",
+                r"mu ~ Deterministic",
+                r"beta ~ N",
+                r"Z ~ N",
+                r"nb_with_p_n ~ NB",
+                r"zip ~ MarginalMixtureRV{inline=False}",
+                r"Y_obs ~ N",
+                r"pot ~ Potential",
+            ],
+            ("latex", True): [
+                r"$\text{alpha} \sim \operatorname{N}(0,~10)$",
+                r"$\text{sigma} \sim \operatorname{N^{+}}(0,~1)$",
+                r"$\text{mu} \sim \operatorname{Deterministic}(f(\text{beta},~\text{alpha}))$",
+                r"$\text{beta} \sim \operatorname{N}(0,~10)$",
+                r"$\text{Z} \sim \operatorname{N}(f(),~f())$",
+                r"$\text{nb_with_p_n} \sim \operatorname{NB}(10,~\text{nbp})$",
+                r"$\text{zip} \sim \text{MarginalMixtureRV{inline=False}}$",
+                r"$\text{Y_obs} \sim \operatorname{N}(\text{mu},~\text{sigma})$",
+                r"$\text{pot} \sim \operatorname{Potential}(f(\text{beta},~\text{alpha}))$",
+            ],
+            ("latex", False): [
+                r"$\text{alpha} \sim \operatorname{N}$",
+                r"$\text{sigma} \sim \operatorname{N^{+}}$",
+                r"$\text{mu} \sim \operatorname{Deterministic}$",
+                r"$\text{beta} \sim \operatorname{N}$",
+                r"$\text{Z} \sim \operatorname{N}$",
+                r"$\text{nb_with_p_n} \sim \operatorname{NB}$",
+                r"$\text{zip} \sim \text{MarginalMixtureRV{inline=False}}$",
+                r"$\text{Y_obs} \sim \operatorname{N}$",
+                r"$\text{pot} \sim \operatorname{Potential}$",
+            ],
+        }
+
+    def test__repr_latex_(self):
+        for distribution, tex in zip(self.distributions, self.expected[("latex", True)]):
+            assert distribution._repr_latex_() == tex
+
+        model_tex = self.model._repr_latex_()
+
+        # make sure each variable is in the model
+        for tex in self.expected[("latex", True)]:
+            for segment in tex.strip("$").split(r"\sim"):
+                assert segment in model_tex
+
+    def test_str_repr(self):
+        for str_format in self.formats:
+            for dist, text in zip(self.distributions, self.expected[str_format]):
+                assert dist.str_repr(*str_format) == text
+
+            model_text = self.model.str_repr(*str_format)
+            for text in self.expected[str_format]:
+                if str_format[0] == "latex":
+                    for segment in text.strip("$").split(r"\sim"):
+                        assert segment in model_text
+                else:
+                    assert text in model_text