Discussion: Dims usage on univariate regression is a bit confusing

[canyon289]

Description of your problem

Using coords or dims in univariate situations necessitates the usage of indexing. This can be confusing as if I have a model that is specified without indexing, and then add dims/coords, an shape mismatch is caused. If the behavior can't be changed can a more useful exception be raised?

Please provide a minimal, self-contained, and reproducible example.

# Generate data and specify coords
df = pd.DataFrame(np.random.normal(size=(10,2)), columns=["x_1", "y"])

COORDS = {
    "slopes": ["x_1"]
}

Doesnt work

with pm.Model(coords=COORDS) as adelie_regression:

    β = pm.Normal("β", sigma=1000, dims=("slopes"))
    
    mu =  β * df["x_1"].values
    pm.Normal(
        "obs",
        mu=mu,
        sigma=1,
        observed=df["y"].values
    )
    idata_regression = pm.sample()

Works

with pm.Model(coords=COORDS) as adelie_regression:

    β = pm.Normal("β", sigma=1000, dims=("slopes"))
    
    mu =  β[0] * df["x_1"].values
    pm.Normal(
        "obs",
        mu=mu,
        sigma=1,
        observed=df["y"].values
    )
    idata_regression = pm.sample()

Also works

with pm.Model() as adelie_regression:

    β = pm.Normal("β", sigma=1000)
    
    mu =  β* df["x_1"].values
    pm.Normal(
        "obs",
        mu=mu,
        sigma=1,
        observed=df["y"].values
    )
    idata_regression = pm.sample()

Exception for first case

Versions and main components

• PyMC/PyMC3 Version: Main at 620546b
• Aesara/Theano Version:

[OriolAbril]

I am not sure there is any sensible way to get around this. dims gets converted to size so what you are asking boils down to having the two options below be equivalent:

β = pm.Normal("β", sigma=1000)
β = pm.Normal("β", sigma=1000, size=(1,))

Which I think will probably mean many more shape issues down the road.

[ricardoV94]

I am not sure there is any sensible way to get around this. dims gets converted to size so what you are asking boils down to having the two options below be equivalent:

β = pm.Normal("β", sigma=1000)
β = pm.Normal("β", sigma=1000, size=(1,))

Which I think will probably mean many more shape issues down the road.

Or alternatively allow dims for scalar distributions (size=())

[OriolAbril]

Or alternatively allow dims for scalar distributions (size=())

Not sure I understand this, could you explain how that would be defined and how the output would look like?

[ricardoV94]

Or alternatively allow dims for scalar distributions (size=())

Not sure I understand this, could you explain how that would be defined and how the output would look like?

I don't have a good suggestion, just wanted to mention the flip side of what you said. We have a way to specify dims for size >= 1, but not for size=() which is what the original problem hinged on.

Not saying it is incoherent, since they have zero dimensions according to the numpy API, therefore dims should be None. We also can't say size=(0,)

I hope this clarifies the original question.

[OriolAbril]

@canyon289 could you explain why do you want dimension names on scalar quantities at a conceptual level? How do they differ from the variable names?

I also ignored the changing the error message previously. I think the most important thing would be to be able to distinguish such cases from "regular" shape mismatches, however, I am not yet sure I understand what "such cases" means, much less if or how they could be distinguished

[canyon289]

To be clear I'm not suggesting that dimension names with scalar quantities has to work, or that I want them implemented. All I'm saying is that the exception raised is confusing because it "feels" like it should work, especially in a workflow where multiple dimensions are being added iteratively.
,
Perhaps we can just raise an exception earlier if we detect it saying that we don't support that behavior so users don't get a dimension mismatch later are confused what's happening. Just a suggestion!

[ricardoV94]

@canyon289 how would you distinguish dim for scalar size and size=1? Also I don't know whether InferenceData allows this either?

[ricardoV94]

Conceptually you are asking to give a name to a dim that does not exist.

assert pm.Normal.dist().ndim == 0

[twiecki]

The current behavior is what I would have expected and is internally consistent. If I specify a dims I get a vector. I'd be more confused if I wrote code that worked on vectors and then change the dims to 1 and it stops working.

The shape errors are indeed not the best, which is a more general problem. Would be good if they mentioned the name of RV at least.

[OriolAbril]

To be clear I'm suggesting that dimension names with scalar quantities has to work, or that I want them implementing.

I read this as

To be clear I'm not suggesting that dimension names with scalar quantities has to work, or that I want them implementing.

Otherwise the rest of the paragraph doesn't make mich sense.

Perhaps we can just raise an exception earlier if we detect it saying that we don't support that behavior so users don't get a dimension mismatch later are confused what's happening.

Here I don't understand what you refer to with "that behaviour".

The shape errors are indeed not the best, which is a more general problem. Would be good if they mentioned the name of RV at least.

That would be nice

---

Also, looking back at the first code snippet, the column in the dataframe is selected with a string, not with a list, so it is a 1d vector of length 10, not a 1x10 2d array. Shouldn't two 1d vectors of lengths 1 and 10 broadcast?

[canyon289]

@OriolAbril my apologies i fixed the original post.
For the behavior its raising a more useful message to the user if dimension names on scalar quantities is not supported.

@canyon289 how would you distinguish dim for scalar size and size=1? Also I don't know whether InferenceData allows this either?

I have no idea, in this situation I'm a user that's confused, not a developer that has answers :D. I don't mean that to be flippant, just want to be clear where my perspective is coming from.

I'm getting the understanding that at an "automagic fix" wouldnt be the best route! That's alright!
Would it be possible to "detect" the use of dimension names on scalar quantities on scalar quantities to tell users this may cause issues shape issues later, either through an exception or logging?

[OriolAbril]

Would it be possible to "detect" the use of dimension names on scalar quantities on scalar quantities to tell users this may cause issues shape issues later, either through an exception or logging?

This issue is proof that this can't be done. You expected the code from the first snippet to be a scalar but I have used that same code several times because I want a 1d array of length 1. I think the best way is probably adding examples of that on the dims/size notebook so that it is clear to users that these two are equivalent:

with pm.Model(coords={"slopes": ["x_1"]}):
    pm.Normal("b", 0, 1, size=(1,))
    pm.Normal("a", 0, 1, dims="slopes")

and to get a scalar neither size nor dims should be used.

[ricardoV94]

Would this be wrong?

with pm.Model(coords={"slopes": ()}):
    a = pm.Normal("a", 0, 1, size=())
    # This is treated as size=(0,), and not size=(), or size=None right now
    b = pm.Normal("b", 0, 1, dims="slopes") 

a.eval(), b.eval()
# array(1.2774636), array([], dtype=float64)

[OriolAbril]

Are size=5 and size=(5,) equivalent?