Repository: pymc-devs/pymc

Issue:
set_data cannot deal with nodes without owners
## Description of your problem

`set_data` seems to have a problem when the new coordinates are nodes without an owner.

**Please provide a minimal, self-contained, and reproducible example.**
```python
import numpy as np
import pymc as pm


with pm.Model() as m:
    m.add_coord("a", np.arange(10), mutable=True)
    pm.MutableData("x", np.zeros(10), dims="a")

m.set_data("x", np.ones(5), coords={"a": np.arange(10, 15)})
```

**Please provide the full traceback.**

<details><summary>Complete error traceback</summary>

```python
AttributeError                            Traceback (most recent call last)
<ipython-input-1-dc6231ed80dd> in <module>
      7     pm.MutableData("x", np.zeros(10), dims="a")
      8 
----> 9 m.set_data("x", np.ones(5), coords={"a": np.arange(10, 15)})

~/repos/pymc3/pymc/model.py in set_data(self, name, values, coords)
   1172                     )
   1173                 else:
-> 1174                     length_belongs_to = length_tensor.owner.inputs[0].owner.inputs[0]
   1175                     if not isinstance(length_belongs_to, SharedVariable):
   1176                         raise ShapeError(

AttributeError: 'NoneType' object has no attribute 'inputs'
```

</details>

**Please provide any additional information below.**

In the traceback, `length_tensor` is `m.dim_lengths["a"]`. The type of this variable is `aesara.tensor.sharedvar.ScalarSharedVariable` and it has no owners. I don't know when the owner logic was introduced, but it doesn't cover the simple case where the new coordinate is an array. Maybe this was intended? In any case, this should be fixed before releasing v4

### Original code

Originally, I was running on an older version with respect to main. The code that I used was this

```python
import numpy as np
import pymc as pm


coords = {"a": np.arange(10)}
with pm.Model(coords=coords) as m:
    pm.MutableData("x", np.zeros(10), dims="a")

m.set_data("x", np.ones(5), coords={"a": np.arange(10, 15)})
```

<details><summary>Which now leads to the following traceback instead</summary>

```python
ShapeError                                Traceback (most recent call last)
<ipython-input-2-d7c5b13f43a6> in <module>
      7     pm.MutableData("x", np.zeros(10), dims="a")
      8 
----> 9 m.set_data("x", np.ones(5), coords={"a": np.arange(10, 15)})

~/repos/pymc3/pymc/model.py in set_data(self, name, values, coords)
   1164                         )
   1165                 if isinstance(length_tensor, TensorConstant):
-> 1166                     raise ShapeError(
   1167                         f"Resizing dimension '{dname}' is impossible, because "
   1168                         "a 'TensorConstant' stores its length. To be able "

ShapeError: Resizing dimension 'a' is impossible, because 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.

```
</details>

I think that it's rather strange not to be able to modify coordinates given via a constant tensor, but to be honest I didn't see when/where it was discussed.

## Versions and main components

* PyMC/PyMC3 Version: 4.0.0b6
* Aesara/Theano Version:
* Python Version:
* Operating system:
* How did you install PyMC/PyMC3: (conda/pip)


