A procedure for sampling paths in a directed acyclic graph Consider a directed acyclic graph $\mathcal{G} = (\mathcal{N},\mathcal{E})$. Assume that the graph is very large (on the order of 10000 nodes and edges). 
Let there be a set of nodes termed starting nodes denoted by $\mathcal{N}_s\subseteq\mathcal{N}$ and a set of nodes termed terminal nodes denoted by $\mathcal{N}_t\subseteq\mathcal{N}$. Consider the set of all paths that start from a node in $\mathcal{N}_s$ and end in a node in $\mathcal{N}_t$. Denote this set by $\mathcal{P}$.
Goal: Derivation of a procedure that allows one to uniformly sample paths from the set $\mathcal{P}$ without having to construct $\mathcal{P}$ explicitly.
Ideas: Of course, listing all such paths is computationally intractable preventing one from just being able to sample from a table of all paths. I am wondering if one can construct a procedure that involves sampling edges of the graph at random (without replacement) until one obtains a path that satisfies the condition that it starts in $\mathcal{N}_s$ and ends in $\mathcal{N}_t$ such that the procedure samples paths uniformly from $\mathcal{P}$. 
I've sketched out such a procedure for a small graph below. The top (green) nodes are the nodes $\mathcal{N}_s$ whereas the bottom (purple) nodes are the nodes $\mathcal{N}_t$. At each step of the procedure a single edge is sampled at random (denoted by the red edges). A check is done to see if it is a path that satisfies the condition, if not, sample a new edge (without replacement) and check again. Repeat until we find a path. If the addition of a single edge results in multiple paths, pick one at random. We add this path to our set of sampled paths, call it $\mathcal{\bar P}$, and repeat the process.

Question: Does such a procedure result in uniformly sampled paths from $\mathcal{P}$? If not, are there procedures that do so without requiring one to construct $\mathcal{P}$?
 A: Here's an example where the output of this procedure is far from uniform:

In the top half of the graph, we have $n$ "diamonds", each consisting of $4$ edges. There are $2^n$ paths from left to right in the top half of the graph: for each diamond, we can choose one of two possible routes to take. There is also one path in the bottom half of the graph consisting of a single edge.
The bottom path has at least a $\frac1{4n+1}$ chance of being chosen because that is the chance we will choose its only edge in the very first step. (And actually, it gets at least $n$ chances to be chosen before any other path can be constructed, so it is the first path to be found a constant fraction of the time.) But we want it to have a $\frac1{2^n+1}$ chance of being chosen. So, for large $n$, its probability of being chosen will be exponentially higher than the chance we choose any of the other paths.

We can use rejection sampling to produce an algorithm that will work for any graph with a little bit of finetuning. But its performance depends heavily on properties of the graph.
Begin with a method we know to be biased. Choose a random starting node in $\mathcal N_s$. Then repeatedly choose a random edge out of the last vertex visited, until we arrive at a vertex in $\mathcal N_t$. This produces a path. (This method gets stuck if it arrives at a vertex not in $\mathcal N_t$ with out-degree $0$, so we should pre-process the graph by iteratively pruning such vertices, which cannot contribute to paths anyway.)
Even though this is biased, we know exactly how likely it is that any given path will be chosen: a path on vertices $v_0, v_1, v_2, \dots, v_k$ with $v_0 \in \mathcal N_s$ and $v_k \in \mathcal N_t$ is chosen with probability
$$
   p = \frac1{|\mathcal N_s|} \cdot \frac1{\deg^+(v_0)} \cdot \frac1{\deg^+(v_1)} \cdots \frac1{\deg^+(v_{k-1})}.
$$
We can also put some lower bound $p_{\text{min}}$ on the probability of choosing the least likely path. A trivial lower bound is $\frac1{|\mathcal N_s|}\prod_{v \in \mathcal N} \frac1{\deg^+(v)}$, but we can do better if we can put an upper bound on the length of a path or on the number of high-degree vertices each path visits. (In the example I drew above, it's easy to see that $p_{\text{min}} = \frac1{2^{n+1}}$.)
So the final step in the algorithm is that once we generate a path with the biased method, we actually return it as our output only with probability $\frac{p_{\text{min}}}{p}$, and start over otherwise. As a result, each path has the same probability of $p \cdot \frac{p_{\text{min}}}{p} = p_{\text{min}}$ of being returned at each iteration, and the sampling is uniform.
The better the lower bound $p_{\text{min}}$, the faster the algorithm will work. The speed is also affected by the distribution: if there are a few very unlikely paths that achieve the lower bound, then they will slow down the algorithm for all the other paths. But if a constant fraction of paths have a probability close to $p_{\text{min}}$ of being returned, then this algorithm will usually require only a few trials.
A: For simplicity I'll assume $\mathcal{N}_s=\{v_s\}$ and $\mathcal{N}_t=\{v_t\}$, but later I generalize the idea to the general case.
The algorithm:
Step 1: Calculate and store the #paths from each node $v \in \mathcal{N}$ to $v_t$. It can be done with a variation of the DFS algorithm, where in each visit we return the sum of the returned successor calls.
nr_paths_to_vt := Mapping(N -> int)  # assume each entry is initialized to 0
dfs(v):
    if v == v_t:
        nr_paths_to_vt[v_t] = 1
        return 1
    nr_paths_to_vt[v] = sum(dfs(s) for s in successors(v))
    return nr_paths_to_vt[v]
call dfs(v_s)

Note: It can be implemented more efficiently by storing intermediate results and avoiding re-visiting already-visited nodes. Here is an implementation example: https://gist.github.com/eladn/0cf26489ecad9101acfa9415dc21a365
Step 2: We are going to iteratively construct a path from $v_s$ to $v_t$, similarly to Misha's suggestion, but with probabilities over the out edges of each node.
p = [v_s] # the constructed path, which starts with the starting node 
while p[-1] != v_t:  # the path haven't reached the terminating node yet
    next_node_distr = [ (node=s, prob=nr_paths_to_vt(s)/nr_paths_to_vt(p[-1])) | for s in successors(p[-1]) ]
    next_node = sample_with_probs(next_node_distr)
    p = p || next_node

Time complexity: Each step requires $O( | \mathcal{N} | + | \mathcal{E} | )$.
Correctness: Can be shown in induction over the graph nodes, where the induction variable $n$ is the distance from the node $v$ to $v_t$.
Generalization to the original case: We can add two new nodes $v_s$ and $v_t$ to the graph, and add edges from $v_s$ to each $v \in \mathcal{N}_s$ and from each $v \in \mathcal{N}_t$ to $v_t$.
