I'm having some troubles with the implementation in pyMC of my probabilistic model.
Note: you can skip directly to the code section, if you're not interested on the use of the model.
The model itself is going to be used to extrapolate information about different the style of a 3D model (and the individual styles of each of its components) using a network with both observed and hidden random variables. Just FYI, the database includes only models of the same type (.: all tables)
I'm following the work done in the paper titles "A Probabilistic Model for Component Based Shape Synthesis", but trying to learn the different parameters and best structure without using EM but giving MCMC a try.
The model can be summarised with this figure (link).
The section 'L' is repeated for every type of components the models are formed (e.g.: if we still think of tables, we'll have 'tops' and 'legs' as sections).
N,D (discrete) and C (continuous) are the observed variables, taking as data some information about adjacency, geometric features, ... .
S and R (discrete) are the hidden variables, trying to represent respectively the style of each component of the model and the global style of the object.
The probability of the whole model (in its simplest form, I'm not considering additional learned link between variables) is factorised as a product of CPDs
$$ p(X) = P(R) \prod_{l \in L}[ P(S_l | R) \, P(N_l|R) \, P(C_l|S_l) \, P(D_l|S_l) ] $$
From what I understood from reading the paper, the conditional probabilities for discrete variables are represented as CPTs, whom parameters are estimated using Dirichlet distributions as priors $$ P(N=n|R=r) = q_{n|r} $$ $$ Q= \{q_{n|r}\} \sim Dirichlet $$
Instead the continuous variables are represented with multivariate normals with Normal-Wishart priors.
Actual start of the code (can skip directly here)
The code I came up with, is basically a mixture of Categorical nodes with Dirichlet Priors (for the discrete variables) and MvNormal (for the continuous ones).
This is a snippet which represents only a "slice" of the model, as I'm considering only a single component type instead of multiple ones.
I think the explanation is pretty straightforward, the value of the hidden categorical variables S and R select the different priors for the observed nodes, which should be then estimated using the data I provide. That should, hopefully, retain the information of possible different styles in said priors.
This doesn't actually happen... In fact, during the sampling, the values of the hidden variables S and R get "stuck" on a single value for all the process, leaving some of the nodes not used and just left in their "uncharacterized" state.
This is a plot of the model created from this code: link
I'm willing to share plots of the traces for the variables or data if needed. Unfortunately Cross Validated doesn't let me add more than 2 links in this post.
n_R = 2 # Number of different global styles
n_S = {1:2} #NUmber of different styles for component type
n_N = {1:5} #Max number of objects for that component
N_data = {} # Data for observed N nodes
C_data = {} # Data for observed C nodes
D_data = {} # Data for observed D nodes
D_feat_num = 35 #Max value for D data
C_feat_len = 14 #Feature lenght for C data
p_R = mc.Dirichlet( name ='p_R', theta = [1.0, ] * n_R)
R = mc.Categorical('R', p = p_R)
p_N1 = [ mc.Dirichlet( name ='p_N1_R%i' % j, theta = [1.0,] * (n_N[1] + 1)) for j in range(n_R) ]
@mc.observed
@mc.stochastic(dtype=int)
def N1(value = N_data[1], R = R, p_N1 = p_N1):
def logp(value, R, p_N1):
return mc.categorical_like(value, np.append(p_N1[R], ( 1.0 - np.sum(p_N1[R]) )))
def random(R, p_N1):
return mc.rcategorical(np.append(p_N1[R], ( 1.0 - np.sum(p_N1[R]) )))
p_S1 = [ mc.Dirichlet( name ='p_S1_R%i' % j, theta = [1.0,] * (n_S[1])) for j in range(n_R) ]
@mc.stochastic(dtype=int)
def S1(value = 0, R = R, p_S1 = p_S1):
def logp(value, R, p_S1):
return mc.categorical_like(value, np.append(p_S1[R], ( 1.0 - np.sum(p_S1[R]) )))
def random(R, p_S1):
return mc.rcategorical(np.append(p_S1[R], ( 1.0 - np.sum(p_S1[R]) )))
p_D1 = [ mc.Dirichlet( name ='p_D1_S%i' % j, theta = [1.0,] * D_feat_num ) for j in range(n_S[1]) ]
@mc.observed
@mc.stochastic(dtype=int)
def D1(value = D_data[1], S1 = S1, p_D1 = p_D1):
def logp(value, S1, p_D1):
return mc.categorical_like(value, np.append(p_D1[S1], ( 1.0 - np.sum(p_D1[S1]) )))
def random(S1, p_D1):
return mc.rcategorical(np.append(p_D1[S1], ( 1.0 - np.sum(p_D1[S1]) )))
tauC1_list = [ mc.Wishart('tauC1_S%i' % j, n = C_feat_len + 1, Tau = np.eye(C_feat_len) , trace = True, plot = False) for j in range(n_S[1])]
muC1_list = [ mc.Normal('muC1_S%i' % j, mu = 0.0, tau = 0.1, size = C_feat_len , trace = True, plot = False) for j in range(n_S[1])]
muC1 = mc.Lambda('muC1' , lambda muC1_list = muC1_list, S1 = S1: muC1_list[S1], trace = True, plot = False)
tauC1 = mc.Lambda('tauC1', lambda tauC1_list = tauC1_list, S1 = S1: tauC1_list[S1], trace = True, plot = False)
C1 = mc.MvNormalCov('C1', mu = muC1, C = tauC1, observed = True, value = C_data[1])
A = locals()
sampler = mc.MCMC(A)
sampler.sample(iter = 10000)
mod = mc.Model(A)
dag(mod, format="png", name="arch_toy_2", path="graph/")
directory = 'graph/prova_toy_2'
if( not os.path.isdir(directory) ):
os.mkdir(directory)
plot(sampler, path=directory)
