In [1]:
import pymc3 as pm
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt

%matplotlib inline
%config InlineBackend.figure_format = 'retina'
%qtconsole --colors=linux
plt.style.use('ggplot')

from matplotlib import gridspec
from theano import tensor as tt
from scipy import stats

Chapter 19 - Number concept development

19.1 Knower-level model for Give-N

$$ \pi \sim \text{Dirichlet}\overbrace{(1,...,1)}^{15} $$ $$ \nu \sim \text{Uniform}(1,1000)$$ $$ z_{i} \sim \text{Categorical}(\frac{1}{6},...,\frac{1}{6})$$

$$ \pi^{\prime}_{ijk} \propto \begin{cases} \pi_k & \text{if $k \gt z_i$} \\ \nu \times \pi_{ijk} & \text{if $k \leq z_i$ and $k = q^g_{ij}$} \\ \frac{1}{\nu} \times \pi_{ijk} & \text{if $k \leq z_i$ and $k \neq q^g_{ij}$} \end{cases} $$
$$ a_{ij}^g \sim \text{Categorical}(\pi^{\prime}_{ij})$$

In [2]:
import scipy.io as sio
matdata = sio.loadmat('data/fc_given.mat')

ns = np.squeeze(np.int64(matdata['ns']))
nz = np.squeeze(np.int64(matdata['nz']))
gnq = np.squeeze(np.int64(matdata['gnq']))
gn = np.squeeze(np.int64(matdata['gn']))
ga = np.squeeze(np.int64(matdata['ga']))
gq = np.squeeze(np.int64(matdata['gq']))

ind5 = np.zeros((nz, gn, gn), dtype=int)
for i in range(nz):
    i1 = i+1
    for j in range(gn):
        j1 = j+1
        for k in range(gn):
            k1 = k+1
            # Will be 1 if Knower-Level is Same or Greater than Answer
            ind1 = int(i1-1 >= k1)
            # Will be 1 for the Possible Answer that Matches the Question
            ind2 = int(k1 == j1)
            # Will be 1 for 0-Knowers
            ind3 = int(i1 == 1)
            # Will be 1 for HN-Knowers
            ind4 = int(i1 == nz)
            ind5[i, j, k] = ind3 + ind4*(2+ind2) + (1-ind4)*(1-ind3)*(ind1*ind2+ind1+1)
            
In [3]:
ind5r = ind5 - 1
ga_obs = np.asarray(ga.flatten()-1, dtype=int)
gq_obs = np.asarray(gq.flatten()-1, dtype=int)
valid_ind = np.where(gq_obs != -1)[0]

with pm.Model() as model1:
    pi = pm.Dirichlet('pi', a=np.ones(gn), shape=gn)
    
    nu = pm.Uniform('nu', lower=1, upper=1000)
    nu_vec = tt.stack([1., 1./nu, nu])
    
    piprime = tt.mul(nu_vec[ind5r], pi)
    npiprime = piprime / tt.sum(piprime, axis=-1, keepdims=True)
    
    zi = pm.Categorical('zi', p=np.ones(nz)/nz, shape=ns)
    zi_vec = tt.repeat(zi, gq.shape[1])
    
    pi_ij = npiprime[zi_vec[valid_ind], gq_obs[valid_ind], :]

    aij = pm.Categorical('aij', p=pi_ij, observed=ga_obs[valid_ind])  #, shape=len(valid_ind))  #

    trace1 = pm.sample(3e3, njobs=2)
    
pm.traceplot(trace1, varnames=['pi', 'nu', 'zi']);
Assigned NUTS to pi_stickbreaking__
Assigned NUTS to nu_interval__
Assigned CategoricalGibbsMetropolis to zi
100%|██████████| 3500/3500.0 [00:46<00:00, 75.70it/s]
In [4]:
def plot_samplerstat(burnin,trace):
    # Sampler statistics
    accept = trace.get_sampler_stats('mean_tree_accept', burn=burnin)
    print('The accept rate is: %.5f' % (accept.mean()))
    diverge = trace.get_sampler_stats('diverging', burn=burnin)
    print('Diverge of the trace')
    print(diverge.nonzero())
    energy = trace.get_sampler_stats('energy', burn=burnin)
    energy_diff = np.diff(energy)
    sns.distplot(energy - energy.mean(), label='energy')
    sns.distplot(energy_diff, label='energy diff')
    plt.legend()
    plt.show()
    
