stan_linear_regression.py

import pystan
import numpy as np

stan_code = """
data {
    int<lower=0> N;
    int<lower=0> K;
    matrix[N,K] x;
    vector[N] y;
}
parameters {
    real alpha;
    vector[K] beta;
    real<lower=0> sigma;
}
model {
    y ~ normal(x * beta + alpha, sigma);
}
"""

sample_x1 = np.random.randn(1000)
sample_x2 = np.random.randn(1000)
sample_y = 2 * sample_x1 + 3 * sample_x2 + np.random.randn(sample_x1.shape[0])

sample_X = np.array(zip(sample_x1, sample_x2))
stan_model = pystan.StanModel(model_code=stan_code)
op = stan_model.sampling(data=dict(x=sample_X, y=sample_y, N=len(sample_x1), K=sample_X.shape[1]), chains=1, iter=5000)
print op