Simulating data from given additive Bayesian networks


#1

Hello,

I need to simulate data from additive Bayesian networks (ABN) with parameters that will be given to Stan with an adjacency matrix. All variables in ABNs are modeled with an appropriate regression model.
Considering the network A --> B --> C where A and C are normally distributed variables and B is binary, the model would look as follows:

A ~ Normal(μ_A, σ^2_A), with μ_A = β_{A,0}

B ~ Bernoulli(π), with logit(π) = β_{B,0} + β_{B,1} * A

C ~ Normal(μ_C, σ^2_C), with μ_C = β_{C,0} + β_{C,1} * B

I got this to work with the following Stan code

data{
int<lower = 1> N;
matrix[N, N] X;
vector[N] sigma;
}
model{
}
generated quantities {
real A;
int B;
real C;
real piB;
A = normal_rng(X[1, 1], sigma[1]);
piB = 1/(1+exp(-(X[2,2]+X[1,2]*A)));
B = bernoulli_rng(piB);
C = normal_rng(X[3,3] + X[2,3]*B, sigma[3]);
}

where N = 3, sigma=(1, 0, 1.5) and X=
(0, 0.8, 0
0, 2, 1.5
0, 0, 0)
are given as data to Stan.

However in the general case, this doesn’t seem to work where we need some form of matrix-vector multiplication. In a purely normal case, it would be something along the lines of this:

functions {
vector delete_element(vector x, int m) {
int N;
vector[num_elements(x)-1] y;
N=num_elements(x);
for (i in 1:(m-1)) {
y[i]=x[i];
}
for (i in (m+1):N) {
y[i-1]=x[i];
}
return y;
}
}
data{
int<lower = 1> N;
matrix[N, N] X;
vector[N] sigma;
}
model {}
generated quantities {
vector[N] nodes;
for (i in 1:N) {
nodes[i] = normal_rng(X[i, i] + dot_product(delete_element(to_vector(X[i, ]), i),
delete_element(nodes, i)), sigma[i]);
}

This doesn’t work, since it includes undefined variables. Even when they are multiplied by zero, Stan gives:

normal_rng: Location parameter is nan, but must be finite!

I have also tried things with specified parameter and model block but could not find a solution so far. Any help would be much appreciated.


#2

Yeah so it looks like in the line:

nodes[i] = normal_rng(X[i, i] + dot_product(delete_element(to_vector(X[i, ]), i), delete_element(nodes, i)), sigma[i]);

You use elements of nodes before they are defined to anything. I think by default Stan fills things with NaNs (someone confirm this?).

Any IEEE 754 operation on a NaN returns a NaN. Did you mean to initialize the nodes with zeros? Something like:

vector[N] nodes = to_vector(0, N);

That help?


#3

I think you want rep_vector instead of to_vector, but otherwise exactly what @bbbales2 says


#4

That helped.
Thank you very much.