<jupyter_start><jupyter_text>We will model the number of deaths as a random sample from a binomial distribution, where $n$ is the number of rats and $p$ the probabbility of a rat dying. We are given $n=5$, but we believve that $p$ may be related to the drug dose $x$. As $x$ increases the number of rats dying seems to increase, and since $p$ is a probability, we use the following model:$y \sim Bin(n,p) \\logit(p) = \alpha + \beta x \\a \sim N(0,5) \\\beta \sim N(0,10)$where we set vague priors for αα and ββ, the parameters for the logistic model.<jupyter_code># define invlogit function
def invlogit(x):
    return pymc.exp(x) / (1 + pymc.exp(x))
# observed data
dose, rats, deaths

# define priors
alpha = pymc.Normal('alpha', mu=0, tau=1.0/5**2)
beta = pymc.Normal('beta', mu=0, tau=1.0/10**2)

# define likelihood
p = pymc.InvLogit('p', alpha + beta*dose)
y = pymc.Binomial('y_obs', n=rats, p=p, value=deaths, observed=True)

# inference
m = pymc.Model([alpha, beta, y])
mc = pymc.MCMC(m)
mc.sample(iter=21000, burn=1000, burn_till_tuned=True, thin=20)
pymc.Matplot.plot(mc)
pprint(mc.stats())
xp = np.linspace(-1, 1, 100)
a = alpha.stats()['mean']
b = beta.stats()['mean']
plt.plot(xp, invlogit(a + b*xp).value)
plt.scatter(dose, deaths/5, s=50);
plt.xlabel('Log dose of drug')
plt.ylabel('Risk of death')
# one should be able to get estimates of the line uncertainity
ilu = np.empty((1000, len(xp)), dtype=float)
for ii, v in enumerate(np.random.random_integers(0, len(alpha.trace[:])-1, 1000)):
    ilu[ii] = invlogit(alpha.trace[:][v] + beta.trace[:][v]*xp).value

ilu
xp = np.linspace(-1, 1, 100)
for v in ilu:
    plt.plot(xp, v, alpha=.01, c='r')


a = alpha.stats()['mean']
b = beta.stats()['mean']
plt.plot(xp, invlogit(a + b*xp).value)

a = alpha.stats()['quantiles'][50]
b = beta.stats()['quantiles'][50]
plt.plot(xp, invlogit(a + b*xp).value, c='y')


plt.scatter(dose, deaths/5, s=50);
plt.xlabel('Log dose of drug')
plt.ylabel('Risk of death')
alpha.stats()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1791936286709806, 'fraction_numerical': 0.030363364858138377, 'mean_word_length': 4.792507204610951, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '446245', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 853, 'n_tokens_neox': 793, 'n_words': 229, 'path': 'Logit/Logit 2.ipynb', 'repo_name': 'balarsen/pymc_learning', 'size': '2009'}