Fitting a basic uniform distribution -- what am I missing?

Modeling
1

I am attempting to fit (what I think is) a very straightforward model, I have y that I believe to be generated from a uniform distribution, so:

y \sim \textrm{Uniform}(\alpha, \beta)

Where \alpha and \beta are unknown. I know the data to be positive, so I defined the bounds of the uniform distribution to be positive_ordered and attempted to put priors on where I believe the bounds to be located.

My model code is below:

data {
    int<lower=0> N; // num rows
    vector[N] y;
}
parameters {
    positive_ordered[2] bounds;
}
model {
    bounds[1] ~ normal(1180, 10);
    bounds[2] ~ normal(1220, 10);
    y ~ uniform(bounds[1], bounds[2]);
}

When attempting to fit this, I get the following error:

Rejecting initial value:
  Log probability evaluates to log(0), i.e. negative infinity.
  Stan can't start sampling from this initial value.

I have attempted to change my priors, use default priors, change the parameter type, etc. but have been unable to find the bug thatā€™s not allowing this model to be fit. I feel like I have a fundamental misunderstanding of what Stan is doing here. My bet is that Iā€™m not specifying my priors correctly.

I would appreciate any help at tracking down this bug, and more importantly, any explanation as to what Iā€™m missing here would be very helpful. Iā€™m new to Stan and sometimes I run into these issues that I find absolutely bewildering that I know are born from ignorance that makes even knowing how to start debugging a challenge.

2

I think you need to set you inits to equal 1 in the model call. Stan generates a random number for starting between -2 to 2.
What is your model call?

3

Thank you for the response!

Iā€™m just using the default model call with pystan, so the line is:

sm.sampling(data={'N': n, 'y': y})

Where sm is the StanModel compiled from the above code. I see there is an init parameter that defaults to init='random'. I imagine may be the right path to go down.

4

Yeah if you can set that to init = 1 and see what happens.

5

No change with init=1 ā€“ digging into the documentation for that parameter, I found:

init : {0, '0', 'random', function returning dict, list of dict}, optional
    Specifies how initial parameter values are chosen: 0 or '0'
    initializes all to be zero on the unconstrained support; 'random'
    generates random initial values; list of size equal to the number
    of chains (`chains`), where the list contains a dict with initial
    parameter values; function returning a dict with initial parameter
    values. The function may take an optional argument `chain_id`.

Trying 0, ā€˜0ā€™ and ā€˜randomā€™ get me nowhere, Iā€™m trying to dig into the option to pass a list of dicts into the sampling. Iā€™ve gotten here:

init_bounds = [
    {'bounds': (1100, 1300)},
    {'bounds': (1170, 1270)},
    {'bounds': (1150, 1250)},
    {'bounds': (1190, 1300)},
]

fit = sm.sampling(data=data, init=init_bounds)

Which gets me closer, it looks like 3 of the 4 chains do properly sample, but the last does not. This feels like Iā€™m doing things in a rather arbitrary fashion, surely there is a better way to approach this?

6

Oh yeah sorry I forgot that you have to set each chain. Can you set each one to have a bound of 0,1. Sorry do most my Stan work in R 'cause ecology.

7

Went ahead and changed the bounds to (0, 1) for each, got the error:

Initialization between (-2, 2) failed after 1 attempts.
 Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.

Suppose it could be a pystan issue? I can switch over to R and give it shot there.

8

Hereā€™s some code that should work

data {
    int<lower=0> N; // num rows
    vector[N] y;
}
parameters {
    real<lower=0,upper=min(y)> lower_bound;
    real<lower=max(y)> upper_bound;
}
model {
    lower_bound ~ normal(1180, 10);
    upper_bound ~ normal(1220, 10);
    y ~ uniform(lower_bound, upper_bound);
}
3 Likes
9

@nhuurre this worked! My question is why? Does the positive_ordered type not encode this exact same information?

1 Like
10

Stanā€™s sampler starts by picking a random parameter value that is consistent with the declared constraint (but only the ones in parameters block; transformed parameters and model are ignored). Because this initial value is completely arbitrary it usually very far from both prior and posterior.
After choosing the initial value the HMC algorithm begins gradient ascend like random walk towards the target distribution.

But this only works if the initial point has a gradient. The uniform distribution rejects a parameter value that puts any y outside of the bounds.

positive_ordered only says that the upper bound is larger than the lower bound. In particular, Stan tries to initialize it to a value like [1,2] which is not consistent with your data.

2 Likes
11

Thank you for the explanation!

So if Iā€™m understanding correctly ā€“ by setting the parameter limits based on the data, y, Stanā€™s sampler is guaranteed to pick an initialization point that is compatible with the uniform model and thus can appropriately compute a gradient?

12

Yes. In general the parameter limits should select the region where probability is nonzero so that the sampler knows not to explore ā€œimpossibleā€ parameter values.

2 Likes