# Fitting simple (binomial) model in PyMC - slow convergence

I am trying to use pymc to estimate parameters for a model. As I am not familiar with the methodology I first generated the data myself and tried to derive the posterior distribution using pymc.

The data consists of a (time-) series of the number of missed payments $D$ for each quarter. For each quarter, the probability of a client missing a payment is a random variable, $p(x,\alpha)$, that depends on two parameters, $\alpha$ which is an unknown constant (in time) that I'm trying to estimate, and $x \sim \mathcal{N}(0,1)$, which is a random variable and drawn for each quarter.

Edit: for clarity, the function $p(x,\alpha)$ below describes the probability of a client to miss a payment. This probability is normally 3% per quarter but fluctuates a bit however, based on a correlation between clients ($\alpha$) and a random number. Please note each quarter is independent from the other. Only $x$ makes this a random variable. I have coded this as a deterministic in PyMC as its value it determined completely by $\alpha$ and $x$.

The actual probability of missing a payment in a quarter is described by the following function:

$$p(x,\alpha) = \Phi \left( \frac{ \Phi^{-1}(rate) - \alpha x }{\sqrt{1-\alpha^2}}\right)$$

where $\Phi$ is the normal cumulative distribution function and $\Phi^{-1}$ its inverse.

To generate a data sample for 50 quarters I have used this code (with the sample $\alpha = 0.37$ ).

sample_alpha = 0.37
const_z = 0.03 # basis prob of missing payment (3%)
nr_clients = 10000
nr_quarters = 50

Xs = np.random.normal(loc=0.0, scale=1.0,size=nr_quarters)

def p(x,a):
return sp.norm.cdf( (sp.norm.ppf(const_z) - sample_alpha*Xs)/(np.sqrt(1-sample_alpha*sample_alpha)) )

sample_mp = np.random.binomial(n=nr_clients, p=p(Xs,sample_alpha))

py.xlabel('quarter')
py.ylabel('nr of missed payments')
py.plot(sample_mp)


Below is shown a random realization of this process. To give more insight: every quarter the random variable $x$ and the actual $\alpha$ produce a certain probability of missing a payment. What is shown below for each quarter is then a drawing from a random binomial variable where the probability of success is this probability of missing a payment.

Now I'm trying to estimate the posterior distribution of $\alpha$ from this realization. I am using the following pymc code:

xfactor = pm.Normal('xfactor', mu=0, tau=1/(1*1), size = sample_mp.shape[0])

alpha = pm.Uniform("alpha", 0.01,0.99)

@pm.deterministic
def P(x=xfactor,a=alpha):
return sp.norm.cdf( (sp.norm.ppf(const_z) - a*x)/(np.sqrt(1-a*a)) )

missed_payments = pm.Binomial('defaults',n=10000,p=P, value=sample_mp,observed=True)
model = pm.Model([xfactor, alpha, P, missed_payments])


And then try to run it to get some feeling of what is going on (I haven't set burn-in and thinning at this point):

mcmc = pm.MCMC(model)
mcmc.sample(30000)
pm.Matplot.plot(mcmc)


The results look like this:

When I increase the nr of samples the plot is similar. I understand from this that $\alpha$ doesn't seem to converge, and the autocorrelation plot looks problematic too, but I am uncertain how to resolve the issue.

• Can you explain your model (esp. sp.norm.cdf part)? What is the parameter alpha and what it governs? – Vladislavs Dovgalecs Feb 12 '17 at 23:31
• Also, since your data is time-series, RandomWalk would be perhaps a better fit (for real-world data). But you're right, this model should fit the synthetic dataset where observations are i.i.d. – Vladislavs Dovgalecs Feb 12 '17 at 23:33
• I think Negative Binomial distribution would fit better the problem. – Vladislavs Dovgalecs Feb 13 '17 at 2:33
• @xeon: I have added the information in the post. I will have a look at RandomWalk; although since the realizations of p are independent from eachother I'm not sure how this would work out. – Letranger2014 Feb 13 '17 at 8:03
• As for the negative binomial: the number of clients is fixed (10,000) and the nr of clients that miss a payment fluctuates per quarter (base rate = 3%). For extra info: alpha here governs an intrinsic correlation between clients, so a higher alpha results in a higher p(x,a), and thus for the same x, a higher alpha means a higher p(x,a). I was thinking this (pymc-devs.github.io/pymc/…) might be of interest but I've been stuck on it for a day or 2 now. Any input is very welcome.. – Letranger2014 Feb 13 '17 at 8:10

Here is my shot at the problem in PyMC3. I can be wrong how the model is built, so please correct me where I am wrong.

The data are 50 observations (50 binomial draws) that are i.i.d. This assumptions is strong one.

I set the true parameter value (p_true=0.37) and set number of Bernoulli trials to 10,000. With these parameters set, I draw (K=50) samples from the Binomial distribution.

import pymc3 as pm
import matplotlib.pyplot as plt

from scipy.stats import binom

p_true = 0.37
n = 10000
K = 50

X = binom.rvs( n=n, p=p_true, size=K )
print( X )

model = pm.Model()

with model:
p = pm.Beta( 'p', alpha=2, beta=2 )
y_obs = pm.Binomial( 'y_obs', p=p, n=n, observed=X )
step = pm.Metropolis()
trace = pm.sample( 10000, step=step, progressbar=True )

pm.traceplot( trace )
plt.show()

• Thanks a lot for your help. I will take some time to study your version. There seems to be 1 misunderstanding which is that the p_true as in your model is not actually fixed: it is also a random variable. In my generation code above you can see that as I first produce the x, that lead to a different p_true for each drawing. Please note that the drawings indeed are independent as the x's are independent. – Letranger2014 Feb 13 '17 at 8:31
• @Letranger2014 Thanks for your feedback! It looks like you DO want to have time-series as your "p" changes from quarter to quarter. Here you would want to use RandomWalk. – Vladislavs Dovgalecs Feb 13 '17 at 9:46
• Indeed, P changes from quarter to quarter. In fact, p is in [0,1] and its pdf can be derived analytically -- it's just somewhat of an ugly function involving Erf and inverse Erf). I therefore opted for sampling x and producing p as a result. I havent succeeded in trying RandomWalk as the feature doesn't seem to be documented very well. – Letranger2014 Feb 15 '17 at 14:22
• Take a look at this example of RandomWalk github.com/pymc-devs/pymc3/blob/… and pymc-devs.github.io/pymc3/notebooks/survival_analysis.html – Vladislavs Dovgalecs Feb 15 '17 at 22:37
• @Letranger2014 I would try replacing the "p = pm.Beta(...)" with "p = pm.GaussianRandomWalk(...)". This would model p at time t conditioned on p at time t-1: prob(p_t | p_t-1) – Vladislavs Dovgalecs Feb 15 '17 at 22:39