burnin = 1000
plot_samplerstat(burnin,trace1)
The accept rate is: 0.83038
Diverge of the trace
(array([], dtype=int64),)
In [5]:
pitr = trace1[burnin:]['pi']
fig = plt.figure(figsize=(18, 6))
plt.bar(np.arange(15)+1, np.mean(pitr, axis=0), align='center')
plt.xlabel('Number')
plt.ylabel('Probability')
plt.xticks(np.arange(16))
plt.xlim([0, 16])
plt.show()
In [6]:
fig = plt.figure(figsize=(20, 12))
gs = gridspec.GridSpec(4, 5)
zitr = trace1[burnin:]['zi']
for i in range(ns):
    ax = plt.subplot(gs[i])
    bartmp = np.unique(zitr[:, i], return_counts=True)
    ax.bar(bartmp[0], bartmp[1])
    ax.set_title('Child %s'%(i+1))
    plt.xlim([-.9, 5.9])
    plt.xticks(np.arange(6), ['P', '1', '2', '3', '4', 'C'])
    
plt.subplot(gs[15])
plt.xlabel('Knower Level')
plt.ylabel('Posterior Mass')
plt.tight_layout()
plt.show()
In [7]:
# Generate posterior prediction for each subject
gqpred = np.tile(np.arange(gn)[np.newaxis, :], (ns, 1)).flatten()

# Generate posterior prediction for each knower level
zpred = np.tile(np.arange(gn)[np.newaxis, :], (nz, 1)).flatten()
z_vect = np.tile(np.arange(nz)[:, np.newaxis], (1, gn)).flatten()

nutr = trace1[burnin:]['nu']
tracelen = nutr.shape[0]
nsample = 500
predga = np.zeros((nsample, len(gqpred)))
predz = np.zeros((nsample, len(zpred)))
randlist = np.random.choice(tracelen, nsample)
caterandom = lambda p: np.asarray([np.random.choice(len(p1), p=p1) for p1 in p])
for i, idx in enumerate(randlist):
    pi1, nu1, zi1 = pitr[idx], nutr[idx], zitr[idx]
    
    nu_vec1 = np.stack([1., 1./nu1, nu1])
    piprime1 = np.multiply(nu_vec1[ind5r], pi1)
    npiprime1 = piprime1 / np.sum(piprime1, axis=-1, keepdims=True)
    
    zi_vec2 = np.repeat(zi1, gn)
    pi_ij_pred = npiprime1[zi_vec2, gqpred, :]
    predga[i, :] = caterandom(pi_ij_pred)

    zi_ij_pred = npiprime1[z_vect, zpred, :]
    predz[i, :] = caterandom(zi_ij_pred)
In [8]:
predga = np.reshape(predga, newshape=(nsample, ns, gn)).astype(int)

fig = plt.figure(figsize=(20, 12))
gs = gridspec.GridSpec(2, 3)
subjlist = np.asarray([15, 2, 4, 3, 10, 20])-1
for i, isbj in enumerate(subjlist):
    mattmp = np.squeeze(predga[:, isbj, :])
    obs = ga[isbj]-1
    qus = gq[isbj]-1
    msk = qus != -1
    img_ = np.zeros((gn, gn))
    for j in range(gn):
        bartmp = np.unique(mattmp[:, j], return_counts=True)
        img_[j, bartmp[0]] = bartmp[1]
    ax = plt.subplot(gs[i])
    ax.imshow(img_.T, cmap='viridis', origin='lower')
    ax.plot(qus[msk], obs[msk], 'o', ms=15, color='r', alpha=.5)
    ax.grid('off')
    ax.set_title('Child %s'%(isbj+1))

plt.subplot(gs[3])
plt.xlabel('Question')
plt.ylabel('Answer')
plt.tight_layout()
plt.show()
In [9]:
predz = np.reshape(predz, newshape=(nsample, nz, gn)).astype(int)

