I need to sample angle uniformly which is simple enough:
parameters {
vector<lower=0, upper=1>[N] angle_raw;
}
transformed parameters {
vector[N] angle;
angle = acos(1-angle_raw); // will be uniformly distributed between 0 and pi/2
}
However, if I need the angle to be ordered, I lose the ability to place constraints. So, following the manual:
parameters {
ordered angle_rawer[N];
}
transformed parameters {
vector<lower=0, upper=1>[N] angle_raw;
vector[N] angle;
angle_raw = inv_logit(angle_rawer);
angle = acos(1-angle_raw); // will be uniformly distributed between 0 and pi/2
}
model {
angle_raw ~ uniform(0,1);
target += log(angle_raw) + log(1 - angle_raw); // jacobian adjustment.
}
Is this the correct way to handle this?
Sorry, no idea what the correct way is, just noticed the following line which might be a problem:
2 Likes
Oops. That was just a typo in the example! Thanks!
Ok, after hunting a bit (and understanding), I think I might have a solution. While I’m not sure if it works in the bigger model I’m building, I have tested a bit and I think it is the right path.
Starting from here: https://groups.google.com/forum/#!topic/stan-users/04GSu-ql3vM
I tested this with a little python code to make sure the general idea worked:
def simplex_sample(n):
"""
sample from an n-dim unit simplex
"""
sample = np.random.uniform(0,1,size=n-1)
sample = np.append(sample,1.)
sample.sort()
sample = np.append(0,sample)
diffs = np.diff(sample)
return diffs
def ordered_interval(n, a = 0, b=1.):
"""
sample an ordered (distinguishable) interval
"""
simplex = simplex_sample(n+1)
return a+ np.cumsum(simplex)[:n]*(b-a)
Which yields:
sample = ordered_interval(n=10000, a = 0, b= 1)
len(np.unique(sample)) # it equals 10000
And then I transform this to variable uniform in angle:
ang = np.rad2deg(np.arccos(1-sample))
Thus, the method should produce an ordered angle uniform in angle.
The stan code should be:
parameters {
simplex angle_rawer[N+1];
}
transformed parameters {
postive_ordered[N] angle_raw;
vector[N] angle;
angle_raw = head( cumulative_sum(angle_rawer), N );
angle = acos(1-angle_raw); // will be uniformly distributed between 0 and pi/2 in angle
}
model {
angle_raw ~ uniform(0,1);
}
And no jacobian adjustments needed (I believe! )
Update: This works great in the bigger model.
