Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
def test_pymc3_convert_dists():
"""Just a basic check that all PyMC3 RVs will convert to and from Theano RVs."""
tt.config.compute_test_value = "ignore"
theano.config.cxx = ""
with pm.Model() as model:
norm_rv = pm.Normal("norm_rv", 0.0, 1.0, observed=1.0)
mvnorm_rv = pm.MvNormal("mvnorm_rv", np.r_[0.0], np.c_[1.0], shape=1, observed=np.r_[1.0])
cauchy_rv = pm.Cauchy("cauchy_rv", 0.0, 1.0, observed=1.0)
halfcauchy_rv = pm.HalfCauchy("halfcauchy_rv", 1.0, observed=1.0)
uniform_rv = pm.Uniform("uniform_rv", observed=1.0)
gamma_rv = pm.Gamma("gamma_rv", 1.0, 1.0, observed=1.0)
invgamma_rv = pm.InverseGamma("invgamma_rv", 1.0, 1.0, observed=1.0)
exp_rv = pm.Exponential("exp_rv", 1.0, observed=1.0)
halfnormal_rv = pm.HalfNormal("halfnormal_rv", 1.0, observed=1.0)
beta_rv = pm.Beta("beta_rv", 2.0, 2.0, observed=1.0)
binomial_rv = pm.Binomial("binomial_rv", 10, 0.5, observed=5)
dirichlet_rv = pm.Dirichlet("dirichlet_rv", np.r_[0.1, 0.1], observed=np.r_[0.1, 0.1])
poisson_rv = pm.Poisson("poisson_rv", 10, observed=5)
bernoulli_rv = pm.Bernoulli("bernoulli_rv", 0.5, observed=0)
betabinomial_rv = pm.BetaBinomial("betabinomial_rv", 0.1, 0.1, 10, observed=5)
categorical_rv = pm.Categorical("categorical_rv", np.r_[0.5, 0.5], observed=1)
multinomial_rv = pm.Multinomial("multinomial_rv", 5, np.r_[0.5, 0.5], observed=np.r_[2])
negbinomial_rv = pm.NegativeBinomial("negbinomial_rv", 10.2, 0.5, observed=5)
# Convert to a Theano `FunctionGraph`
fgraph = model_graph(model)
def fit(self, X, Y, n_samples=10000, tune_steps=1000, n_jobs=4):
with pm.Model() as self.model:
# Priors
std = pm.Uniform("std", 0, self.sps, testval=X.std())
beta = pm.StudentT("beta", mu=0, lam=self.sps, nu=self.nu)
alpha = pm.StudentT("alpha", mu=0, lam=self.sps, nu=self.nu, testval=Y.mean())
# Deterministic model
mean = pm.Deterministic("mean", alpha + beta * X)
# Posterior distribution
obs = pm.Normal("obs", mu=mean, sd=std, observed=Y)
## Run MCMC
# Find search start value with maximum a posterior estimation
start = pm.find_MAP()
# sample posterior distribution for latent variables
trace = pm.sample(n_samples, njobs=n_jobs, tune=tune_steps, start=start)
# Recover posterior samples
self.burned_trace = trace[int(n_samples / 2):]
input_true_py = np.zeros(input_true.getSize())
for d in range(input_true.getSize()):
input_true_py[d] = input_true.get(d)
# add errors to measurements.
np.random.seed(0)
measurements = np.zeros((numTimeSteps, numMeasurements))
for m in range(numMeasurements):
for t in range(numTimeSteps):
measurements[t, m] = trueValues[t] + \
(np.random.randn(1) * sigma_true)
# The probabilistic model
with pm.Model() as prob_model:
# Priors for unknown model parameters
all_params = pm.Uniform('parameter', lower=0.5, upper=1.5, shape=dim)
# sigma = pm.HalfNormal('sigma', sd=0.1)
sigma = pm.Lognormal('sigma', mu=-1, sd=1)
# Forward model
ode_sol = my_ODEop(all_params)
forward = ode_sol
# Likelihood
Y_obs = []
for m in range(numMeasurements):
Y_obs.append(pm.Normal('Y_obs_%i' % m, mu=forward,
sd=sigma, observed=measurements[:, m]))
if stepType == 'Metropolis':
step = pm.Metropolis()
elif stepType == 'NUTS':
def get_garch_model():
r = np.array([28, 8, -3, 7, -1, 1, 18, 12], dtype=np.float64)
sigma1 = np.array([15, 10, 16, 11, 9, 11, 10, 18], dtype=np.float64)
alpha0 = np.array([10, 10, 16, 8, 9, 11, 12, 18], dtype=np.float64)
shape = r.shape
with Model() as garch:
alpha1 = Uniform('alpha1', 0., 1., shape=shape)
beta1 = Uniform('beta1', 0., 1 - alpha1, shape=shape)
mu = Normal('mu', mu=0., sigma=100., shape=shape)
theta = tt.sqrt(alpha0 + alpha1 * tt.pow(r - mu, 2) +
beta1 * tt.pow(sigma1, 2))
Normal('obs', mu, sigma=theta, observed=r)
return garch
return p_grid, posterior
points = 20
w, n = 6, 9
p_grid, posterior = posterior_grid_approx(points, w, n)
plt.figure()
plt.plot(p_grid, posterior, 'o-', label='success = {}\ntosses = {}'.format(w, n))
plt.xlabel('probability of water', fontsize=14)
plt.ylabel('posterior probability', fontsize=14)
plt.title('{} points'.format(points))
plt.legend(loc=0);
plt.show()
data = np.repeat((0, 1), (3, 6))
with pm.Model() as normal_aproximation:
p = pm.Uniform('p', 0, 1)
w = pm.Binomial('w', n=len(data), p=p, observed=data.sum())
mean_q = pm.find_MAP()
std_q = ((1/pm.find_hessian(mean_q, vars=[p]))**0.5)[0]
mean_q['p'], std_q
norm = stats.norm(mean_q, std_q)
prob = .89
z = stats.norm.ppf([(1-prob)/2, (1+prob)/2])
pi = mean_q['p'] + std_q * z
pi
# analytical calculation
w, n = 6, 9
x = np.linspace(0, 1, 100)
if self.correction_name in problem_config.hierarchicals:
nhierarchs = len(self.get_unique_stations())
param = problem_config.hierarchicals[self.correction_name]
logger.info(
'Estimating time shift for each station...')
kwargs = dict(
name=self.correction_name,
shape=nhierarchs,
lower=num.repeat(param.lower, nhierarchs),
upper=num.repeat(param.upper, nhierarchs),
testval=num.repeat(param.testvalue, nhierarchs),
transform=None,
dtype=tconfig.floatX)
try:
station_corrs_rv = Uniform(**kwargs)
except TypeError:
kwargs.pop('name')
station_corrs_rv = Uniform.dist(**kwargs)
self.hierarchicals[self.correction_name] = station_corrs_rv
else:
nhierarchs = 0
' hyperparameter(s): %s!' % hp_name)
if not num.array_equal(hyperpar.lower, hyperpar.upper):
dimension = hyperpar.dimension * ndata
kwargs = dict(
name=hyperpar.name,
shape=dimension,
lower=num.repeat(hyperpar.lower, ndata),
upper=num.repeat(hyperpar.upper, ndata),
testval=num.repeat(hyperpar.testvalue, ndata),
dtype=tconfig.floatX,
transform=None)
try:
hyperparams[hp_name] = Uniform(**kwargs)
except TypeError:
kwargs.pop('name')
hyperparams[hp_name] = Uniform.dist(**kwargs)
modelinit = False
n_hyp += dimension
else:
logger.info(
'not solving for %s, got fixed at %s' % (
hyperpar.name,
utility.list_to_str(hyperpar.lower.flatten())))
hyperparams[hyperpar.name] = hyperpar.lower
if len(hyperparameters) > 0:
z = np.array([45, 63, 58, 64, 58, 63, 51, 60, 59, 47, 63, 61, 60, 51, 59, 45,
61, 59, 60, 58, 63, 56, 63, 64, 64, 60, 64, 62, 49, 64, 64, 58, 64, 52, 64, 64,
64, 62, 64, 61, 59, 59, 55, 62, 51, 58, 55, 54, 59, 57, 58, 60, 54, 42, 59, 57,
59, 53, 53, 42, 59, 57, 29, 36, 51, 64, 60, 54, 54, 38, 61, 60, 61, 60, 62, 55,
38, 43, 58, 60, 44, 44, 32, 56, 43, 36, 38, 48, 32, 40, 40, 34, 45, 42, 41, 32,
48, 36, 29, 37, 53, 55, 50, 47, 46, 44, 50, 56, 58, 42, 58, 54, 57, 54, 51, 49,
52, 51, 49, 51, 46, 46, 42, 49, 46, 56, 42, 53, 55, 51, 55, 49, 53, 55, 40, 46,
56, 47, 54, 54, 42, 34, 35, 41, 48, 46, 39, 55, 30, 49, 27, 51, 41, 36, 45, 41,
53, 32, 43, 33])
condition = np.repeat([0,1,2,3], nSubj)
# Specify the model in PyMC
with pm.Model() as model:
# define the hyper-hyperparameters for kappa
mean_gamma = pm.Uniform('mean_gamma', 0, 30)
sd_gamma = pm.Uniform('sd_gamma', 0, 30)
s_kappa = mean_gamma**2/sd_gamma**2
r_kappa = mean_gamma/sd_gamma**2
# define the hyperparameters
kappa = pm.Gamma('kappa', s_kappa, r_kappa)
mu = pm.Beta('mu', 1, 1, shape=ncond)
# define the prior
theta = pm.Beta('theta', mu[condition] * kappa, (1 - mu[condition]) * kappa, shape=len(z))
# define the likelihood
y = pm.Binomial('y', p=theta, n=N, observed=z)
trace = pm.sample(2000, tune=1000)
## Check the results.
burnin = 0 # posterior samples to discard
## Print summary for each trace
#pm.df_summary(trace[burnin:])
# define the hyperpriors
a1_SD_unabs = pm.StudentT('a1_SD_unabs', mu=0, lam=0.001, nu=1)
a1_SD = abs(a1_SD_unabs) + 0.1
a1tau = 1 / a1_SD**2
a2_SD_unabs = pm.StudentT('a2_SD_unabs', mu=0, lam=0.001, nu=1)
a2_SD = abs(a2_SD_unabs) + 0.1
a2tau = 1 / a2_SD**2
a1a2_SD_unabs = pm.StudentT('a1a2_SD_unabs', mu=0, lam=0.001, nu=1)
a1a2_SD = abs(a1a2_SD_unabs) + 0.1
a1a2tau = 1 / a1a2_SD**2
# define the priors
sigma = pm.Uniform('sigma', 0, 10) # y values are assumed to be standardized
tau = 1 / sigma**2
a0 = pm.Normal('a0', mu=0, tau=0.001) # y values are assumed to be standardized
a1 = pm.Normal('a1', mu=0 , tau=a1tau, shape=Nx1Lvl)
a2 = pm.Normal('a2', mu=0 , tau=a2tau, shape=Nx2Lvl)
a1a2 = pm.Normal('a1a2', mu=0 , tau=a1a2tau, shape=[Nx1Lvl, Nx2Lvl])
b1 = pm.Deterministic('b1', a1 - tt.mean(a1))
b2 = pm.Deterministic('b2', a2 - tt.mean(a2))
b1b2 = pm.Deterministic('b1b2', a1a2 - tt.mean(a1a2))
mu = a0 + b1[x1] + b2[x2] + b1b2[x1, x2]
# define the likelihood
yl = pm.Normal('yl', mu=mu, tau=tau, observed=z)
100, 99])
y = pd.DataFrame({
'value': np.r_[drug, placebo],
'group': np.r_[['drug']*len(drug), ['placebo']*len(placebo)]
})
y_mean = y.value.mean()
y_std = y.value.std() * 2
sigma_low = 1
sigma_high = 10
with pm.Model():
group1_mean = pm.Normal('group1_mean', y_mean, sd=y_std)
group2_mean = pm.Normal('group2_mean', y_mean, sd=y_std)
group1_std = pm.Uniform('group1_std', lower=sigma_low, upper=sigma_high)
group2_std = pm.Uniform('group2_std', lower=sigma_low, upper=sigma_high)
lambda_1 = group1_std**-2
lambda_2 = group2_std**-2
nu = pm.Exponential('ν_minus_one', 1/29.) + 1
pm.StudentT('drug', nu=nu, mu=group1_mean, lam=lambda_1, observed=drug)
pm.StudentT('placebo', nu=nu, mu=group2_mean, lam=lambda_2, observed=placebo)
diff_of_means = pm.Deterministic('difference of means', group1_mean - group2_mean)
pm.Deterministic('difference of stds', group1_std - group2_std)
pm.Deterministic(
'effect size', diff_of_means / np.sqrt((group1_std**2 + group2_std**2) / 2))
pm.sample(draws=20000, cores=4, chains=4,
progressbar=False, compute_convergence_checks=False)