Jacobian adjustment thing

jacobian.py

#%%

import matplotlib.pyplot as plt
import pystan
import numpy

model_code = """
data {
  int<lower=1> N;
  vector[N] y;
}

parameters {
  real a;
  real<lower = 0.0> b;
}

model {
  y ~ normal(a, b);
}
"""

sm1 = pystan.StanModel(model_code = model_code)

#%%

model_code = """
data {
  int<lower=1> N; // Number of single samples
  vector[N] y;
}

parameters {
  real a;
  real<lower=0.0> b;
}

model {
  vector[N] yh;
  
  yh = (y - a) / b;
  
  yh ~ normal(0.0, 1.0);
  
  target += N * log(abs(1 / b));
}
"""

sm2 = pystan.StanModel(model_code = model_code)

#%%
N = 100

y = 3.0 + 5.0 * numpy.random.randn(N)

print "First model"
print sm1.sampling({'N' : N, 'y' : y})

print "Second model"
print sm2.sampling({'N' : N, 'y' : y})