How is Bayesian Network specified? I want to generate a random Bayesian Network (with about 500 variables), and then pull about 500 samples from it. What is the standard "recipe" for a Bayesian Network that will let me do that. Will starting off with $V^2$ binary matrix representing the network's connections be enough?
 A: You can do this with the bnlearn R package. Creating a random DAG is  straightforward. Assigning random parameters to local (at each node) not so much, as I am not aware of a bnlearn function that does this. Thus, I provide below two functions to do this in a simple way. See below for all code and additional comments.
Install the package if you do not have it and then load it
# install.packages('bnlearn')
library(bnlearn)

Once you have bnlearn loaded, you can proceed. First, you create a random DAG with random.graph(). A number of algorithms are available for this and you can choose yours with the method argument.
nodes <- paste0('X', 1:500, sep = '')
g <- random.graph(nodes)

Once you have the DAG you need to decide the distributions of the local nodes.
In bnlearn, these can be either categorical or Gaussian. You can combine categorical and Gaussian nodes in a single network, with the constraint that the Gaussian ones cannot be parents of categorical ones.
Here I will assume all 500 local distributions are Gaussians. Each conditional local distribution then corresponds to a linear regression of the varaible given its parents. I provide two functions to fit the local distributions with random parameters. 
get_random_node_coefficients(): A function to assign random coefficients to a local (conditional) distribution of a node given its parents in the DAG. It simply samples the coefficients and the standard deviation of the linear regression from a standard normal, keeping the absolute value of the sampled standard deviation.
get_random_node_coefficients <- function(node, dag) { 
  parents <- parents(dag, node) 
  intercept <- rnorm(1)
  coef_parents <- rnorm(length(parents))
  coef <- c(intercept, coef_parents)
  names(coef) <- c("(Intercept)", parents) 
  # sd must be non-negative
  sd <- abs(rnorm(1))
  list(coef = coef, sd = sd)
} 

get_random_coefficients(): A function to get the coefficients of all nodes in an adequate format.
get_random_coefficients <- function(g) {
 coefs <- lapply(nodes(g), get_random_node_coefficients, g)
 names(coefs) <- nodes(g)  
 coefs
}  

I now use the above defined functions to specify the coefficients of the local linear regressions for the all nodes:
coefs <- get_random_coefficients(g)
g.fit <- custom.fit(g, dist = coefs)

Finally, you use rnb() to sample from the network. You need to create a sort of a dummy data object to tell rnb() the form of the data frame you want your results in.
data <- data.frame(matrix(numeric(0), nrow = 0, ncol = length(nodes))) 
sample <- rbn(g.fit, n = 500, data) 

And you have the sample.
head(sample)  
hist(sample$X1)

