Let the distinct values be $$X = \{x_1, x_2, \ldots, x_n\}.$$ Let the sample size be $m\, (=3).$ The ordered samples of distinct values are the subset
$$\Omega = \{(y_1, y_2, \ldots, y_m) \in X^m\mid y_1\lt y_2\lt \ldots \lt y_m\}.$$
The question seeks a probability distribution on $\Omega$ with probabilities
$$p(i_1,i_2,\ldots, i_m) = \Pr((y_1,\ldots, y_m) = (x_{i_1}, \ldots, x_{i_m}))$$
with given marginal distributions. That is, let $1 \le k \le m$ indicate the order of the value and let $P_k$ be its marginal distribution. For each $1 \le i \le n,$ the question specifies the probabilities
$$P_{k,i} = \Pr(\text{the } k^{\text{th}} \text{ smallest value is }x_i) = \sum p(i_1,\ldots, i_{k-1},\, i,\, i_{k+1}, \ldots, i_m)\tag{1}$$
where the sum is over all the $i_j,$ $j=1, \ldots, m$ except for $j=k.$
This can be expressed as a linear program for the values $p(i_1, i_2, \ldots, i_m)$ by maximizing their sum subject to the equality constraints $(1)$ and the inequality constraints $p(i_1,i_2, \ldots, i_m) \ge 0.$ It is possible there are no solutions or multiple solutions.
If a solution exists, then sampling is straightforward because the $p(i_1,\ldots,i_m)$ specify a discrete distribution on a finite set and there are efficient algorithms to sample discrete distributions. In R, for instance, use the built-in function sample or sample.int (after indexing the nonzero values of $p(i_1,\ldots,i_m)$ from $1$ to whatever).
I found an implementation of the Simplex algorithm in the boot package for R. According to it, there is no solution to the data in the question when the samples must consist of distinct values. But if we relax this condition and not insist they be distinct, there is a solution. It is given by the objects X (representing $\Omega$) and Prob in the R code:
> (cbind(X, Prob))
Var1 Var2 Var3 Prob
2 2 2 4 0.03
12 4 4 4 0.07
27 3 3 5 0.11
32 4 4 5 0.18
33 1 5 5 0.06
36 4 5 5 0.05
46 2 3 6 0.24
55 3 5 6 0.04
57 1 6 6 0.12
61 1 2 7 0.07
78 2 6 7 0.03
The first column assigns abstract identifiers to the data. The next three columns give one possible ordered $m$-element sample $(x_{i_1}, x_{i_2}, \ldots, x_{i_m})\in \Omega.$ The final column gives its probability $p(i_1,i_2,\ldots, i_m).$
For instance, the question requires the value $x_1 = 1$ to appear $25\%$ of the time. That is, $1$ should be the smallest value with a probability of $0.25.$ In the outcomes labeled 33, 57, and 61 $1$ is the smallest value and these outcomes contribute $0.06 + 0.12 + 0.07 = 0.25$ of the probability, as intended.
Here is a jittered scatterplot matrix of a thousand draws from this distribution:
The darkness of each square is proportional to the number of times it appeared in this simulation. The order restriction is visibly manifest: in the upper triangle of plots, all dots lie on or below the diagonal line $x=y,$ while in the lower triangle all dots lie on or above the diagonal.
To verify this is correct, I ran an experiment in which $167$ triples were sampled from this distribution and compared to the desired probabilities using a chi-squared test. That produced a p-value. Upon repeating this experiment for a total of 10,000 runs, the distribution of the p-values was essentially uniform (varying by an amount explainable by chance variation).
The following R script shows how to obtain the $P_{k,i}$ from a natural representation of the problem and includes the scatterplot matrix and the simulation code. As testimony to the efficiency of the algorithm, all the work done by this script (which includes generating a half million random variables) takes less than a second.
#
# Specify the problem.
#
D <- list(data.frame(value=1:4, prob=c(25,30,15,30)),
data.frame(value=2:6, prob=c(10,35,25,15,15)),
data.frame(value=4:7, prob=c(10,40,40,10)))
distinct <- FALSE # Enforce distinct values in ordered samples?
set.seed(17) # Optional: creates reproducible results
#
# Create a data frame of all possible ordered samples.
#
X <- do.call(expand.grid, lapply(D, function(X) X$value))
if (isTRUE(distinct)) {
X$Valid <- apply(X, 1, function(x) all(sort(x) == x) & length(unique(x)) == length(x))
} else {
X$Valid <- apply(X, 1, function(x) all(sort(x) == x))
}
X <- subset(X, X$Valid==TRUE)
X$Valid <- NULL
names(D) <- varnames <- names(X)
#
# Find a feasible set of probabilities (a joint distribution) corresponding
# to the observati0ns of `X`.
#
# This solves the linear program
# maximize sum(Prob)
# given A3 %*% Prob == b3, Prob >= 0.
# (Due to a quirk in `simplex`, the first equality is relaxed to an inequality.)
# Solutions have the given marginal probabilities (as expressed by `A3` and `b3`).
#
library(boot) # simplex()
A3 <- do.call(cbind, lapply(varnames, function(s)
sapply(D[[s]]$value, function(i) X[[s]] == i)
))
A3 <- ifelse(t(A3), 1, 0)
b3 <- unlist(lapply(D, function(Y) Y$prob / sum(Y$prob)))
# b3 <- b3 / sum(b3)
obj <- simplex(a=rep(1, nrow(X)), A1=A3, b1=b3, maxi=TRUE)
if(obj$solved <= 0) stop("No feasible solution found.")
if(sum((A3 %*% obj$soln - b3)^2) >= 1e-15) {
warning("The `simplex` solution appears incorrect.")
}
# -- Limit the sample space to the nonzero probabilities
Prob <- zapsmall(obj$soln)
Prob <- Prob / sum(Prob)
X <- X[Prob > 0, ]
Prob <- Prob[Prob > 0]
(cbind(X, Prob)) # Display the distribution
#
# Sample from this distribution repeatedly.
#
n.sample <- min(1e4, ceiling(5 / min(Prob)))
df <- length(Prob) - 1
xlim <- range(unlist(X)) + c(-1,1)
p.value <- replicate(1e4, {
#
# Sample from the multivariate distribution.
#
i <- sample.int(nrow(X), n.sample, replace=TRUE, prob=Prob)
#
# Compare the sample frequencies to the intended frequencies.
#
obs <- tabulate(i, nrow(X))
exp <- Prob * n.sample
chisq <- sum((obs - exp)^2 / exp)
pchisq(chisq, df, lower.tail=FALSE)
})
#
# Plot the chi-squared p-values that tested the sample frequencies.
#
hist(p.value, freq=FALSE)
abline(h=1, col="Red")
#
# Plot a sample.
#
Y <- X[sample.int(nrow(X), 1e3, replace=TRUE, prob=Prob), ]
pairs(Y[1:min(1e3, nrow(Y)), varnames] + runif(nrow(Y)*length(names(Y)), -0.15, 0.15),
xlim=xlim, ylim=xlim,
pch=16, col="#00000002")
