Make two random variables equal

Questions v5
1

Hi!

I have a quite hard question that have haunting me for months.

I’m building a complicated model which looks like this:

image
image1116×1004 96.2 KB

In the red circle, I make the random variables “ARR” and “LARR” equal by soft zero difference:

    zerodiff = pm.Normal.dist(mu=0, sigma=0.01)
    pm.Potential(f"zerodiff_potential", pm.logp(zerodiff, ARR - LARR))

My questions are:

  1. ARR and LARR should be exactly the same. But it’s hard to handle the hierarchy, so I draw them with two different pathways and make them equal. Is this the best way to deal with these problems?
  2. The soft difference seems to make the chain harder to converge, any suggestions?

Thanks a lot!

Yangkang

2

LARR doesn’t have any parent nodes (except for the hierarchical parameters), so why not merge the two variables and have T_{\mu_{LARRTG}} and T_{\mu_{LARR}} inherit from \mu_{ARR}?

3

Thanks for the timely reply! Yeah I haven’t thought about it… There are some difficulties because the 𝜇_𝐴𝑅𝑅 and 𝜇_𝐿𝐴𝑅𝑅 are of different shape, but I will try on that way. There is another “zero difference” in the model:

image
image1810×990 151 KB

Where the Tmin and TM are made equal. In this case the TM is modeled with spatial autoregression parameters. Any suggestions on this?

Thanks!

4

If \mu_{ARR} and \mu_{LARR} are not the same shape, how can they be equal?

5

The mu_LARR is the mean of LARR across years, and LARR are drawn from it using year indexing. On the other hand, mu_ARR is the full size.

6

Can you use a deterministic for LARR then? Keep ARR as-is, and use LARR = pm.Deterministic(ARR[year_indices])?

7

Yeah I think so. Also could be more complex.

    ### mu ARR
    mu_ARR = pm.Deterministic('mu_ARR', beta_IAR[year_index_obs] + gamma_IAR[cell_index_obs] + \
                                          phi_IAR[cell_index_obs] * sigma_phi_IAR)
    sigma_ARR = pm.HalfNormal('sigma_ARR', sigma=100)
    
    ### ARR
    ARR = pm.Normal('ARR', mu=mu_ARR, sigma=sigma_ARR) ## Actually ARRij

The ARR is a combination of year effect and cell effect (spatial effect). I think I can remove the beta_IAR to get a mean_ARR across year, and then expand it into the full size by adding the year effect.

1 Like
8

@jessegrabowski I think I’m facing another problem in the model. Basically like “How to do regression with two RVs”. If Y and X are both RVs to be estimated, is there a way to also estimate the slope and intercept parameters at the same time?

9

Yeah, it just means that the value of Y will depend on X, alpha, and beta. Nothing stops you from writing a “regression prior” like this:

    X = pm.Normal('X', 0, 1, size=(100,3 ))
    alpha = pm.Normal('alpha', 0, 1)
    beta = pm.Normal('beta', 0, 1, size=3)
    y = pm.Normal('y', alpha + X @ beta, 1)

You could then go on to use y in further computation. It might not be identified though, because there are infinite combinations of X, alpha, and beta that produce the same y. It’s hard to say more on that without domain knowledge, but just be aware that you might have to use quite informative priors on something like this.

10

Thanks for the example. But what if the X and y are already modeled?

X = pm.Normal('X', mu = hyper_x_mu, 1, shape=(10))
y = pm.Normal('y', mu = hyper_y_mu, sigma = hyper_y_sigma, shape=(10))

Where hyper_x_mu, hyper_y_mu and hyper_y_sigma are RVs from upstream models.

Should I do this… or?

residuals = pm.Normal('residuals', mu=0, sigma=100, shape=(10))
alpha = pm.Normal('alpha', 0, 1)
beta = pm.Normal('beta', 0, 1)
should_be_zero = ((y + residuals) - alpha) / beta - 1 
#use potential to constrain 'should_be_zero' to zero?
11

I guess I would try to avoid defining a model that way. If you want y to be a function of X, link them directly via regression. Something like y = pm.Normal('y', mu = hyper_y_mu + beta_y * X, sigma = hyper_y_sigma).

My advice would be that any time you run into a situation where you want to impose some kind of constraint, instead think about how you can re-write the problem to directly define what you want. Reach for the hard constraints only in the cases where it’s really absolutely impossible to just directly write down what you want.

1 Like