Discussion:
codecov[bot]: # [Codecov](https://codecov.io/gh/pymc-devs/pymc/pull/5817?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pymc-devs) Report
> Merging [#5817](https://codecov.io/gh/pymc-devs/pymc/pull/5817?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pymc-devs) (dd0a661) into [main](https://codecov.io/gh/pymc-devs/pymc/commit/4d2f3a8d5f99fc773bc0a1aad6d51e4f3f047698?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pymc-devs) (4d2f3a8) will **decrease** coverage by `11.69%`.
> The diff coverage is `28.57%`.

[![Impacted file tree graph](https://codecov.io/gh/pymc-devs/pymc/pull/5817/graphs/tree.svg?width=650&height=150&src=pr&token=JFuXtOJ4Cb&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pymc-devs)](https://codecov.io/gh/pymc-devs/pymc/pull/5817?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pymc-devs)

```diff
@@             Coverage Diff             @@
##             main    #5817       +/-   ##
===========================================
- Coverage   89.39%   77.70%   -11.70%     
===========================================
  Files          74       74               
  Lines       13767    13781       +14     
===========================================
- Hits        12307    10708     -1599     
- Misses       1460     3073     +1613     
```


| [Impacted Files](https://codecov.io/gh/pymc-devs/pymc/pull/5817?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pymc-devs) | Coverage Δ | |
|---|---|---|
| [pymc/model.py](https://codecov.io/gh/pymc-devs/pymc/pull/5817/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pymc-devs#diff-cHltYy9tb2RlbC5weQ==) | `69.75% <28.57%> (-16.78%)` | :arrow_down: |
| [pymc/model\_graph.py](https://codecov.io/gh/pymc-devs/pymc/pull/5817/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pymc-devs#diff-cHltYy9tb2RlbF9ncmFwaC5weQ==) | `17.60% <0.00%> (-76.00%)` | :arrow_down: |
| [pymc/func\_utils.py](https://codecov.io/gh/pymc-devs/pymc/pull/5817/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pymc-devs#diff-cHltYy9mdW5jX3V0aWxzLnB5) | `22.50% <0.00%> (-65.00%)` | :arrow_down: |
| [pymc/printing.py](https://codecov.io/gh/pymc-devs/pymc/pull/5817/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pymc-devs#diff-cHltYy9wcmludGluZy5weQ==) | `21.78% <0.00%> (-64.36%)` | :arrow_down: |
| [pymc/distributions/bound.py](https://codecov.io/gh/pymc-devs/pymc/pull/5817/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pymc-devs#diff-cHltYy9kaXN0cmlidXRpb25zL2JvdW5kLnB5) | `45.54% <0.00%> (-54.46%)` | :arrow_down: |
| [pymc/variational/updates.py](https://codecov.io/gh/pymc-devs/pymc/pull/5817/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pymc-devs#diff-cHltYy92YXJpYXRpb25hbC91cGRhdGVzLnB5) | `37.93% <0.00%> (-54.19%)` | :arrow_down: |
| [pymc/distributions/censored.py](https://codecov.io/gh/pymc-devs/pymc/pull/5817/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pymc-devs#diff-cHltYy9kaXN0cmlidXRpb25zL2NlbnNvcmVkLnB5) | `44.18% <0.00%> (-48.84%)` | :arrow_down: |
| [pymc/distributions/timeseries.py](https://codecov.io/gh/pymc-devs/pymc/pull/5817/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pymc-devs#diff-cHltYy9kaXN0cmlidXRpb25zL3RpbWVzZXJpZXMucHk=) | `30.63% <0.00%> (-46.84%)` | :arrow_down: |
| [pymc/exceptions.py](https://codecov.io/gh/pymc-devs/pymc/pull/5817/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pymc-devs#diff-cHltYy9leGNlcHRpb25zLnB5) | `58.62% <0.00%> (-41.38%)` | :arrow_down: |
| [pymc/step\_methods/hmc/quadpotential.py](https://codecov.io/gh/pymc-devs/pymc/pull/5817/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pymc-devs#diff-cHltYy9zdGVwX21ldGhvZHMvaG1jL3F1YWRwb3RlbnRpYWwucHk=) | `43.64% <0.00%> (-36.91%)` | :arrow_down: |
| ... and [28 more](https://codecov.io/gh/pymc-devs/pymc/pull/5817/diff?src=pr&el=tree-more&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pymc-devs) | |

AlexAndorra: Thanks for the proposed fix @michaelosthege ! I agree with @lucianopaz that it'd be great if the user didn't have to set_dim before they set_data
michaelosthege: I'll add the test suggested above and look into adding the `coords` kwarg to `set_data`.
michaelosthege: @lucianopaz please review the new test cases - they should govern which usage patterns are now supported.

The only thing where I'm hestiant is introducing mutable dims through the `pm.Model` constructor.
For that we could introduce a new kwarg, but I didn't want to make that decision already.
Options:
+ `pm.Model(coords_mutable=..., coords=...)`
+ `pm.Model(coords_mutable=..., coords_immutable=...)`
+ more?

Please continue on this branch as you see fit. Let's get it over the finish line today
ricardoV94: This is a vague hand waving comment, so apologies for that. I feel the internal code to manage dims is become a bit too convoluted for its own good. Suggestion:

1. Have a `coords` and `mutable_coords` kwargs in `pm.Model`. 
2. Remove defining news dims / coords with variables
3. Downsides?


michaelosthege: The idea of defining `coords` through `pm.Data` came from the `export_index_as_coords` kwarg.

Thinking this further, with #5796 we could do this:

```
population = xarray.DataArray(
    "population",
    [10, 200],
    dims="city",
    coords=dict(city=["Atown", "Bvil"]),
)

with pm.Model() as m:
    pop = pm.MutableData("population", population)
    assert "city" in m.coords # True
```

Note that this would work for multi-dimensional `xarray.DataArray` objects but not for `pd.DataFrame` with multi-index.

For wroking with GPs this would be really useful, because one has to carry data and grid coordinates around quite a lot..

-----

I agree that the code is a little convoluted - it would probably benefit from a little extraction refactoring.

But the downsides of removing `coords` from `pm.Data` would be that we'd have to break the `export_index_as_coords` feature and we couldn't do what the above example does.

If we add a kwarg to the `Model` signature, I would prefer `coords_mutable` because of how things are sorted by code completion tools.
michaelosthege: Any decision regarding the default `(im)mutable` of `pm.Model(coords)` will build on top of the commits made in this branch.

Since the commits here don't affect the thing we're still in discussion about, can we go ahead and merge this in the interest of closing #5812 and increasing overall safety?
twiecki: What's missing here, can we merge?

Patch:
diff --git a/docs/source/contributing/developer_guide_implementing_distribution.md b/docs/source/contributing/developer_guide_implementing_distribution.md
index 39920a389..32fa23640 100644
--- a/docs/source/contributing/developer_guide_implementing_distribution.md
+++ b/docs/source/contributing/developer_guide_implementing_distribution.md
@@ -1,10 +1,10 @@
 # Implementing a Distribution
 
-This guide provides an overview on how to implement a distribution for version 4 of PyMC.
+This guide provides an overview on how to implement a distribution for PyMC version `>=4.0.0`.
 It is designed for developers who wish to add a new distribution to the library.
-Users will not be aware of all this complexity and should instead make use of helper methods such as `~pymc.distributions.DensityDist`.
+Users will not be aware of all this complexity and should instead make use of helper methods such as `~pymc.DensityDist`.
 
-PyMC {class}`~pymc.distributions.Distribution` builds on top of Aesara's {class}`~aesara.tensor.random.op.RandomVariable`, and implements `logp`, `logcdf` and `moment` methods as well as other initialization and validation helpers.
+PyMC {class}`~pymc.Distribution` builds on top of Aesara's {class}`~aesara.tensor.random.op.RandomVariable`, and implements `logp`, `logcdf` and `moment` methods as well as other initialization and validation helpers.
 Most notably `shape/dims` kwargs, alternative parametrizations, and default `transforms`.
 
 Here is a summary check-list of the steps needed to implement a new distribution.
@@ -23,8 +23,8 @@ This guide does not attempt to explain the rationale behind the `Distributions`
 {class}`~aesara.tensor.random.op.RandomVariable` are responsible for implementing the random sampling methods, which in version 3 of PyMC used to be one of the standard `Distribution` methods, alongside `logp` and `logcdf`.
 The `RandomVariable` is also responsible for parameter broadcasting and shape inference.
 
-Before creating a new `RandomVariable` make sure that it is not offered in the {mod}`Numpy library <numpy.random>`.
-If it is, it should be added to the {doc}`Aesara library <aesara>` first and then imported into the PyMC library.
+Before creating a new `RandomVariable` make sure that it is not already offered in the {mod}`NumPy library <numpy.random>`.
+If it is, it should be added to the {doc}`Aesara library <aesara:index>` first and then imported into the PyMC library.
 
 In addition, it might not always be necessary to implement a new `RandomVariable`.
 For example if the new `Distribution` is just a special parametrization of an existing `Distribution`.
@@ -61,14 +61,14 @@ class BlahRV(RandomVariable):
 
     # If you want to add a custom signature and default values for the
     # parameters, do it like this. Otherwise this can be left out.
-    def __call__(self, loc=0.0, scale=1.0, size=None, **kwargs) -> TensorVariable:
-        return super().__call__(loc, scale, size=size, **kwargs)
+    def __call__(self, loc=0.0, scale=1.0, **kwargs) -> TensorVariable:
+        return super().__call__(loc, scale, **kwargs)
 
     # This is the Python code that produces samples.  Its signature will always
     # start with a NumPy `RandomState` object, then the distribution
     # parameters, and, finally, the size.
     #
-    # This is effectively the PyMC v4.x replacement for `Distribution.random`.
+    # This is effectively the PyMC >=4.0 replacement for `Distribution.random`.
     @classmethod
     def rng_fn(
         cls,
@@ -86,12 +86,12 @@ blah = BlahRV()
 
 Some important things to keep in mind:
 
-1. Everything inside the `rng_fn` method is pure Python code (as are the inputs) and should not make use of other `Aesara` symbolic ops. The random method should make use of the `rng` which is a Numpy    {class}`~numpy.random.RandomState`, so that samples are reproducible.
-1. The `size` argument (together with the inputs shape) are the only way for the user to specify non-default `RandomVariable` dimensions. The `rng_fn` will have to take this into consideration for correct output. `size` is the specification used by `Numpy` and `Scipy` and works like PyMC `shape` for univariate distributions, but is different for multivariate distributions. Unfortunately there is no general reference documenting how `size` ought to work for multivariate distributions. This [discussion](https://github.com/numpy/numpy/issues/17669) may be helpful to get more context.
-1. `Aesara` tries to infer the output shape of the `RandomVariable` (given a user-specified size) by introspection of the `ndim_supp` and `ndim_params` attributes. However, the default method may not work for more complex distributions. In that case, custom `_supp_shape_from_params` (and less probably, `_infer_shape`) should also be implemented in the new `RandomVariable` class. One simple example is seen in the {class}`~pymc.distributions.multivariate.DirichletMultinomialRV` where it was necessary to specify the `rep_param_idx` so that the `default_supp_shape_from_params` helper method could do its job. In more complex cases, it may not suffice to use this default helper. This could hapen for instance if the argument values determined the support shape of the distribution, as happens in the `~pymc.distribution.multivarite._LKJCholeskyCovRV`.
+1. Everything inside the `rng_fn` method is pure Python code (as are the inputs) and should not make use of other `Aesara` symbolic ops. The random method should make use of the `rng` which is a NumPy {class}`~numpy.random.RandomState`, so that samples are reproducible.
+1. Non-default `RandomVariable` dimensions will end up in the `rng_fn` via the `size` kwarg. The `rng_fn` will have to take this into consideration for correct output. `size` is the specification used by NumPy and SciPy and works like PyMC `shape` for univariate distributions, but is different for multivariate distributions. For multivariate distributions the __`size` excludes the `ndim_supp` support dimensions__, whereas the __`shape` of the resulting `TensorVariabe` or `ndarray` includes the support dimensions__. This [discussion](https://github.com/numpy/numpy/issues/17669) may be helpful to get more context.
+1. `Aesara` tries to infer the output shape of the `RandomVariable` (given a user-specified size) by introspection of the `ndim_supp` and `ndim_params` attributes. However, the default method may not work for more complex distributions. In that case, custom `_supp_shape_from_params` (and less probably, `_infer_shape`) should also be implemented in the new `RandomVariable` class. One simple example is seen in the {class}`~pymc.DirichletMultinomialRV` where it was necessary to specify the `rep_param_idx` so that the `default_supp_shape_from_params` helper method can do its job. In more complex cases, it may not suffice to use this default helper. This could happen for instance if the argument values determined the support shape of the distribution, as happens in the `~pymc.distributions.multivarite._LKJCholeskyCovRV`.
 1. It's okay to use the `rng_fn` `classmethods` of other Aesara and PyMC `RandomVariables` inside the new `rng_fn`. For example if you are implementing a negative HalfNormal `RandomVariable`, your `rng_fn` can simply return `- halfnormal.rng_fn(rng, scale, size)`.
 
-*Note: In addition to `size`, the PyMC API also provides `shape` and `dims` as alternatives to define a distribution dimensionality, but this is taken care of by {class}`~pymc.distributions.Distribution`, and should not require any extra changes.*
+*Note: In addition to `size`, the PyMC API also provides `shape` and `dims` as alternatives to define a distribution dimensionality, but this is taken care of by {class}`~pymc.Distribution`, and should not require any extra changes.*
 
 For a quick test that your new `RandomVariable` `Op` is working, you can call the `Op` with the necessary parameters and then call `eval()` on the returned object:
 
@@ -115,13 +115,13 @@ blah([0, 0], [1, 2], size=(10, 2)).eval()
 
 ## 2. Inheriting from a PyMC base `Distribution` class
 
-After implementing the new `RandomVariable` `Op`, it's time to make use of it in a new PyMC {class}`pymc.distributions.Distribution`.
-PyMC 4.x works in a very {term}`functional <Functional Programming>` way, and the `distribution` classes are there mostly to facilitate porting the `PyMC3` v3.x code to the new `PyMC` v4.x version, add PyMC API features and keep related methods organized together.
+After implementing the new `RandomVariable` `Op`, it's time to make use of it in a new PyMC {class}`~pymc.Distribution`.
+PyMC `>=4.0.0` works in a very {term}`functional <Functional Programming>` way, and the `distribution` classes are there mostly to facilitate porting the `PyMC3` v3.x code to PyMC `>=4.0.0`, add PyMC API features and keep related methods organized together.
 In practice, they take care of:
 
-1. Linking ({term}`Dispatching`) a rv_op class with the corresponding `moment`, `logp` and `logcdf` methods.
+1. Linking ({term}`Dispatching`) an rv_op class with the corresponding `moment`, `logp` and `logcdf` methods.
 1. Defining a standard transformation (for continuous distributions) that converts a bounded variable domain (e.g., positive line) to an unbounded domain (i.e., the real line), which many samplers prefer.
-1. Validating the parametrization of a distribution and converting non-symbolic inputs (i.e., numeric literals or numpy arrays) to symbolic variables.
+1. Validating the parametrization of a distribution and converting non-symbolic inputs (i.e., numeric literals or NumPy arrays) to symbolic variables.
 1. Converting multiple alternative parametrizations to the standard parametrization that the `RandomVariable` is defined in terms of.
 
 Here is how the example continues:
@@ -136,16 +136,16 @@ from pymc.distributions.dist_math import check_parameters
 # Subclassing `PositiveContinuous` will dispatch a default `log` transformation
 class Blah(PositiveContinuous):
     # This will be used by the metaclass `DistributionMeta` to dispatch the
-    # class `logp` and `logcdf` methods to the `blah` `op`
+    # class `logp` and `logcdf` methods to the `blah` `Op` defined in the last line of the code above.
     rv_op = blah
 
-    # dist() is responsible for returning an instance of the rv_op. We pass
-    # the standard parametrizations to super().dist
+    # dist() is responsible for returning an instance of the rv_op.
+    # We pass the standard parametrizations to super().dist
     @classmethod
     def dist(cls, param1, param2=None, alt_param2=None, **kwargs):
         param1 = at.as_tensor_variable(intX(param1))
         if param2 is not None and alt_param2 is not None:
-            raise ValueError('Only one of param2 and alt_param2 is allowed')
+            raise ValueError("Only one of param2 and alt_param2 is allowed.")
         if alt_param2 is not None:
             param2 = 1 / alt_param2
         param2 = at.as_tensor_variable(floatX(param2))
@@ -155,15 +155,16 @@ class Blah(PositiveContinuous):
         return super().dist([param1, param2], **kwargs)
 
     # moment returns a symbolic expression for the stable moment from which to start sampling
-    # the variable, given the implicit `rv`, `size` and `param1` ... `paramN`
+    # the variable, given the implicit `rv`, `size` and `param1` ... `paramN`.
+    # This is typically a "representative" point such as the the mean or mode.
     def moment(rv, size, param1, param2):
         moment, _ = at.broadcast_arrays(param1, param2)
         if not rv_size_is_none(size):
             moment = at.full(size, moment)
         return moment
 
-    # Logp returns a symbolic expression for the logp evaluation of the variable
-    # given the `value` of the variable and the parameters `param1` ... `paramN`
+    # Logp returns a symbolic expression for the elementwise log-pdf or log-pmf evaluation
+    # of the variable given the `value` of the variable and the parameters `param1` ... `paramN`.
     def logp(value, param1, param2):
         logp_expression = value * (param1 + at.log(param2))
 
@@ -175,14 +176,14 @@ class Blah(PositiveContinuous):
         )
 
         # We use `check_parameters` for parameter validation. After the default expression,
-        # multiple comma-separated symbolic conditions can be added. Whenever
-        # a bound is invalidated, the returned expression raises an error with the message
-        # defined in the optional `msg` keyword argument.
+        # multiple comma-separated symbolic conditions can be added.
+        # Whenever a bound is invalidated, the returned expression raises an error
+        # with the message defined in the optional `msg` keyword argument.
         return check_parameters(
             logp_expression,
             param2 >= 0,
             msg="param2 >= 0",
-    )
+        )
 
     # logcdf works the same way as logp. For bounded variables, it is expected to return
     # `-inf` for values below the domain start and `0` for values above the domain end.
@@ -193,12 +194,12 @@ class Blah(PositiveContinuous):
 
 Some notes:
 
-1. A distribution should at the very least inherit from {class}`~pymc.distributions.Discrete` or {class}`~pymc.distributions.Continuous`. For the latter, more specific subclasses exist: `PositiveContinuous`, `UnitContinuous`, `BoundedContinuous`, `CircularContinuous`, `SimplexContinuous`, which specify default transformations for the variables. If you need to specify a one-time custom transform you can also create a `_default_transform` dispatch function as is done for the {class}`~pymc.distributions.multivariate.LKJCholeskyCov`.
-1. If a distribution does not have a corresponding `random` implementation, a `RandomVariable` should still be created that raises a `NotImplementedError`. This is the case for the {class}`~pymc.distributions.continuous.Flat`. In this case it will be necessary to provide a `moment` method.
-1. As mentioned above, `PyMC` v4.x works in a very {term}`functional <Functional Programming>` way, and all the information that is needed in the `logp` and `logcdf` methods is expected to be "carried" via the `RandomVariable` inputs. You may pass numerical arguments that are not strictly needed for the `rng_fn` method but are used in the `logp` and `logcdf` methods. Just keep in mind whether this affects the correct shape inference behavior of the `RandomVariable`. If specialized non-numeric information is needed you might need to define your custom`_logp` and `_logcdf` {term}`Dispatching` functions, but this should be done as a last resort.
+1. A distribution should at the very least inherit from {class}`~pymc.Discrete` or {class}`~pymc.Continuous`. For the latter, more specific subclasses exist: `PositiveContinuous`, `UnitContinuous`, `BoundedContinuous`, `CircularContinuous`, `SimplexContinuous`, which specify default transformations for the variables. If you need to specify a one-time custom transform you can also create a `_default_transform` dispatch function as is done for the {class}`~pymc.distributions.multivariate.LKJCholeskyCov`.
+1. If a distribution does not have a corresponding `rng_fn` implementation, a `RandomVariable` should still be created to raise a `NotImplementedError`. This is, for example, the case in {class}`~pymc.distributions.continuous.Flat`. In this case it will be necessary to provide a `moment` method, because without a `rng_fn`, PyMC can't fall back to a random draw to use as an initial point for MCMC.
+1. As mentioned above, PyMC `>=4.0.0` works in a very {term}`functional <Functional Programming>` way, and all the information that is needed in the `logp` and `logcdf` methods is expected to be "carried" via the `RandomVariable` inputs. You may pass numerical arguments that are not strictly needed for the `rng_fn` method but are used in the `logp` and `logcdf` methods. Just keep in mind whether this affects the correct shape inference behavior of the `RandomVariable`. If specialized non-numeric information is needed you might need to define your custom`_logp` and `_logcdf` {term}`Dispatching` functions, but this should be done as a last resort.
 1. The `logcdf` method is not a requirement, but it's a nice plus!
-1. Currently only one moment is supported in the `moment` method, and probably the "higher-order" one is the most useful (that is `mean` > `median` > `mode`)... You might need to truncate the moment if you are dealing with a discrete distribution.
-1. When creating the `moment` method, we have to be careful with `size != None` and broadcast properly when some parameters that are not used in the moment may nevertheless inform about the shape of the distribution. E.g. `pm.Normal.dist(mu=0, sigma=np.arange(1, 6))` returns a moment of `[mu, mu, mu, mu, mu]`.
+1. Currently, only one moment is supported in the `moment` method, and probably the "higher-order" one is the most useful (that is `mean` > `median` > `mode`)... You might need to truncate the moment if you are dealing with a discrete distribution.
+1. When creating the `moment` method, be careful with `size != None` and broadcast properly also based on parameters that are not necessarily used to calculate the moment. For example, the `sigma` in `pm.Normal.dist(mu=0, sigma=np.arange(1, 6))` is irrelevant for the moment, but may nevertheless inform about the shape. In this case, the `moment` should return `[mu, mu, mu, mu, mu]`.
 
 For a quick check that things are working you can try the following:
 
@@ -292,11 +293,11 @@ You can see an example in {class}`~pymc.test.test_random.TestDirichlet`.
 ### Notes on `check_pymcs_draws_match_reference` test
 
 The `check_pymcs_draws_match_reference` is a very simple test for the equality of draws from the `RandomVariable` and the exact same python function, given the same inputs and random seed.
-A small number (`size=15`) is checked. This is not supposed to be a test for the correctness of the random generator.
-The latter kind of test (if warranted) can be performed with the aid of `pymc_random` and `pymc_random_discrete` methods in the same test file, which will perform an expensive statistical comparison between the RandomVariable `rng_fn` and a reference Python function.
-This kind of test only makes sense if there is a good independent generator reference (i.e., not just the same composition of numpy / scipy python calls that is done inside `rng_fn`).
+A small number (`size=15`) is checked. This is not supposed to be a test for the correctness of the random number generator.
+The latter kind of test (if warranted) can be performed with the aid of `pymc_random` and `pymc_random_discrete` methods in the same test file, which will perform an expensive statistical comparison between the `RandomVariable.rng_fn` and a reference Python function.
+This kind of test only makes sense if there is a good independent generator reference (i.e., not just the same composition of NumPy / SciPy calls that is done inside `rng_fn`).
 
-Finally, when your `rng_fn` is doing something more than just calling a `numpy` or `scipy` method, you will need to setup an equivalent seeded function with which to compare for the exact draws (instead of relying on `seeded_[scipy|numpy]_distribution_builder`).
+Finally, when your `rng_fn` is doing something more than just calling a NumPy or SciPy method, you will need to set up an equivalent seeded function with which to compare for the exact draws (instead of relying on `seeded_[scipy|numpy]_distribution_builder`).
 You can find an example in {class}`~pymc.tests.test_distributions_random.TestWeibull`, whose `rng_fn` returns `beta * np.random.weibull(alpha, size=size)`.
 
 
@@ -342,13 +343,13 @@ def test_blah(self):
 
 ```
 
-These methods will perform a grid evaluation on the combinations of domain and paramdomains values, and check that the pymc methods and the reference functions match.
+These methods will perform a grid evaluation on the combinations of `domain` and `paramdomains` values, and check that the PyMC methods and the reference functions match.
 There are a couple of details worth keeping in mind:
 
-1. By default the first and last values (edges) of the `Domain` are not compared (they are used for other things). If it is important to test the edge of the `Domain`, the edge values can be repeated. This is done by the `Bool`: `Bool = Domain([0, 0, 1, 1], "int64")`
-1. There are some default domains (such as `R` and `Rplus`) that you can use for testing your new distribution, but it's also perfectly fine to create your own domains inside the test function if there is a good reason for it (e.g., when the default values lead too many extreme unlikely combinations that are not very informative about the correctness of the implementation).
-1. By default, a random subset of 100 `param` x `paramdomain` combinations is tested, in order to keep the test runtime under control. When testing your shiny new distribution, you can temporarily set `n_samples=-1` to force all combinations to be tested. This is important to avoid the your `PR` leading to surprising failures in future runs whenever some bad combinations of parameters are randomly tested.
-1. On Github all tests are run twice, under the `aesara.config.floatX` flags of `"float64"` and `"float32"`. However, the reference Python functions will run in a pure "float64" environment, which means the reference and the PyMC results can diverge quite a lot (e.g., underflowing to `-np.inf` for extreme parameters). You should therefore make sure you test locally in both regimes. A quick and dirty way of doing this is to temporariliy add `aesara.config.floatX = "float32"` at the very top of file, immediately after `import aesara`. Remember to set `n_samples=-1` as well to test all combinations. The test output will show what exact parameter values lead to a failure. If you are confident that your implementation is correct, you may opt to tweak the decimal precision with `select_by_precision`, or adjust the tested `Domain` values. In extreme cases, you can mark the test with a conditional `xfail` (if only one of the sub-methods is failing, they should be separated, so that the `xfail` is as narrow as possible):
+1. By default, the first and last values (edges) of the `Domain` are not compared (they are used for other things). If it is important to test the edge of the `Domain`, the edge values can be repeated. This is done by the `Bool`: `Bool = Domain([0, 0, 1, 1], "int64")`
+3. There are some default domains (such as `R` and `Rplus`) that you can use for testing your new distribution, but it's also perfectly fine to create your own domains inside the test function if there is a good reason for it (e.g., when the default values lead too many extreme unlikely combinations that are not very informative about the correctness of the implementation).
+4. By default, a random subset of 100 `param` x `paramdomain` combinations is tested, to keep the test runtime under control. When testing your shiny new distribution, you can temporarily set `n_samples=-1` to force all combinations to be tested. This is important to avoid your `PR` leading to surprising failures in future runs whenever some bad combinations of parameters are randomly tested.
+5. On GitHub some tests run twice, under the `aesara.config.floatX` flags of `"float64"` and `"float32"`. However, the reference Python functions will run in a pure "float64" environment, which means the reference and the PyMC results can diverge quite a lot (e.g., underflowing to `-np.inf` for extreme parameters). You should therefore make sure you test locally in both regimes. A quick and dirty way of doing this is to temporarily add `aesara.config.floatX = "float32"` at the very top of file, immediately after `import aesara`. Remember to set `n_samples=-1` as well to test all combinations. The test output will show what exact parameter values lead to a failure. If you are confident that your implementation is correct, you may opt to tweak the decimal precision with `select_by_precision`, or adjust the tested `Domain` values. In extreme cases, you can mark the test with a conditional `xfail` (if only one of the sub-methods is failing, they should be separated, so that the `xfail` is as narrow as possible):
 
 ```python
 
diff --git a/docs/source/contributing/python_style.md b/docs/source/contributing/python_style.md
index b707d07a5..bf8a9b289 100644
--- a/docs/source/contributing/python_style.md
+++ b/docs/source/contributing/python_style.md
@@ -47,3 +47,31 @@ or, if you just want to manually run them on a subset of files,
 ```bash
 pre-commit run --files <file_1> <file_2> ... <file_n>
 ```
+
+## Gotchas & Troubleshooting
+__Pre-commit runs on staged files__
+
+If you have some `git` changes staged and other unstaged, the `pre-commit` will only run on the staged files.
+
+__Pre-commit repeatedly complains about the same formatting changes__
+
+Check the unstaged changes (see previous point).
+
+__Whitespace changes in the `environment-dev-*.yml` files__
+
+On Windows, there are some bugs in pre-commit hooks that can lead to changes in some environment YAML files.
+Until this is fixed upstream, you should __ignore these changes__.
+To actually make the commit, deactivate the automated `pre-commit` with `pre-commit uninstall` and make sure to run it manually with `pre-commit run --all`.
+
+__Failures in the `mypy` step__
+
+We are running static type checks with `mypy` to continuously improve the reliability and type safety of the PyMC codebase.
+However, there are many files with unresolved type problems, which is why we are allowing some files to fail the `mypy` check.
+
+If you are seeing the `mypy` step complain, chances are that you are in one of the following two situations:
+* 😕 Your changes introduced type problems in a file that was previously free of type problems.
+* 🥳 Your changes fixed type problems.
+
+In any case __read the logging output of the `mypy` hook__, because it contains the instructions how to proceed.
+
+You can also run the `mypy` check manually with `python scripts/run_mypy.py [--verbose]`.
diff --git a/docs/source/contributing/review_pr_pymc_examples.md b/docs/source/contributing/review_pr_pymc_examples.md
index 215b4c2dc..caf1cb557 100644
--- a/docs/source/contributing/review_pr_pymc_examples.md
+++ b/docs/source/contributing/review_pr_pymc_examples.md
@@ -14,7 +14,7 @@ The most important guideline is the following: **When you aren't _completely_
 sure about something communicate it, ask around and (re)read the documentation**.
 :::
 
-pymc-examples is a huge collection of notebooks, covering many fields
+[pymc-examples](https://github.com/pymc-devs/pymc-examples/) is a huge collection of notebooks, covering many fields
 and applications of PyMC and Bayesian statistics. In addition,
 its HTML generation infrastructure has been improved by specialized
 people who were also paid so they could dedicate themselves to this task.
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)


Test patch:
diff --git a/docs/source/contributing/running_the_test_suite.md b/docs/source/contributing/running_the_test_suite.md
index 60aeb8608..e415b1d71 100644
--- a/docs/source/contributing/running_the_test_suite.md
+++ b/docs/source/contributing/running_the_test_suite.md
@@ -1,15 +1,47 @@
 (running_the_test_suite)=
 # Running the test suite
+The first step to run tests is the installation of additional dependencies that are needed for testing:
 
-TODO: explain steps to run test suite. This is a how to guide, so assume readers know how to install
-things and so on, at most mention what dependencies are needed.
+```bash
+pip install -r requirements-dev.txt
+```
+
+The PyMC test suite uses `pytest` as the testing framework.
+If you are unfamiliar with `pytest`, check out [this short video series](https://calmcode.io/pytest/).
+
+With the optional dependencies installed, you can start running tests.
+Below are some example of how you might want to run certain parts of the test suite.
+
+```{attention}
+Running the entire test suite will take hours.
+Therefore, we recommend to run just specific tests that target the parts of the codebase you're working on.
+```
+
+To run all tests from a single file:
+```bash
+pytest -v pymc/tests/test_model.py
+```
+
+```{tip}
+The `-v` flag is short-hand for `--verbose` and prints the names of the test cases that are currently running.
+```
+
+Often, you'll want to focus on just a few test cases first.
+By using the `-k` flag, you can filter for test cases that match a certain pattern.
+For example, the following command runs all test cases from `test_model.py` that have "coord" in their name:
 
 ```bash
-pip install pytest pytest-cov coverage
+pytest -v pymc/tests/test_model.py -k coord
+```
+
 
-# To run a subset of tests
-pytest --verbose pymc/tests/<name of test>.py
+To get a coverage report, you can pass `--cov=pymc`, optionally with `--cov-report term-missing` to get a printout of the line numbers that were visited by the invoked tests.
+Note that because you are not running the entire test suite, the coverage will be terrible.
+But you can still watch for specific line numbers of the code that you're working on.
 
-# To get a coverage report
-pytest --verbose --cov=pymc --cov-report term-missing pymc/tests/<name of test>.py
+```bash
+pytest -v --cov=pymc --cov-report term-missing pymc/tests/<name of test>.py
 ```
+
+When you are reasonably confident about the changes you made, you can push the changes and open a pull request.
+Our GitHub Actions pipeline will run the entire test suite and if there are failures you can go back and run these tests on your local machine.
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: