Hi everyone, I’ve been working through the *Bayesian Methods for Hackers* book by Cameron Davidson-Pilon but instead of PyMC I’m writing everything in R and Stan.

Things have gone rather smoothly so far but I came across an interesting example that I’m having trouble vectorizing. I believe I have the model specified correctly without vectorization but here’s the info so far:

The original model in PyMC2 is below:

```
def euclidean_distance(x, y):
return np.sqrt(((x - y) ** 2).sum(axis=1))
def f_distance(gxy_pos, halo_pos, c):
# foo_position should be a 2-d numpy array
return np.maximum(euclidean_distance(gxy_pos, halo_pos), c)[:, None]
def tangential_distance(glxy_position, halo_position):
# foo_position should be a 2-d numpy array
delta = glxy_position - halo_position
t = (2 * np.arctan(delta[:, 1] / delta[:, 0]))[:, None]
return np.concatenate([-np.cos(t), -np.sin(t)], axis=1)
import pymc as pm
# set the size of the halo's mass
mass_large = pm.Uniform("mass_large", 40, 180, trace=False)
# set the initial prior position of the halos, it's a 2-d Uniform dist.
halo_position = pm.Uniform("halo_position", 0, 4200, size=(1, 2))
@pm.deterministic
def mean(mass=mass_large, h_pos=halo_position, glx_pos=data[:, :2]):
return mass / f_distance(glx_pos, h_pos, 240) *\
tangential_distance(glx_pos, h_pos)
ellpty = pm.Normal("ellipcity", mean, 1. / 0.05, observed=True,
value=data[:, 2:])
mcmc = pm.MCMC([ellpty, mean, halo_position, mass_large])
map_ = pm.MAP([ellpty, mean, halo_position, mass_large])
map_.fit()
mcmc.sample(200000, 140000, 3)
```

Here is my un-vectorized version of the model written in Stan:

```
functions {
real f_distance(row_vector gxy_pos, vector halo_pos, real c) {
return fmax(distance(gxy_pos, halo_pos), c);
}
vector tangential_distance(row_vector glxy_position, vector halo_position) {
real xdiff = glxy_position[1] - halo_position[1];
real ydiff = glxy_position[2] - halo_position[2];
real t = (2 * atan(ydiff / xdiff));
vector[2] tdist;
tdist[1] = -cos(t);
tdist[2] = -sin(t);
return tdist;
}
}
data {
int<lower=0> n;
matrix[n, 2] cart_pos; // x,y coordinates of galaxy position
matrix[n, 2] ellip_pos; // a measure of ellipticity
}
parameters {
real<lower=40.0,upper=180.0> exp_mass_large;
vector<lower=0,upper=4200.0>[2] halo_position;
}
transformed parameters {
real mass_large = log(exp_mass_large); // one large halo
}
model {
vector[2] mu;
for (i in 1:n) {
mu = mass_large / f_distance(cart_pos[i,:], halo_position, 240.0) *
tangential_distance(cart_pos[i,:], halo_position);
ellip_pos[i,1] ~ normal(mu[1], 0.05); // x-coordinate
ellip_pos[i,2] ~ normal(mu[2], 0.05); // y-coordinate
}
}
```

The way I wrote the model above appears to sample okay. There is definitely some bi-modal behavior with the marginal distributions that I didn’t see in the book but overall the results look similar. However I have attempted (and failed) to write the above program in a vectorized way similar to the book. I’ve been stumbling with using `vector`

data types or arrays of `vectors`

and when to use element-wise multiplication/division (`.*`

, `./`

) or even if functions like `fmax`

can be applied element-wise.

Any pointers would be appreciated!