fig = plt.figure(figsize=(20, 12))
gs = gridspec.GridSpec(2, 3)
knowertype = ('PN-Knower', 'One-Knower', 'Two-Knower',
              'Three-Knower', 'Four-Knower', 'CP-Knower')
for i in range(nz):
    mattmp = np.squeeze(predz[:, i, :])
    img_ = np.zeros((gn, gn))
    for j in range(gn):
        bartmp = np.unique(mattmp[:, j], return_counts=True)
        img_[j, bartmp[0]] = bartmp[1]
    ax = plt.subplot(gs[i])
    ax.imshow(img_.T, cmap='viridis', origin='lower')
    ax.grid('off')
    ax.set_title(knowertype[i])

plt.subplot(gs[3])
plt.xlabel('Question')
plt.ylabel('Answer')
plt.tight_layout()
plt.show()

19.2 Knower-level model for Fast-Cards

$$ \pi \sim \text{Dirichlet}\overbrace{(1,...,1)}^{50} $$ $$ \nu \sim \text{Uniform}(1,1000)$$ $$ z_{i} \sim \text{Categorical}(\frac{1}{6},...,\frac{1}{6})$$

$$ \pi^{\prime}_{ijk} \propto \begin{cases} \pi_k & \text{if $k \gt z_i$} \\ \nu \times \pi_{ijk} & \text{if $k \leq z_i$ and $k = q^f_{ij}$} \\ \frac{1}{\nu} \times \pi_{ijk} & \text{if $k \leq z_i$ and $k \neq q^f_{ij}$} \end{cases} $$
$$ a_{ij}^f \sim \text{Categorical}(\pi^{\prime}_{ij})$$

In [10]:
fnq = np.squeeze(np.int64(matdata['fnq']))
fa = np.squeeze(np.int64(matdata['fa']))
fq = np.squeeze(np.int64(matdata['fq']))
fn = 50

find5 = np.zeros((nz, gn, fn), dtype=int)
for i in range(nz):
    i1 = i+1
    for j in range(gn):
        j1 = j+1
        for k in range(fn):
            k1 = k+1
            # Will be 1 if Knower-Level is Same or Greater than Answer
            find1 = int(i1-1 >= k1)
            # Will be 1 for the Possible Answer that Matches the Question
            find2 = int(k1 == j1)
            # Will be 1 for 0-Knowers
            find3 = int(i1 == 1)
            # Will be 1 for HN-Knowers
            find4 = int(i1 == nz)
            find5[i, j, k] = find3 + find4*(2+find2) + (1-find4)*(1-find3)*(find1*find2+find1+1)
In [11]:
find5r = find5 - 1
fa_obs = np.asarray(fa.flatten()-1, dtype=int)
fq_obs = np.asarray(fq.flatten()-1, dtype=int)
valid_ind2 = np.where(fq_obs != -1)[0]

with pm.Model() as model2:
    pi = pm.Dirichlet('pi', a=np.ones(fn), shape=fn)
    
    nu = pm.Uniform('nu', lower=1, upper=1000)
    nu_vec = tt.stack([1., 1./nu, nu])
    
    piprime = tt.mul(nu_vec[find5r], pi)
    npiprime = piprime / tt.sum(piprime, axis=-1, keepdims=True)
    
    zi = pm.Categorical('zi', p=np.ones(nz)/nz, shape=ns)
    zi_vec = tt.repeat(zi, fq.shape[1])
    
    pi_ij = npiprime[zi_vec[valid_ind2], fq_obs[valid_ind2], :]

    aij = pm.Categorical('aij', p=pi_ij, observed=fa_obs[valid_ind2])  #, shape=len(valid_ind))  #

    trace2 = pm.sample(3e3, njobs=2)
    
pm.traceplot(trace2, varnames=['pi', 'nu', 'zi']);
Assigned NUTS to pi_stickbreaking__
Assigned NUTS to nu_interval__
Assigned CategoricalGibbsMetropolis to zi
100%|██████████| 3500/3500.0 [03:03<00:00, 18.64it/s]