diff --git a/RELEASE-NOTES.md b/RELEASE-NOTES.md
index 7a73c65b72..b7918cb30c 100644
--- a/RELEASE-NOTES.md
+++ b/RELEASE-NOTES.md
@@ -8,11 +8,18 @@
 
 ### New Features
 - 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 [#4625](https://github.com/pymc-devs/pymc3/pull/4625)):
+  - With `shape` the length of dimensions must be given numerically or as scalar Aesara `Variables`. A `SpecifyShape` `Op` is added automatically unless `Ellipsis` is used. Using `shape` restricts 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 creates new dimensions in addition to what is implied by RV parameters.
+  - An `Ellipsis` (`...`) in the last position of `shape` or `dims` can be used as short-hand notation for implied dimensions.
 - ...
 
 ### Maintenance
 - Remove float128 dtype support (see [#4514](https://github.com/pymc-devs/pymc3/pull/4514)).
 - Logp method of `Uniform` and `DiscreteUniform` no longer depends on `pymc3.distributions.dist_math.bound` for proper evaluation (see [#4541](https://github.com/pymc-devs/pymc3/pull/4541)).
+- `Model.RV_dims` and `Model.coords` are now read-only properties. To modify the `coords` dictionary use `Model.add_coord`. Also `dims` or coordinate values that are `None` will be auto-completed (see [#4625](https://github.com/pymc-devs/pymc3/pull/4625)).
+- The length of `dims` in the model is now tracked symbolically through `Model.dim_lengths` (see [#4625](https://github.com/pymc-devs/pymc3/pull/4625)).
 - ...
 
 ## PyMC3 3.11.2 (14 March 2021)
diff --git a/pymc3/aesaraf.py b/pymc3/aesaraf.py
index e30248b841..5a99849fd2 100644
--- a/pymc3/aesaraf.py
+++ b/pymc3/aesaraf.py
@@ -45,6 +45,7 @@
 from aesara.sandbox.rng_mrg import MRG_RandomStream as RandomStream
 from aesara.tensor.elemwise import Elemwise
 from aesara.tensor.random.op import RandomVariable
+from aesara.tensor.shape import SpecifyShape
 from aesara.tensor.sharedvar import SharedVariable
 from aesara.tensor.subtensor import AdvancedIncSubtensor, AdvancedIncSubtensor1
 from aesara.tensor.var import TensorVariable
@@ -146,6 +147,8 @@ def change_rv_size(
         Expand the existing size by `new_size`.
 
     """
+    if isinstance(rv_var.owner.op, SpecifyShape):
+        rv_var = rv_var.owner.inputs[0]
     rv_node = rv_var.owner
     rng, size, dtype, *dist_params = rv_node.inputs
     name = rv_var.name
diff --git a/pymc3/backends/arviz.py b/pymc3/backends/arviz.py
index 97a0c13613..8a3f7b46cc 100644
--- a/pymc3/backends/arviz.py
+++ b/pymc3/backends/arviz.py
@@ -162,10 +162,7 @@ def __init__(
         self.trace = trace
 
         # this permits us to get the model from command-line argument or from with model:
-        try:
-            self.model = modelcontext(model)
-        except TypeError:
-            self.model = None
+        self.model = modelcontext(model)
 
         self.attrs = None
         if trace is not None:
@@ -223,10 +220,14 @@ def arbitrary_element(dct: Dict[Any, np.ndarray]) -> np.ndarray:
         self.coords = {} if coords is None else coords
         if hasattr(self.model, "coords"):
             self.coords = {**self.model.coords, **self.coords}
+        self.coords = {key: value for key, value in self.coords.items() if value is not None}
 
         self.dims = {} if dims is None else dims
         if hasattr(self.model, "RV_dims"):
-            model_dims = {k: list(v) for k, v in self.model.RV_dims.items()}
+            model_dims = {
+                var_name: [dim for dim in dims if dim is not None]
+                for var_name, dims in self.model.RV_dims.items()
+            }
             self.dims = {**model_dims, **self.dims}
 
         self.density_dist_obs = density_dist_obs
diff --git a/pymc3/data.py b/pymc3/data.py
index 846e8272b7..06dfb2766b 100644
--- a/pymc3/data.py
+++ b/pymc3/data.py
@@ -19,7 +19,7 @@
 import urllib.request
 
 from copy import copy
-from typing import Any, Dict, List
+from typing import Any, Dict, List, Sequence
 
 import aesara
 import aesara.tensor as at
@@ -502,7 +502,7 @@ class Data:
     >>> for data_vals in observed_data:
     ...     with model:
     ...         # Switch out the observed dataset
-    ...         pm.set_data({'data': data_vals})
+    ...         model.set_data('data', data_vals)
     ...         traces.append(pm.sample())
 
     To set the value of the data container variable, check out
@@ -543,6 +543,11 @@ def __new__(self, name, value, *, dims=None, export_index_as_coords=False):
 
         if export_index_as_coords:
             model.add_coords(coords)
+        elif dims:
+            # Register new dimension lengths
+            for d, dname in enumerate(dims):
+                if not dname in model.dim_lengths:
+                    model.add_coord(dname, values=None, length=shared_object.shape[d])
 
         # To draw the node for this variable in the graphviz Digraph we need
         # its shape.
@@ -562,7 +567,7 @@ def __new__(self, name, value, *, dims=None, export_index_as_coords=False):
         return shared_object
 
     @staticmethod
-    def set_coords(model, value, dims=None):
+    def set_coords(model, value, dims=None) -> Dict[str, Sequence]:
         coords = {}
 
         # If value is a df or a series, we interpret the index as coords:
diff --git a/pymc3/distributions/distribution.py b/pymc3/distributions/distribution.py
index 4960592c8c..3fba3cf7a2 100644
--- a/pymc3/distributions/distribution.py
+++ b/pymc3/distributions/distribution.py
@@ -20,21 +20,18 @@
 
 from abc import ABCMeta
 from copy import copy
-from typing import TYPE_CHECKING
+from typing import Any, Optional, Sequence, Tuple, Union
 
+import aesara
+import aesara.tensor as at
 import dill
 
+from aesara.graph.basic import Variable
 from aesara.tensor.random.op import RandomVariable
+from aesara.tensor.shape import SpecifyShape, specify_shape
 
+from pymc3.aesaraf import change_rv_size, pandas_to_array
 from pymc3.distributions import _logcdf, _logp
-
-if TYPE_CHECKING:
-    from typing import Optional, Callable
-
-import aesara
-import aesara.graph.basic
-import aesara.tensor as at
-
 from pymc3.util import UNSET, get_repr_for_variable
 from pymc3.vartypes import string_types
 
@@ -52,6 +49,10 @@
 
 PLATFORM = sys.platform
 
+Shape = Union[int, Sequence[Union[str, type(Ellipsis)]], Variable]
+Dims = Union[str, Sequence[Union[str, None, type(Ellipsis)]]]
+Size = Union[int, Tuple[int, ...]]
+
 
 class _Unpickling:
     pass
@@ -122,13 +123,111 @@ def logcdf(op, var, rvs_to_values, *dist_params, **kwargs):
         return new_cls
 
 
+def _valid_ellipsis_position(items: Union[None, Shape, Dims, Size]) -> bool:
+    if items is not None and not isinstance(items, Variable) and Ellipsis in items:
+        if any(i == Ellipsis for i in items[:-1]):
+            return False
+    return True
+
+
+def _validate_shape_dims_size(
+    shape: Any = None, dims: Any = None, size: Any = None
+) -> Tuple[Optional[Shape], Optional[Dims], Optional[Size]]:
+    # Raise on unsupported parametrization
+    if shape is not None and dims is not None:
+        raise ValueError("Passing both `shape` ({shape}) and `dims` ({dims}) is not supported!")
+    if dims is not None and size is not None:
+        raise ValueError("Passing both `dims` ({dims}) and `size` ({size}) is not supported!")
+    if shape is not None and size is not None:
+        raise ValueError("Passing both `shape` ({shape}) and `size` ({size}) is not supported!")
+
+    # Raise on invalid types
+    if not isinstance(shape, (type(None), int, list, tuple, Variable)):
+        raise ValueError("The `shape` parameter must be an int, list or tuple.")
+    if not isinstance(dims, (type(None), str, list, tuple)):
+        raise ValueError("The `dims` parameter must be a str, list or tuple.")
+    if not isinstance(size, (type(None), int, list, tuple)):
+        raise ValueError("The `size` parameter must be an int, list or tuple.")
+
+    # Auto-convert non-tupled parameters
+    if isinstance(shape, int):
+        shape = (shape,)
+    if isinstance(dims, str):
+        dims = (dims,)
+    if isinstance(size, int):
+        size = (size,)
+
+    # Convert to actual tuples
+    if not isinstance(shape, (type(None), tuple, Variable)):
+        shape = tuple(shape)
+    if not isinstance(dims, (type(None), tuple)):
+        dims = tuple(dims)
+    if not isinstance(size, (type(None), tuple)):
+        size = tuple(size)
+
+    if not _valid_ellipsis_position(shape):
+        raise ValueError(
+            f"Ellipsis in `shape` may only appear in the last position. Actual: {shape}"
+        )
+    if not _valid_ellipsis_position(dims):
+        raise ValueError(f"Ellipsis in `dims` may only appear in the last position. Actual: {dims}")
+    if size is not None and Ellipsis in size:
+        raise ValueError("The `size` parameter cannot contain an Ellipsis. Actual: {size}")
+    return shape, dims, size
+
+
 class Distribution(metaclass=DistributionMeta):
     """Statistical distribution"""
 
     rv_class = None
     rv_op = None
 
-    def __new__(cls, name, *args, **kwargs):
+    def __new__(
+        cls,
+        name: str,
+        *args,
+        rng=None,
+        dims: Optional[Dims] = None,
+        testval=None,
+        observed=None,
+        total_size=None,
+        transform=UNSET,
+        **kwargs,
+    ) -> RandomVariable:
+        """Adds a RandomVariable corresponding to a PyMC3 distribution to the current model.
+
+        Note that all remaining kwargs must be compatible with ``.dist()``
+
+        Parameters
+        ----------
+        cls : type
+            A PyMC3 distribution.
+        name : str
+            Name for the new model variable.
+        rng : optional
+            Random number generator to use with the RandomVariable.
+        dims : tuple, optional
+            A tuple of dimension names known to the model.
+        testval : optional
+            Test value to be attached to the output RV.
+            Must match its shape exactly.
+        observed : optional
+            Observed data to be passed when registering the random variable in the model.
+            See ``Model.register_rv``.
+        total_size : float, optional
+            See ``Model.register_rv``.
+        transform : optional
+            See ``Model.register_rv``.
+        **kwargs
+            Keyword arguments that will be forwarded to ``.dist()``.
+            Most prominently: ``shape`` and ``size``
+
+        Returns
+        -------
+        rv : RandomVariable
+            The created RV, registered in the Model.
+        """
+
         try:
             from pymc3.model import Model
 
@@ -141,40 +240,150 @@ def __new__(cls, name, *args, **kwargs):
                 "for a standalone distribution."
             )
 
-        rng = kwargs.pop("rng", None)
-
-        if rng is None:
-            rng = model.default_rng
-
         if not isinstance(name, string_types):
             raise TypeError(f"Name needs to be a string but got: {name}")
 
-        data = kwargs.pop("observed", None)
+        if rng is None:
+            rng = model.default_rng
 
-        total_size = kwargs.pop("total_size", None)
+        _, dims, _ = _validate_shape_dims_size(dims=dims)
+        resize = None
+
+        # Create the RV without specifying testval, because the testval may have a shape
+        # that only matches after replicating with a size implied by dims (see below).
+        rv_out = cls.dist(*args, rng=rng, testval=None, **kwargs)
+        n_implied = rv_out.ndim
+
+        # The `.dist()` can wrap automatically with a SpecifyShape Op which brings informative
+        # error messages earlier in model construction.
+        # Here, however, the underyling RV must be used - a new SpecifyShape Op can be added at the end.
+        assert_shape = None
+        if isinstance(rv_out.owner.op, SpecifyShape):
+            rv_out, assert_shape = rv_out.owner.inputs
+
+        # `dims` are only available with this API, because `.dist()` can be used
+        # without a modelcontext and dims are not tracked at the Aesara level.
+        if dims is not None:
+            if Ellipsis in dims:
+                # Auto-complete the dims tuple to the full length
+                dims = (*dims[:-1], *[None] * rv_out.ndim)
+
+            n_resize = len(dims) - n_implied
+
+            # All resize dims must be known already (numerically or symbolically).
+            unknown_resize_dims = set(dims[:n_resize]) - set(model.dim_lengths)
+            if unknown_resize_dims:
+                raise KeyError(
+                    f"Dimensions {unknown_resize_dims} are unknown to the model and cannot be used to specify a `size`."
+                )
 
-        dims = kwargs.pop("dims", None)
+            # The numeric/symbolic resize tuple can be created using model.RV_dim_lengths
+            resize = tuple(model.dim_lengths[dname] for dname in dims[:n_resize])
+        elif observed is not None:
+            if not hasattr(observed, "shape"):
+                observed = pandas_to_array(observed)
+            n_resize = observed.ndim - n_implied
+            resize = tuple(observed.shape[d] for d in range(n_resize))
+
+        if resize:
+            # A batch size was specified through `dims`, or implied by `observed`.
+            rv_out = change_rv_size(rv_var=rv_out, new_size=resize, expand=True)
+
+        if dims is not None:
+            # Now that we have a handle on the output RV, we can register named implied dimensions that
+            # were not yet known to the model, such that they can be used for size further downstream.
+            for di, dname in enumerate(dims[n_resize:]):
+                if not dname in model.dim_lengths:
+                    model.add_coord(dname, values=None, length=rv_out.shape[n_resize + di])
 
-        if "shape" in kwargs:
-            raise DeprecationWarning("The `shape` keyword is deprecated; use `size`.")
+        if testval is not None:
+            # Assigning the testval earlier causes trouble because the RV may not be created with the final shape already.
+            rv_out.tag.test_value = testval
 
-        transform = kwargs.pop("transform", UNSET)
+        rv_registered = model.register_rv(
+            rv_out, name, observed, total_size, dims=dims, transform=transform
+        )
 
-        rv_out = cls.dist(*args, rng=rng, **kwargs)
+        # Wrapping in specify_shape now does not break transforms:
+        if assert_shape is not None:
+            rv_registered = specify_shape(rv_registered, assert_shape)
 
-        return model.register_rv(rv_out, name, data, total_size, dims=dims, transform=transform)
+        return rv_registered
 
     @classmethod
-    def dist(cls, dist_params, **kwargs):
+    def dist(
+        cls,
+        dist_params,
+        *,
+        shape: Optional[Shape] = None,
+        size: Optional[Size] = None,
+        testval=None,
+        **kwargs,
+    ) -> RandomVariable:
+        """Creates a RandomVariable corresponding to the `cls` distribution.
+
+        Parameters
+        ----------
+        dist_params
+        shape : tuple, optional
+            A tuple of sizes for each dimension of the new RV.
+
+            Ellipsis (...) may be used in the last position of the tuple,
+            and automatically expand to the shape implied by RV inputs.
+
+            Without Ellipsis, a `SpecifyShape` Op is automatically applied,
+            constraining this model variable to exactly the specified shape.
+        size : int, tuple, Variable, optional
+            A scalar or tuple for replicating the RV in addition
+            to its implied shape/dimensionality.
+        testval : optional
+            Test value to be attached to the output RV.
+            Must match its shape exactly.
 
-        testval = kwargs.pop("testval", None)
+        Returns
+        -------
+        rv : RandomVariable
+            The created RV.
+        """
+        if "dims" in kwargs:
+            raise NotImplementedError("The use of a `.dist(dims=...)` API is not yet supported.")
+
+        shape, _, size = _validate_shape_dims_size(shape=shape, size=size)
+        assert_shape = None
+
+        # Create the RV without specifying size or testval.
+        # The size will be expanded later (if necessary) and only then the testval fits.
+        rv_native = cls.rv_op(*dist_params, size=None, **kwargs)
+
+        if shape is None and size is None:
+            size = ()
+        elif shape is not None:
+            # SpecifyShape is automatically applied for symbolic and non-Ellipsis shapes
+            if isinstance(shape, Variable):
+                assert_shape = shape
+                size = ()
+            else:
+                if Ellipsis in shape:
+                    size = tuple(shape[:-1])
+                else:
+                    size = tuple(shape[: len(shape) - rv_native.ndim])
+                    assert_shape = shape
+        # no-op conditions:
+        # `elif size is not None` (User already specified how to expand the RV)
+        # `else` (Unreachable)
+
+        if size:
+            rv_out = change_rv_size(rv_var=rv_native, new_size=size, expand=True)
+        else:
+            rv_out = rv_native
 
-        rv_var = cls.rv_op(*dist_params, **kwargs)
+        if assert_shape is not None:
+            rv_out = specify_shape(rv_out, shape=assert_shape)
 
         if testval is not None:
-            rv_var.tag.test_value = testval
+            rv_out.tag.test_value = testval
 
-        return rv_var
+        return rv_out
 
     def _distr_parameters_for_repr(self):
         """Return the names of the parameters for this distribution (e.g. "mu"
diff --git a/pymc3/model.py b/pymc3/model.py
index b4478108c0..77c3e18899 100644
--- a/pymc3/model.py
+++ b/pymc3/model.py
@@ -18,10 +18,20 @@
 import warnings
 
 from sys import modules
-from typing import TYPE_CHECKING, Any, Dict, List, Optional, Type, TypeVar, Union
+from typing import (
+    TYPE_CHECKING,
+    Any,
+    Dict,
+    List,
+    Optional,
+    Sequence,
+    Tuple,
+    Type,
+    TypeVar,
+    Union,
+)
 
 import aesara
-import aesara.graph.basic
 import aesara.sparse as sparse
 import aesara.tensor as at
 import numpy as np
@@ -32,11 +42,11 @@
 from aesara.graph.basic import Constant, Variable, graph_inputs
 from aesara.graph.fg import FunctionGraph, MissingInputError
 from aesara.tensor.random.opt import local_subtensor_rv_lift
+from aesara.tensor.sharedvar import ScalarSharedVariable
 from aesara.tensor.var import TensorVariable
 from pandas import Series
 
 from pymc3.aesaraf import (
-    change_rv_size,
     gradient,
     hessian,
     inputvars,
@@ -46,7 +56,7 @@
 from pymc3.blocking import DictToArrayBijection, RaveledVars
 from pymc3.data import GenTensorVariable, Minibatch
 from pymc3.distributions import logp_transform, logpt, logpt_sum
-from pymc3.exceptions import ImputationWarning, SamplingError
+from pymc3.exceptions import ImputationWarning, SamplingError, ShapeError
 from pymc3.math import flatten_list
 from pymc3.util import UNSET, WithMemoization, get_var_name, treedict, treelist
 from pymc3.vartypes import continuous_types, discrete_types, typefilter
@@ -606,8 +616,9 @@ def __new__(cls, *args, **kwargs):
 
     def __init__(self, name="", model=None, aesara_config=None, coords=None, check_bounds=True):
         self.name = name
-        self.coords = {}
-        self.RV_dims = {}
+        self._coords = {}
+        self._RV_dims = {}
+        self._dim_lengths = {}
         self.add_coords(coords)
         self.check_bounds = check_bounds
 
@@ -826,6 +837,27 @@ def basic_RVs(self):
         """
         return self.free_RVs + self.observed_RVs
 
+    @property
+    def RV_dims(self) -> Dict[str, Tuple[Union[str, None], ...]]:
+        """Tuples of dimension names for specific model variables.
+
+        Entries in the tuples may be ``None``, if the RV dimension was not given a name.
+        """
+        return self._RV_dims
+
+    @property
+    def coords(self) -> Dict[str, Union[Sequence, None]]:
+        """Coordinate values for model dimensions."""
+        return self._coords
+
+    @property
+    def dim_lengths(self) -> Dict[str, Tuple[Variable, ...]]:
+        """The symbolic lengths of dimensions in the model.
+
+        The values are typically instances of ``TensorVariable`` or ``ScalarSharedVariable``.
+        """
+        return self._dim_lengths
+
     @property
     def unobserved_RVs(self):
         """List of all random variables, including deterministic ones.
@@ -913,20 +945,144 @@ def shape_from_dims(self, dims):
             shape.extend(np.shape(self.coords[dim]))
         return tuple(shape)
 
-    def add_coords(self, coords):
+    def add_coord(
+        self,
+        name: str,
+        values: Optional[Sequence] = None,
+        *,
+        length: Optional[Variable] = None,
+    ):
+        """Registers a dimension coordinate with the model.
+
+        Parameters
+        ----------
+        name : str
+            Name of the dimension.
+            Forbidden: {"chain", "draw", "__sample__"}
+        values : optional, array-like
+            Coordinate values or ``None`` (for auto-numbering).
+            If ``None`` is passed, a ``length`` must be specified.
+        length : optional, scalar
+            A symbolic scalar of the dimensions length.
+            Defaults to ``aesara.shared(len(values))``.
+        """
+        if name in {"draw", "chain", "__sample__"}:
+            raise ValueError(
+                "Dimensions can not be named `draw`, `chain` or `__sample__`, "
+                "as those are reserved for use in `InferenceData`."
+            )
+        if values is None and length is None:
+            raise ValueError(
+                f"Either `values` or `length` must be specified for the '{name}' dimension."
+            )
+        if 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 name in self.coords:
+            if not values.equals(self.coords[name]):
+                raise ValueError(f"Duplicate and incompatiple coordinate: {name}.")
+        else:
+            self._coords[name] = values
+            self._dim_lengths[name] = length or aesara.shared(len(values))
+
+    def add_coords(
+        self,
+        coords: Dict[str, Optional[Sequence]],
+        *,
+        lengths: Optional[Dict[str, Union[Variable, None]]] = None,
+    ):
+        """Vectorized version of ``Model.add_coord``."""
         if coords is None:
             return
+        lengths = lengths or {}
 
-        for name in coords:
-            if name in {"draw", "chain"}:
-                raise ValueError(
-                    "Dimensions can not be named `draw` or `chain`, as they are reserved for the sampler's outputs."
+        for name, values in coords.items():
+            self.add_coord(name, values, length=lengths.get(name, None))
+
+    def set_data(
+        self,
+        name: str,
+        values: Dict[str, Optional[Sequence]],
+        coords: Optional[Dict[str, Sequence]] = None,
+    ):
+        """Changes the values of a data variable in the model.
+
+        In contrast to pm.Data().set_value, this method can also
+        update the corresponding coordinates.
+
+        Parameters
+        ----------
+        name : str
+            Name of a shared variable in the model.
+        values : array-like
+            New values for the shared variable.
+        coords : optional, dict
+            New coordinate values for dimensions of the shared variable.
+            Must be provided for all named dimensions that change in length
+            and already have coordinate values.
+        """
+        shared_object = self[name]
+        if not isinstance(shared_object, SharedVariable):
+            raise TypeError(
+                f"The variable `{name}` must be a `SharedVariable` (e.g. `pymc3.Data`) allow updating. "
+                f"The current type is: {type(shared_object)}"
+            )
+
+        if isinstance(values, list):
+            values = np.array(values)
+        values = pandas_to_array(values)
+        dims = self.RV_dims.get(name, None) or ()
+        coords = coords or {}
+
+        if values.ndim != shared_object.ndim:
+            raise ValueError(
+                f"New values for '{name}' must have {shared_object.ndim} dimensions, just like the original."
+            )
+
+        for d, dname in enumerate(dims):
+            length_tensor = self.dim_lengths[dname]
+            old_length = length_tensor.eval()
+            new_length = values.shape[d]
+            original_coords = self.coords.get(dname, None)
+            new_coords = coords.get(dname, None)
+
+            length_changed = new_length != old_length
+
+            # Reject resizing if we already know that it would create shape problems.
+            # NOTE: If there are multiple pm.Data containers sharing this dim, but the user only
+            #       changes the values for one of them, they will run into shape problems nonetheless.
+            length_belongs_to = length_tensor.owner.inputs[0].owner.inputs[0]
+            if not isinstance(length_belongs_to, SharedVariable) and length_changed:
+                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. "
+                    f"Check if the dimension was defined implicitly before the shared variable '{name}' was created, "
+                    f"for example by a model variable.",
+                    actual=new_length,
+                    expected=old_length,
                 )
-            if name in self.coords:
-                if not coords[name].equals(self.coords[name]):
-                    raise ValueError("Duplicate and incompatiple coordinate: %s." % name)
-            else:
-                self.coords[name] = coords[name]
+            if original_coords is not None and length_changed:
+                if length_changed and new_coords is None:
+                    raise ValueError(
+                        f"The '{name}' variable already had {len(original_coords)} coord values defined for"
+                        f"its {dname} dimension. With the new values this dimension changes to length "
+                        f"{new_length}, so new coord values for the {dname} dimension are required."
+                    )
+            if new_coords is not None:
+                # Update the registered coord values (also if they were None)
+                if len(new_coords) != new_length:
+                    raise ShapeError(
+                        f"Length of new coordinate values for dimension '{dname}' does not match the provided values.",
+                        actual=len(new_coords),
+                        expected=new_length,
+                    )
+                self._coords[dname] = new_coords
+            if isinstance(length_tensor, ScalarSharedVariable) and new_length != old_length:
+                # Updating the shared variable resizes dependent nodes that use this dimension for their `size`.
+                length_tensor.set_value(new_length)
+
+        shared_object.set_value(values)
 
     def register_rv(self, rv_var, name, data=None, total_size=None, dims=None, transform=UNSET):
         """Register an (un)observed random variable with the model.
@@ -963,7 +1119,10 @@ def register_rv(self, rv_var, name, data=None, total_size=None, dims=None, trans
                 and not isinstance(data, (GenTensorVariable, Minibatch))
                 and data.owner is not None
             ):
-                raise TypeError("Observed data cannot consist of symbolic variables.")
+                raise TypeError(
+                    "Variables that depend on other nodes cannot be used for observed data."
+                    f"The data variable was: {data}"
+                )
 
             data = pandas_to_array(data)
 
@@ -983,6 +1142,7 @@ def make_obs_var(
         ==========
         rv_var
             The random variable that is observed.
+            Its dimensionality must be compatible with the data already.
         data
             The observed data.
         dims: tuple
@@ -994,21 +1154,10 @@ def make_obs_var(
         name = rv_var.name
         data = pandas_to_array(data).astype(rv_var.dtype)
 
-        # The shapes of the observed random variable and its data might not
-        # match.  We need need to update the observed random variable's `size`
-        # (i.e. number of samples) so that it matches the data.
-
-        # Setting `size` produces a random variable with shape `size +
-        # support_shape`, where `len(support_shape) == op.ndim_supp`, we need
-        # to disregard the last `op.ndim_supp`-many dimensions when we
-        # determine the appropriate `size` value from `data.shape`.
-        ndim_supp = rv_var.owner.op.ndim_supp
-        if ndim_supp > 0:
-            new_size = data.shape[:-ndim_supp]
-        else:
-            new_size = data.shape
-
-        rv_var = change_rv_size(rv_var, new_size)
+        if data.ndim != rv_var.ndim:
+            raise ShapeError(
+                "Dimensionality of data and RV don't match.", actual=data.ndim, expected=rv_var.ndim
+            )
 
         if aesara.config.compute_test_value != "off":
             test_value = getattr(rv_var.tag, "test_value", None)
@@ -1132,7 +1281,7 @@ def create_value_var(self, rv_var: TensorVariable, transform: Any) -> TensorVari
 
         return value_var
 
-    def add_random_variable(self, var, dims=None):
+    def add_random_variable(self, var, dims: Optional[Tuple[Union[str, None], ...]] = None):
         """Add a random variable to the named variables of the model."""
         if self.named_vars.tree_contains(var.name):
             raise ValueError(f"Variable name {var.name} already exists.")
@@ -1140,8 +1289,8 @@ def add_random_variable(self, var, dims=None):
         if dims is not None:
             if isinstance(dims, str):
                 dims = (dims,)
-            assert all(dim in self.coords for dim in dims)
-            self.RV_dims[var.name] = dims
+            assert all(dim in self.coords or dim is None for dim in dims)
+            self._RV_dims[var.name] = dims
 
         self.named_vars[var.name] = var
         if not hasattr(self, self.name_of(var.name)):
@@ -1500,18 +1649,7 @@ def set_data(new_data, model=None):
     model = modelcontext(model)
 
     for variable_name, new_value in new_data.items():
-        if isinstance(model[variable_name], SharedVariable):
-            if isinstance(new_value, list):
-                new_value = np.array(new_value)
-            model[variable_name].set_value(pandas_to_array(new_value))
-        else:
-            message = (
-                "The variable `{}` must be defined as `pymc3."
-                "Data` inside the model to allow updating. The "
-                "current type is: "
-                "{}.".format(variable_name, type(model[variable_name]))
-            )
-            raise TypeError(message)
+        model.set_data(variable_name, new_value)
 
 
 def fn(outs, mode=None, model=None, *args, **kwargs):
diff --git a/pymc3/tests/sampler_fixtures.py b/pymc3/tests/sampler_fixtures.py
index 814ed616b7..b4f5cc6cff 100644
--- a/pymc3/tests/sampler_fixtures.py
+++ b/pymc3/tests/sampler_fixtures.py
@@ -92,7 +92,7 @@ class BetaBinomialFixture(KnownCDF):
     @classmethod
     def make_model(cls):
         with pm.Model() as model:
-            p = pm.Beta("p", [0.5, 0.5, 1.0], [0.5, 0.5, 1.0], size=3)
+            p = pm.Beta("p", [0.5, 0.5, 1.0], [0.5, 0.5, 1.0])
             pm.Binomial("y", p=p, n=[4, 12, 9], observed=[1, 2, 9])
         return model
 
@@ -155,13 +155,15 @@ def setup_class(cls):
 
     def test_neff(self):
         if hasattr(self, "min_n_eff"):
-            idata = to_inference_data(self.trace[self.burn :])
+            with self.model:
+                idata = to_inference_data(self.trace[self.burn :])
             n_eff = az.ess(idata)
             for var in n_eff:
                 npt.assert_array_less(self.min_n_eff, n_eff[var])
 
     def test_Rhat(self):
-        idata = to_inference_data(self.trace[self.burn :])
+        with self.model:
+            idata = to_inference_data(self.trace[self.burn :])
         rhat = az.rhat(idata)
         for var in rhat:
             npt.assert_allclose(rhat[var], 1, rtol=0.01)
diff --git a/pymc3/tests/test_data_container.py b/pymc3/tests/test_data_container.py
index 88a1432d48..1ca14088cc 100644
--- a/pymc3/tests/test_data_container.py
+++ b/pymc3/tests/test_data_container.py
@@ -17,10 +17,14 @@
 import pytest
 
 from aesara import shared
+from aesara.tensor.sharedvar import ScalarSharedVariable
+from aesara.tensor.var import TensorVariable
 
 import pymc3 as pm
 
+from pymc3.aesaraf import floatX
 from pymc3.distributions import logpt
+from pymc3.exceptions import ShapeError
 from pymc3.tests.helpers import SeededTest
 
 
@@ -158,22 +162,42 @@ def test_shared_data_as_rv_input(self):
         """
         with pm.Model() as m:
             x = pm.Data("x", [1.0, 2.0, 3.0])
-            _ = pm.Normal("y", mu=x, size=3)
-            trace = pm.sample(
-                chains=1, return_inferencedata=False, compute_convergence_checks=False
+            assert x.eval().shape == (3,)
+            y = pm.Normal("y", mu=x, size=2)
+            assert y.eval().shape == (2, 3)
+            idata = pm.sample(
+                chains=1,
+                tune=500,
+                draws=550,
+                return_inferencedata=True,
+                compute_convergence_checks=False,
             )
+        samples = idata.posterior["y"]
+        assert samples.shape == (1, 550, 2, 3)
 
         np.testing.assert_allclose(np.array([1.0, 2.0, 3.0]), x.get_value(), atol=1e-1)
-        np.testing.assert_allclose(np.array([1.0, 2.0, 3.0]), trace["y"].mean(0), atol=1e-1)
+        np.testing.assert_allclose(
+            np.array([1.0, 2.0, 3.0]), samples.mean(("chain", "draw", "y_dim_0")), atol=1e-1
+        )
 
         with m:
             pm.set_data({"x": np.array([2.0, 4.0, 6.0])})
-            trace = pm.sample(
-                chains=1, return_inferencedata=False, compute_convergence_checks=False
+            assert x.eval().shape == (3,)
+            assert y.eval().shape == (2, 3)
+            idata = pm.sample(
+                chains=1,
+                tune=500,
+                draws=620,
+                return_inferencedata=True,
+                compute_convergence_checks=False,
             )
+        samples = idata.posterior["y"]
+        assert samples.shape == (1, 620, 2, 3)
 
         np.testing.assert_allclose(np.array([2.0, 4.0, 6.0]), x.get_value(), atol=1e-1)
-        np.testing.assert_allclose(np.array([2.0, 4.0, 6.0]), trace["y"].mean(0), atol=1e-1)
+        np.testing.assert_allclose(
+            np.array([2.0, 4.0, 6.0]), samples.mean(("chain", "draw", "y_dim_0")), atol=1e-1
+        )
 
     def test_shared_scalar_as_rv_input(self):
         # See https://github.com/pymc-devs/pymc3/issues/3139
@@ -216,7 +240,7 @@ def test_set_data_to_non_data_container_variables(self):
             )
         with pytest.raises(TypeError) as error:
             pm.set_data({"beta": [1.1, 2.2, 3.3]}, model=model)
-        error.match("defined as `pymc3.Data` inside the model")
+        error.match("The variable `beta` must be a `SharedVariable`")
 
     @pytest.mark.xfail(reason="Depends on ModelGraph")
     def test_model_to_graphviz_for_model_with_data_container(self):
@@ -272,9 +296,47 @@ def test_explicit_coords(self):
 
         assert "rows" in pmodel.coords
         assert pmodel.coords["rows"] == ["R1", "R2", "R3", "R4", "R5"]
+        assert "rows" in pmodel.dim_lengths
+        assert isinstance(pmodel.dim_lengths["rows"], ScalarSharedVariable)
+        assert pmodel.dim_lengths["rows"].eval() == 5
         assert "columns" in pmodel.coords
         assert pmodel.coords["columns"] == ["C1", "C2", "C3", "C4", "C5", "C6", "C7"]
         assert pmodel.RV_dims == {"observations": ("rows", "columns")}
+        assert "columns" in pmodel.dim_lengths
+        assert isinstance(pmodel.dim_lengths["columns"], ScalarSharedVariable)
+        assert pmodel.dim_lengths["columns"].eval() == 7
+
+    def test_symbolic_coords(self):
+        """
+        In v4 dimensions can be created without passing coordinate values.
+        Their lengths are then automatically linked to the corresponding Tensor dimension.
+        """
+        with pm.Model() as pmodel:
+            intensity = pm.Data("intensity", np.ones((2, 3)), dims=("row", "column"))
+            assert "row" in pmodel.dim_lengths
+            assert "column" in pmodel.dim_lengths
+            assert isinstance(pmodel.dim_lengths["row"], TensorVariable)
+            assert isinstance(pmodel.dim_lengths["column"], TensorVariable)
+            assert pmodel.dim_lengths["row"].eval() == 2
+            assert pmodel.dim_lengths["column"].eval() == 3
+
+            intensity.set_value(floatX(np.ones((4, 5))))
+            assert pmodel.dim_lengths["row"].eval() == 4
+            assert pmodel.dim_lengths["column"].eval() == 5
+
+    def test_no_resize_of_implied_dimensions(self):
+        with pm.Model() as pmodel:
+            # Imply a dimension through RV params
+            pm.Normal("n", mu=[1, 2, 3], dims="city")
+            # _Use_ the dimension for a data variable
+            inhabitants = pm.Data("inhabitants", [100, 200, 300], dims="city")
+
+            # Attempting to re-size the dimension through the data variable would
+            # cause shape problems in InferenceData conversion, because the RV remains (3,).
+            with pytest.raises(
+                ShapeError, match="was initialized from 'n' which is not a shared variable"
+            ):
+                pmodel.set_data("inhabitants", [1, 2, 3, 4])
 
     def test_implicit_coords_series(self):
         ser_sales = pd.Series(
diff --git a/pymc3/tests/test_distributions_random.py b/pymc3/tests/test_distributions_random.py
index 8146c132d7..0d13bc5720 100644
--- a/pymc3/tests/test_distributions_random.py
+++ b/pymc3/tests/test_distributions_random.py
@@ -174,12 +174,7 @@ def get_random_variable(self, shape, with_vector_params=False, name=None):
                         # in the test case parametrization "None" means "no specified (default)"
                         return self.distribution(name, transform=None, **params)
                     else:
-                        ndim_supp = self.distribution.rv_op.ndim_supp
-                        if ndim_supp == 0:
-                            size = shape
-                        else:
-                            size = shape[:-ndim_supp]
-                        return self.distribution(name, size=size, transform=None, **params)
+                        return self.distribution(name, shape=shape, transform=None, **params)
                 except TypeError:
                     if np.sum(np.atleast_1d(shape)) == 0:
                         pytest.skip("Timeseries must have positive shape")
@@ -188,10 +183,9 @@ def get_random_variable(self, shape, with_vector_params=False, name=None):
         @staticmethod
         def sample_random_variable(random_variable, size):
             """ Draws samples from a RandomVariable using its .random() method. """
-            if size is None:
-                return random_variable.eval()
-            else:
-                return change_rv_size(random_variable, size, expand=True).eval()
+            if size:
+                random_variable = change_rv_size(random_variable, size, expand=True)
+            return random_variable.eval()
 
         @pytest.mark.parametrize("size", [None, (), 1, (1,), 5, (4, 5)], ids=str)
         @pytest.mark.parametrize("shape", [None, ()], ids=str)
diff --git a/pymc3/tests/test_logp.py b/pymc3/tests/test_logp.py
index aea9db1fdc..f53c640a8f 100644
--- a/pymc3/tests/test_logp.py
+++ b/pymc3/tests/test_logp.py
@@ -70,7 +70,7 @@ def test_logpt_basic():
 
 
 @pytest.mark.parametrize(
-    "indices, size",
+    "indices, shape",
     [
         (slice(0, 2), 5),
         (np.r_[True, True, False, False, True], 5),
@@ -78,15 +78,15 @@ def test_logpt_basic():
         ((np.array([0, 1, 4]), np.array([0, 1, 4])), (5, 5)),
     ],
 )
-def test_logpt_incsubtensor(indices, size):
+def test_logpt_incsubtensor(indices, shape):
     """Make sure we can compute a log-likelihood for ``Y[idx] = data`` where ``Y`` is univariate."""
 
-    mu = floatX(np.power(10, np.arange(np.prod(size)))).reshape(size)
+    mu = floatX(np.power(10, np.arange(np.prod(shape)))).reshape(shape)
     data = mu[indices]
     sigma = 0.001
     rng = aesara.shared(np.random.RandomState(232), borrow=True)
 
-    a = Normal.dist(mu, sigma, size=size, rng=rng)
+    a = Normal.dist(mu, sigma, rng=rng)
     a.name = "a"
 
     a_idx = at.set_subtensor(a[indices], data)
diff --git a/pymc3/tests/test_model.py b/pymc3/tests/test_model.py
index 66cca73169..a155702ba1 100644
--- a/pymc3/tests/test_model.py
+++ b/pymc3/tests/test_model.py
@@ -461,33 +461,29 @@ def test_make_obs_var():
     # Create a fake model and fake distribution to be used for the test
     fake_model = pm.Model()
     with fake_model:
-        fake_distribution = pm.Normal.dist(mu=0, sigma=1)
+        fake_distribution = pm.Normal.dist(mu=0, sigma=1, size=(3, 3))
         # Create the testval attribute simply for the sake of model testing
         fake_distribution.name = input_name
 
     # Check function behavior using the various inputs
+    # dense, sparse: Ensure that the missing values are appropriately set to None
+    # masked: a deterministic variable is returned
+
     dense_output = fake_model.make_obs_var(fake_distribution, dense_input, None, None)
+    assert dense_output == fake_distribution
+    assert isinstance(dense_output.tag.observations, TensorConstant)
     del fake_model.named_vars[fake_distribution.name]
+
     sparse_output = fake_model.make_obs_var(fake_distribution, sparse_input, None, None)
+    assert sparse_output == fake_distribution
+    assert sparse.basic._is_sparse_variable(sparse_output.tag.observations)
     del fake_model.named_vars[fake_distribution.name]
+
+    # Here the RandomVariable is split into observed/imputed and a Deterministic is returned
     masked_output = fake_model.make_obs_var(fake_distribution, masked_array_input, None, None)
+    assert masked_output != fake_distribution
     assert not isinstance(masked_output, RandomVariable)
-
-    # Ensure that the missing values are appropriately set to None
-    for func_output in [dense_output, sparse_output]:
-        assert isinstance(func_output.owner.op, RandomVariable)
-
-    # Ensure that the Aesara variable names are correctly set.
-    # Note that the output for masked inputs do not have their names set
-    # to the passed value.
-    for func_output in [dense_output, sparse_output]:
-        assert func_output.name == input_name
-
-    # Ensure the that returned functions are all of the correct type
-    assert isinstance(dense_output.tag.observations, TensorConstant)
-    assert sparse.basic._is_sparse_variable(sparse_output.tag.observations)
-
-    # Masked output is something weird. Just ensure it has missing values
+    # Ensure it has missing values
     assert {"testing_inputs_missing"} == {v.name for v in fake_model.vars}
     assert {"testing_inputs", "testing_inputs_observed"} == {
         v.name for v in fake_model.observed_RVs
diff --git a/pymc3/tests/test_ode.py b/pymc3/tests/test_ode.py
index 94dfb0dd6f..d2e2e25b90 100644
--- a/pymc3/tests/test_ode.py
+++ b/pymc3/tests/test_ode.py
@@ -168,6 +168,7 @@ def ode_func_5(y, t, p):
         np.testing.assert_array_equal(np.ravel(model5_sens_ic), model5._sens_ic)
 
 
+@pytest.mark.xfail(reason="See https://github.com/pymc-devs/pymc3/issues/4652.")
 def test_logp_scalar_ode():
     """Test the computation of the log probability for these models"""
 
diff --git a/pymc3/tests/test_sampling.py b/pymc3/tests/test_sampling.py
index 71b8a9b8dc..217cf96010 100644
--- a/pymc3/tests/test_sampling.py
+++ b/pymc3/tests/test_sampling.py
@@ -213,7 +213,7 @@ def test_return_inferencedata(self, monkeypatch):
                 return_inferencedata=True,
                 discard_tuned_samples=True,
                 idata_kwargs={"prior": prior},
-                random_seed=-1
+                random_seed=-1,
             )
             assert "prior" in result
             assert isinstance(result, InferenceData)
@@ -385,11 +385,10 @@ def test_shared_named(self):
                 "theta0",
                 mu=np.atleast_2d(0),
                 tau=np.atleast_2d(1e20),
-                size=(1, 1),
                 testval=np.atleast_2d(0),
             )
             theta = pm.Normal(
-                "theta", mu=at.dot(G_var, theta0), tau=np.atleast_2d(1e20), size=(1, 1)
+                "theta", mu=at.dot(G_var, theta0), tau=np.atleast_2d(1e20), shape=(1, 1)
             )
             res = theta.eval()
             assert np.isclose(res, 0.0)
@@ -401,11 +400,10 @@ def test_shared_unnamed(self):
                 "theta0",
                 mu=np.atleast_2d(0),
                 tau=np.atleast_2d(1e20),
-                size=(1, 1),
                 testval=np.atleast_2d(0),
             )
             theta = pm.Normal(
-                "theta", mu=at.dot(G_var, theta0), tau=np.atleast_2d(1e20), size=(1, 1)
+                "theta", mu=at.dot(G_var, theta0), tau=np.atleast_2d(1e20), shape=(1, 1)
             )
             res = theta.eval()
             assert np.isclose(res, 0.0)
@@ -417,11 +415,10 @@ def test_constant_named(self):
                 "theta0",
                 mu=np.atleast_2d(0),
                 tau=np.atleast_2d(1e20),
-                size=(1, 1),
                 testval=np.atleast_2d(0),
             )
             theta = pm.Normal(
-                "theta", mu=at.dot(G_var, theta0), tau=np.atleast_2d(1e20), size=(1, 1)
+                "theta", mu=at.dot(G_var, theta0), tau=np.atleast_2d(1e20), shape=(1, 1)
             )
 
             res = theta.eval()
@@ -936,14 +933,13 @@ def test_ignores_observed(self):
         npt.assert_array_almost_equal(prior["positive_mu"], np.abs(prior["mu"]), decimal=4)
 
     def test_respects_shape(self):
-        for shape in (2, (2,), (10, 2), (10, 10)):
+        for shape in ((2,), (10, 2), (10, 10)):
             with pm.Model():
-                mu = pm.Gamma("mu", 3, 1, size=1)
-                goals = pm.Poisson("goals", mu, size=shape)
+                mu = pm.Gamma("mu", 3, 1)
+                goals = pm.Poisson("goals", mu, shape=shape)
+                assert goals.eval().shape == shape, f"Current shape setting: {shape}"
                 trace1 = pm.sample_prior_predictive(10, var_names=["mu", "mu", "goals"])
                 trace2 = pm.sample_prior_predictive(10, var_names=["mu", "goals"])
-            if shape == 2:  # want to test shape as an int
-                shape = (2,)
             assert trace1["goals"].shape == (10,) + shape
             assert trace2["goals"].shape == (10,) + shape
 
@@ -971,14 +967,14 @@ def test_multivariate2(self):
 
     def test_layers(self):
         with pm.Model() as model:
-            a = pm.Uniform("a", lower=0, upper=1, size=10)
-            b = pm.Binomial("b", n=1, p=a, size=10)
+            a = pm.Uniform("a", lower=0, upper=1, size=5)
+            b = pm.Binomial("b", n=1, p=a, size=7)
 
         model.default_rng.get_value(borrow=True).seed(232093)
 
         b_sampler = aesara.function([], b)
         avg = np.stack([b_sampler() for i in range(10000)]).mean(0)
-        npt.assert_array_almost_equal(avg, 0.5 * np.ones((10,)), decimal=2)
+        npt.assert_array_almost_equal(avg, 0.5 * np.ones((7, 5)), decimal=2)
 
     def test_transformed(self):
         n = 18
@@ -1087,16 +1083,15 @@ def test_sample_from_xarray_prior(self, point_list_arg_bug_fixture):
 
         with pmodel:
             prior = pm.sample_prior_predictive(samples=20)
-
-        idat = pm.to_inference_data(trace, prior=prior)
+            idat = pm.to_inference_data(trace, prior=prior)
 
         with pmodel:
             pp = pm.sample_posterior_predictive(idat.prior, var_names=["d"])
 
     def test_sample_from_xarray_posterior(self, point_list_arg_bug_fixture):
         pmodel, trace = point_list_arg_bug_fixture
-        idat = pm.to_inference_data(trace)
         with pmodel:
+            idat = pm.to_inference_data(trace)
             pp = pm.sample_posterior_predictive(idat.posterior, var_names=["d"])
 
 
diff --git a/pymc3/tests/test_shape_handling.py b/pymc3/tests/test_shape_handling.py
index 37c0619322..3c93721f15 100644
--- a/pymc3/tests/test_shape_handling.py
+++ b/pymc3/tests/test_shape_handling.py
@@ -11,7 +11,7 @@
 #   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 #   See the License for the specific language governing permissions and
 #   limitations under the License.
-
+import aesara
 import numpy as np
 import pytest
 
@@ -19,6 +19,7 @@
 
 import pymc3 as pm
 
+from pymc3.distributions.distribution import _validate_shape_dims_size
 from pymc3.distributions.shape_utils import (
     broadcast_dist_samples_shape,
     broadcast_dist_samples_to,
@@ -219,3 +220,197 @@ def test_sample_generate_values(fixture_model, fixture_sizes):
         prior = pm.sample_prior_predictive(samples=fixture_sizes)
         for rv in RVs:
             assert prior[rv.name].shape == size + tuple(rv.distribution.shape)
+
+
+class TestShapeDimsSize:
+    @pytest.mark.parametrize("param_shape", [(), (3,)])
+    @pytest.mark.parametrize("batch_shape", [(), (3,)])
+    @pytest.mark.parametrize(
+        "parametrization",
+        [
+            "implicit",
+            "shape",
+            "shape...",
+            "dims",
+            "dims...",
+            "size",
+        ],
+    )
+    def test_param_and_batch_shape_combos(
+        self, param_shape: tuple, batch_shape: tuple, parametrization: str
+    ):
+        coords = {}
+        param_dims = []
+        batch_dims = []
+
+        # Create coordinates corresponding to the parameter shape
+        for d in param_shape:
+            dname = f"param_dim_{d}"
+            coords[dname] = [f"c_{i}" for i in range(d)]
+            param_dims.append(dname)
+        assert len(param_dims) == len(param_shape)
+        # Create coordinates corresponding to the batch shape
+        for d in batch_shape:
+            dname = f"batch_dim_{d}"
+            coords[dname] = [f"c_{i}" for i in range(d)]
+            batch_dims.append(dname)
+        assert len(batch_dims) == len(batch_shape)
+
+        with pm.Model(coords=coords) as pmodel:
+            mu = aesara.shared(np.random.normal(size=param_shape))
+
+            with pytest.warns(None):
+                if parametrization == "implicit":
+                    rv = pm.Normal("rv", mu=mu).shape == param_shape
+                else:
+                    if parametrization == "shape":
+                        rv = pm.Normal("rv", mu=mu, shape=batch_shape + param_shape)
+                        assert rv.eval().shape == batch_shape + param_shape
+                    elif parametrization == "shape...":
+                        rv = pm.Normal("rv", mu=mu, shape=(*batch_shape, ...))
+                        assert rv.eval().shape == batch_shape + param_shape
+                    elif parametrization == "dims":
+                        rv = pm.Normal("rv", mu=mu, dims=batch_dims + param_dims)
+                        assert rv.eval().shape == batch_shape + param_shape
+                    elif parametrization == "dims...":
+                        rv = pm.Normal("rv", mu=mu, dims=(*batch_dims, ...))
+                        n_size = len(batch_shape)
+                        n_implied = len(param_shape)
+                        ndim = n_size + n_implied
+                        assert len(pmodel.RV_dims["rv"]) == ndim, pmodel.RV_dims
+                        assert len(pmodel.RV_dims["rv"][:n_size]) == len(batch_dims)
+                        assert len(pmodel.RV_dims["rv"][n_size:]) == len(param_dims)
+                        if n_implied > 0:
+                            assert pmodel.RV_dims["rv"][-1] is None
+                    elif parametrization == "size":
+                        rv = pm.Normal("rv", mu=mu, size=batch_shape)
+                        assert rv.eval().shape == batch_shape + param_shape
+                    else:
+                        raise NotImplementedError("Invalid test case parametrization.")
+
+    def test_define_dims_on_the_fly(self):
+        with pm.Model() as pmodel:
+            agedata = aesara.shared(np.array([10, 20, 30]))
+
+            # Associate the "patient" dim with an implied dimension
+            age = pm.Normal("age", agedata, dims=("patient",))
+            assert "patient" in pmodel.dim_lengths
+            assert pmodel.dim_lengths["patient"].eval() == 3
+
+            # Use the dim to replicate a new RV
+            effect = pm.Normal("effect", 0, dims=("patient",))
+            assert effect.ndim == 1
+            assert effect.eval().shape == (3,)
+
+            # Now change the length of the implied dimension
+            agedata.set_value([1, 2, 3, 4])
+            # The change should propagate all the way through
+            assert effect.eval().shape == (4,)
+
+    @pytest.mark.xfail(reason="Simultaneous use of size and dims is not implemented")
+    def test_data_defined_size_dimension_can_register_dimname(self):
+        with pm.Model() as pmodel:
+            x = pm.Data("x", [[1, 2, 3, 4]], dims=("first", "second"))
+            assert "first" in pmodel.dim_lengths
+            assert "second" in pmodel.dim_lengths
+            # two dimensions are implied; a "third" dimension is created
+            y = pm.Normal("y", mu=x, size=2, dims=("third", "first", "second"))
+            assert "third" in pmodel.dim_lengths
+            assert y.eval().shape() == (2, 1, 4)
+
+    def test_can_resize_data_defined_size(self):
+        with pm.Model() as pmodel:
+            x = pm.Data("x", [[1, 2, 3, 4]], dims=("first", "second"))
+            y = pm.Normal("y", mu=0, dims=("first", "second"))
+            z = pm.Normal("z", mu=y, observed=np.ones((1, 4)))
+            assert x.eval().shape == (1, 4)
+            assert y.eval().shape == (1, 4)
+            assert z.eval().shape == (1, 4)
+            assert "first" in pmodel.dim_lengths
+            assert "second" in pmodel.dim_lengths
+            pmodel.set_data("x", [[1, 2], [3, 4], [5, 6]])
+            assert x.eval().shape == (3, 2)
+            assert y.eval().shape == (3, 2)
+            assert z.eval().shape == (3, 2)
+
+    @pytest.mark.xfail(reason="See https://github.com/pymc-devs/pymc3/issues/4652.")
+    def test_observed_with_column_vector(self):
+        with pm.Model() as model:
+            pm.Normal("x1", mu=0, sd=1, observed=np.random.normal(size=(3, 4)))
+            model.logp()
+            pm.Normal("x2", mu=0, sd=1, observed=np.random.normal(size=(3, 1)))
+            model.logp()
+
+    def test_dist_api_works(self):
+        mu = aesara.shared(np.array([1, 2, 3]))
+        with pytest.raises(NotImplementedError, match="API is not yet supported"):
+            pm.Normal.dist(mu=mu, dims=("town",))
+        assert pm.Normal.dist(mu=mu, shape=(3,)).eval().shape == (3,)
+        assert pm.Normal.dist(mu=mu, shape=(5, 3)).eval().shape == (5, 3)
+        assert pm.Normal.dist(mu=mu, shape=(7, ...)).eval().shape == (7, 3)
+        assert pm.Normal.dist(mu=mu, size=(4,)).eval().shape == (4, 3)
+
+    def test_auto_assert_shape(self):
+        with pytest.raises(AssertionError, match="will never match"):
+            pm.Normal.dist(mu=[1, 2], shape=[])
+
+        mu = at.vector(name="mu_input")
+        rv = pm.Normal.dist(mu=mu, shape=[3, 4])
+        f = aesara.function([mu], rv, mode=aesara.Mode("py"))
+        assert f([1, 2, 3, 4]).shape == (3, 4)
+
+        with pytest.raises(AssertionError, match=r"Got shape \(3, 2\), expected \(3, 4\)."):
+            f([1, 2])
+
+        # The `shape` can be symbolic!
+        s = at.vector(dtype="int32")
+        rv = pm.Uniform.dist(2, [4, 5], shape=s)
+        f = aesara.function([s], rv, mode=aesara.Mode("py"))
+        f(
+            [
+                2,
+            ]
+        )
+        with pytest.raises(
+            AssertionError,
+            match=r"Got 1 dimensions \(shape \(2,\)\), expected 2 dimensions with shape \(3, 4\).",
+        ):
+            f([3, 4])
+        with pytest.raises(
+            AssertionError,
+            match=r"Got 1 dimensions \(shape \(2,\)\), expected 0 dimensions with shape \(\).",
+        ):
+            f([])
+        pass
+
+    def test_lazy_flavors(self):
+
+        _validate_shape_dims_size(shape=5)
+        _validate_shape_dims_size(dims="town")
+        _validate_shape_dims_size(size=7)
+
+        assert pm.Uniform.dist(2, [4, 5], size=[3, 4]).eval().shape == (3, 4, 2)
+        assert pm.Uniform.dist(2, [4, 5], shape=[3, 2]).eval().shape == (3, 2)
+        with pm.Model(coords=dict(town=["Greifswald", "Madrid"])):
+            assert pm.Normal("n2", mu=[1, 2], dims=("town",)).eval().shape == (2,)
+
+    def test_invalid_flavors(self):
+        # redundant parametrizations
+        with pytest.raises(ValueError, match="Passing both"):
+            _validate_shape_dims_size(shape=(2,), dims=("town",))
+        with pytest.raises(ValueError, match="Passing both"):
+            _validate_shape_dims_size(dims=("town",), size=(2,))
+        with pytest.raises(ValueError, match="Passing both"):
+            _validate_shape_dims_size(shape=(3,), size=(3,))
+
+        # invalid, but not necessarly rare
+        with pytest.raises(ValueError, match="must be an int, list or tuple"):
+            _validate_shape_dims_size(size="notasize")
+
+        # invalid ellipsis positions
+        with pytest.raises(ValueError, match="may only appear in the last position"):
+            _validate_shape_dims_size(shape=(3, ..., 2))
+        with pytest.raises(ValueError, match="may only appear in the last position"):
+            _validate_shape_dims_size(dims=(..., "town"))
+        with pytest.raises(ValueError, match="cannot contain"):
+            _validate_shape_dims_size(size=(3, ...))
diff --git a/pymc3/tests/test_step.py b/pymc3/tests/test_step.py
index fd02139879..6ef93eb5d1 100644
--- a/pymc3/tests/test_step.py
+++ b/pymc3/tests/test_step.py
@@ -982,7 +982,7 @@ def test_linalg(self, caplog):
             a = Normal("a", size=2, testval=floatX(np.zeros(2)))
             a = at.switch(a > 0, np.inf, a)
             b = at.slinalg.solve(floatX(np.eye(2)), a)
-            Normal("c", mu=b, size=2, testval=floatX(np.r_[0.0, 0.0]))
+            Normal("c", mu=b, shape=(2,), testval=floatX(np.r_[0.0, 0.0]))
             caplog.clear()
             trace = sample(20, init=None, tune=5, chains=2)
             warns = [msg.msg for msg in caplog.records]
diff --git a/pymc3/tests/test_transforms.py b/pymc3/tests/test_transforms.py
index fd32d8b9b6..e040e6e244 100644
--- a/pymc3/tests/test_transforms.py
+++ b/pymc3/tests/test_transforms.py
@@ -274,11 +274,11 @@ def test_chain_jacob_det():
 
 
 class TestElementWiseLogp(SeededTest):
-    def build_model(self, distfam, params, size, transform, testval=None):
+    def build_model(self, distfam, params, *, size=None, shape=None, transform=None, testval=None):
         if testval is not None:
             testval = pm.floatX(testval)
         with pm.Model() as m:
-            distfam("x", size=size, transform=transform, testval=testval, **params)
+            distfam("x", size=size, shape=shape, transform=transform, testval=testval, **params)
         return m
 
     def check_transform_elementwise_logp(self, model):
@@ -328,32 +328,34 @@ def test_half_normal(self, sd, size):
         model = self.build_model(pm.HalfNormal, {"sd": sd}, size=size, transform=tr.log)
         self.check_transform_elementwise_logp(model)
 
-    @pytest.mark.parametrize("lam,size", [(2.5, 2), (5.0, (2, 3)), (np.ones(3), (4, 3))])
+    @pytest.mark.parametrize("lam,size", [(2.5, 2), (5.0, (2, 3)), (np.ones(3), (4, 5))])
     def test_exponential(self, lam, size):
         model = self.build_model(pm.Exponential, {"lam": lam}, size=size, transform=tr.log)
         self.check_transform_elementwise_logp(model)
 
     @pytest.mark.parametrize(
-        "a,b,size",
+        "a,b,shape",
         [
             (1.0, 1.0, 2),
             (0.5, 0.5, (2, 3)),
             (np.ones(3), np.ones(3), (4, 3)),
         ],
     )
-    def test_beta(self, a, b, size):
-        model = self.build_model(pm.Beta, {"alpha": a, "beta": b}, size=size, transform=tr.logodds)
+    def test_beta(self, a, b, shape):
+        model = self.build_model(
+            pm.Beta, {"alpha": a, "beta": b}, shape=shape, transform=tr.logodds
+        )
         self.check_transform_elementwise_logp(model)
 
     @pytest.mark.parametrize(
-        "lower,upper,size",
+        "lower,upper,shape",
         [
             (0.0, 1.0, 2),
             (0.5, 5.5, (2, 3)),
             (pm.floatX(np.zeros(3)), pm.floatX(np.ones(3)), (4, 3)),
         ],
     )
-    def test_uniform(self, lower, upper, size):
+    def test_uniform(self, lower, upper, shape):
         def transform_params(rv_var):
             _, _, _, lower, upper = rv_var.owner.inputs
             lower = at.as_tensor_variable(lower) if lower is not None else None
@@ -362,25 +364,25 @@ def transform_params(rv_var):
 
         interval = tr.Interval(transform_params)
         model = self.build_model(
-            pm.Uniform, {"lower": lower, "upper": upper}, size=size, transform=interval
+            pm.Uniform, {"lower": lower, "upper": upper}, shape=shape, transform=interval
         )
         self.check_transform_elementwise_logp(model)
 
     @pytest.mark.parametrize(
-        "mu,kappa,size", [(0.0, 1.0, 2), (-0.5, 5.5, (2, 3)), (np.zeros(3), np.ones(3), (4, 3))]
+        "mu,kappa,shape", [(0.0, 1.0, 2), (-0.5, 5.5, (2, 3)), (np.zeros(3), np.ones(3), (4, 3))]
     )
     @pytest.mark.xfail(reason="Distribution not refactored yet")
-    def test_vonmises(self, mu, kappa, size):
+    def test_vonmises(self, mu, kappa, shape):
         model = self.build_model(
-            pm.VonMises, {"mu": mu, "kappa": kappa}, size=size, transform=tr.circular
+            pm.VonMises, {"mu": mu, "kappa": kappa}, shape=shape, transform=tr.circular
         )
         self.check_transform_elementwise_logp(model)
 
     @pytest.mark.parametrize(
-        "a,size", [(np.ones(2), None), (np.ones((2, 3)) * 0.5, None), (np.ones(3), (4,))]
+        "a,shape", [(np.ones(2), None), (np.ones((2, 3)) * 0.5, None), (np.ones(3), (4,))]
     )
-    def test_dirichlet(self, a, size):
-        model = self.build_model(pm.Dirichlet, {"a": a}, size=size, transform=tr.stick_breaking)
+    def test_dirichlet(self, a, shape):
+        model = self.build_model(pm.Dirichlet, {"a": a}, shape=shape, transform=tr.stick_breaking)
         self.check_vectortransform_elementwise_logp(model, vect_opt=1)
 
     def test_normal_ordered(self):
@@ -394,59 +396,59 @@ def test_normal_ordered(self):
         self.check_vectortransform_elementwise_logp(model, vect_opt=0)
 
     @pytest.mark.parametrize(
-        "sd,size",
+        "sd,shape",
         [
             (2.5, (2,)),
             (np.ones(3), (4, 3)),
         ],
     )
     @pytest.mark.xfail(condition=(aesara.config.floatX == "float32"), reason="Fails on float32")
-    def test_half_normal_ordered(self, sd, size):
-        testval = np.sort(np.abs(np.random.randn(*size)))
+    def test_half_normal_ordered(self, sd, shape):
+        testval = np.sort(np.abs(np.random.randn(*shape)))
         model = self.build_model(
             pm.HalfNormal,
             {"sd": sd},
-            size=size,
+            shape=shape,
             testval=testval,
             transform=tr.Chain([tr.log, tr.ordered]),
         )
         self.check_vectortransform_elementwise_logp(model, vect_opt=0)
 
-    @pytest.mark.parametrize("lam,size", [(2.5, (2,)), (np.ones(3), (4, 3))])
-    def test_exponential_ordered(self, lam, size):
-        testval = np.sort(np.abs(np.random.randn(*size)))
+    @pytest.mark.parametrize("lam,shape", [(2.5, (2,)), (np.ones(3), (4, 3))])
+    def test_exponential_ordered(self, lam, shape):
+        testval = np.sort(np.abs(np.random.randn(*shape)))
         model = self.build_model(
             pm.Exponential,
             {"lam": lam},
-            size=size,
+            shape=shape,
             testval=testval,
             transform=tr.Chain([tr.log, tr.ordered]),
         )
         self.check_vectortransform_elementwise_logp(model, vect_opt=0)
 
     @pytest.mark.parametrize(
-        "a,b,size",
+        "a,b,shape",
         [
             (1.0, 1.0, (2,)),
             (np.ones(3), np.ones(3), (4, 3)),
         ],
     )
-    def test_beta_ordered(self, a, b, size):
-        testval = np.sort(np.abs(np.random.rand(*size)))
+    def test_beta_ordered(self, a, b, shape):
+        testval = np.sort(np.abs(np.random.rand(*shape)))
         model = self.build_model(
             pm.Beta,
             {"alpha": a, "beta": b},
-            size=size,
+            shape=shape,
             testval=testval,
             transform=tr.Chain([tr.logodds, tr.ordered]),
         )
         self.check_vectortransform_elementwise_logp(model, vect_opt=0)
 
     @pytest.mark.parametrize(
-        "lower,upper,size",
+        "lower,upper,shape",
         [(0.0, 1.0, (2,)), (pm.floatX(np.zeros(3)), pm.floatX(np.ones(3)), (4, 3))],
     )
-    def test_uniform_ordered(self, lower, upper, size):
+    def test_uniform_ordered(self, lower, upper, shape):
         def transform_params(rv_var):
             _, _, _, lower, upper = rv_var.owner.inputs
             lower = at.as_tensor_variable(lower) if lower is not None else None
@@ -455,43 +457,45 @@ def transform_params(rv_var):
 
         interval = tr.Interval(transform_params)
 
-        testval = np.sort(np.abs(np.random.rand(*size)))
+        testval = np.sort(np.abs(np.random.rand(*shape)))
         model = self.build_model(
             pm.Uniform,
             {"lower": lower, "upper": upper},
-            size=size,
+            shape=shape,
             testval=testval,
             transform=tr.Chain([interval, tr.ordered]),
         )
         self.check_vectortransform_elementwise_logp(model, vect_opt=1)
 
-    @pytest.mark.parametrize("mu,kappa,size", [(0.0, 1.0, (2,)), (np.zeros(3), np.ones(3), (4, 3))])
+    @pytest.mark.parametrize(
+        "mu,kappa,shape", [(0.0, 1.0, (2,)), (np.zeros(3), np.ones(3), (4, 3))]
+    )
     @pytest.mark.xfail(reason="Distribution not refactored yet")
-    def test_vonmises_ordered(self, mu, kappa, size):
-        testval = np.sort(np.abs(np.random.rand(*size)))
+    def test_vonmises_ordered(self, mu, kappa, shape):
+        testval = np.sort(np.abs(np.random.rand(*shape)))
         model = self.build_model(
             pm.VonMises,
             {"mu": mu, "kappa": kappa},
-            size=size,
+            shape=shape,
             testval=testval,
             transform=tr.Chain([tr.circular, tr.ordered]),
         )
         self.check_vectortransform_elementwise_logp(model, vect_opt=0)
 
     @pytest.mark.parametrize(
-        "lower,upper,size,transform",
+        "lower,upper,shape,transform",
         [
             (0.0, 1.0, (2,), tr.stick_breaking),
             (0.5, 5.5, (2, 3), tr.stick_breaking),
             (np.zeros(3), np.ones(3), (4, 3), tr.Chain([tr.sum_to_1, tr.logodds])),
         ],
     )
-    def test_uniform_other(self, lower, upper, size, transform):
-        testval = np.ones(size) / size[-1]
+    def test_uniform_other(self, lower, upper, shape, transform):
+        testval = np.ones(shape) / shape[-1]
         model = self.build_model(
             pm.Uniform,
             {"lower": lower, "upper": upper},
-            size=size,
+            shape=shape,
             testval=testval,
             transform=transform,
         )