Right, so IIUC in a way you are interested in the function that maps the values group of parameters (node1_i,node2_i,node3_i) to the ouput value ouput_i right? I would start with a linear regression, something like: y = pm.Normal('observed', tt.dot(np.asarray([node1, node2, node3]), beta), sd=sigma, observed=output), with beta and sigma as the free parameters.