I am working on a hierarchical bayesian model, which can be seen as logistic regression. Specifically, there are 4 variables in my model, which are xa, xb, pa, and pb. The sigmoid function in my model is
P(A) = \dfrac{1}{1+\exp\left(-\epsilon \big(u(x_a)-u(x_b) + \Phi(\psi(p_a) - \psi(p_b))\big)\right)},
where
u(x) = \dfrac{1}{\alpha}(1-\exp(-\alpha x))
\psi(p) = p^\gamma
\Phi(\psi(pa) - \psi(pb)) = \dfrac{\rho}{\beta}\left(1-\exp\Big(-\beta\big(\frac{\psi(pa) - \psi(pb)}{\theta}\big)^\theta\Big)\right).
The intervals for each parameter are \epsilon \in [0, +\infty), \alpha \in (0, +\infty), \gamma \in (0,1], \rho \in (0, +\infty), \beta \in (0, +\infty), \theta \in (1, +\infty).
My code is
import pymc3 as pm
import theano.tensor as tt
# The u() function
def utility(alpha, x):
value = (1 - tt.exp(-1 * alpha * x) ) / alpha
return value
# The \psi() function
def spf(gamma, op):
sp = op ** gamma
return sp
# The \Phi() function
def atr_diff(rho, beta, theta, spf_diff):
value = rho/beta * (1 - tt.exp(-1 * beta * (spf_diff / theta) ** theta) )
return value
with pm.Model() as model1:
mu_gamma = pm.Normal('mu_gamma', mu=0, sigma=1)
sigma_gamma = pm.Uniform('sigma_gamma', 0.01, 40)
gamma_phi_s_offset = pm.Normal('gamma_phi_s_offset', mu=0, sigma=1, shape=n_subject)
gamma_phi_s = pm.Deterministic('gamma_phi_s', mu_gamma + gamma_phi_s_offset * sigma_gamma)
gamma_log_s = pm.Normal.dist(mu=0, sigma=1).logcdf(gamma_phi_s)
gamma_s = pm.Deterministic('gamma_s', tt.exp(gamma_log_s))
mu_alpha = pm.Uniform('mu_alpha', -5, 5)
sigma_alpha = pm.Uniform('sigma_alpha', 0.01, 40)
alpha_s = pm.Lognormal('alpha_s', mu=mu_alpha, sigma=sigma_alpha, shape=n_subject)
mu_rho = pm.Normal('mu_rho', mu=0, sigma=10)
sigma_rho = pm.Uniform('sigma_rho', 0.01, 100)
rho_s = pm.Lognormal('rho_s', mu=mu_rho, sigma=sigma_rho, shape=n_subject)
mu_beta = pm.Normal('mu_beta', mu=0, sigma=10)
sigma_beta = pm.Uniform('sigma_beta', 0.01, 100)
beta_s = pm.Lognormal('beta_s', mu=mu_beta, sigma=sigma_beta, shape=n_subject)
mu_theta = pm.Normal('mu_theta', mu=0, sigma=10)
sigma_theta = pm.Uniform('sigma_theta', 0.01, 1000)
theta_s = pm.Lognormal('theta_s', mu=mu_theta, sigma=sigma_theta, shape=n_subject)
mu_epsilon = pm.Normal('mu_epsilon', mu=0, sigma=10)
sigma_epsilon = pm.Uniform('sigma_epsilon', 0.01, 100)
epsilon_s = pm.Lognormal('epsilon_s', mu=mu_epsilon, sigma=sigma_epsilon, shape=n_subject)
utility_diff = utility(alpha_s[subject_idx], xa) - utility(alpha_s[subject_idx], xb)
sprob_diff = spf(gamma_s[subject_idx], pa) - spf(gamma_s[subject_idx], pb)
diff = utility_diff + atr_diff(rho_s[subject_idx], beta_s[subject_idx], theta_s[subject_idx]+1, sprob_diff)
p = pm.math.invlogit(epsilon_s[subject_idx] * diff)
likelihood = pm.Bernoulli('likelihood', p=p, observed=choice)
The result of model1.check_test_point() is
mu_gamma -0.92
sigma_gamma_interval__ -1.39
gamma_phi_s_offset -2.76
mu_alpha_interval__ -1.39
sigma_alpha_interval__ -1.39
alpha_s_log__ -11.74
mu_rho -3.22
sigma_rho_interval__ -1.39
rho_s_log__ -14.49
mu_beta -3.22
sigma_beta_interval__ -1.39
beta_s_log__ -14.49
mu_theta -3.22
sigma_theta_interval__ -1.39
theta_s_log__ -21.40
mu_epsilon -3.22
sigma_epsilon_interval__ -1.39
epsilon_s_log__ -14.49
likelihood -208.95
Name: Log-probability of test_point, dtype: float64
There is no -inf or nan in this result, but I have got bad initial energy when I used NUTS. This is an essential part of my Master thesis, and please help me. THANK YOU!