Write a loop to select all combination of variable values generating positive equation values in R

Question

I have the following four equations (a,b,c,d), with several different variables (x,t,v,w,n,f). My goal would be to try and find all sets of variable values that would generate all positive (and non-zero) numbers for equations(a,b,c,d). A regular loop would just go through each number of the sequence generated and systematically check if it generates a positive value or not. I want it to pick up random numbers from each sequence and test it against the others in R. For example (x=8, t = 2.1,v=13,w=1,n=10,f=1) is a possible set of combinations.

Please do not suggest analytically solving for these and then finding out values. These are simply representations of equations I'm dealing with. The equations I have are quite complex, and more than 15 variables.

#Equations
a <- x * t - 2*x
b <- v - x^2 
c <- x - w*t - t*t 
d <- (n - f)/t

x <- seq(from = 0.0001, to = 1000, by = 0.1)
t <- seq(from = 0.0001, to = 1000, by = 0.1)
v <- seq(from = 0.0001, to = 1000, by = 0.1)
w <- seq(from = 0.0001, to = 1000, by = 0.1)
n <- seq(from = 0.0001, to = 1000, by = 0.1)
f <- seq(from = 0.0001, to = 1000, by = 0.1)

Mischiefmaker

Answer 1

For a start, it might be better to organize your equations and your probe values into lists:

set.seed(1222)

values <- list(x = x, t = t, v = v, w = w, n = n, f = f)

eqs <- list(
  a = expression(x * t - 2 * x),
  b = expression(v - x^2), 
  c = expression(x - w*t - t*t), 
  d = expression((n - f)/t)
)

Then we can define a number of samples to take randomly from each probe vector:

samples <- 3
values.sampled <- lapply(values, sample, samples)

$x
[1] 642.3001 563.1001 221.3001

$t
[1] 583.9001 279.0001 749.1001

$v
[1] 446.6001 106.7001   0.7001

$w
[1] 636.0001 208.8001 525.5001

$n
[1] 559.8001  28.4001 239.0001

$f
[1] 640.4001 612.5001 790.1001

We can then iterate over each stored equation, evaluating the equation within the "sampled" environment:

results <- sapply(eqs, eval, envir = values.sampled)

            a          b         c          d
[1,] 373754.5 -412102.82 -711657.5 -0.1380373
[2,] 155978.8 -316975.02 -135533.2 -2.0935476
[3,] 165333.3  -48973.03 -954581.8 -0.7356827

From there you can remove any value that is 0 or less:

results[results <= 0] <- NA

Estrade

For the last bit would this work? resultsfinal <- ifelse(results > 0, results, NULL) – Mischiefmaker

Do I need to set.seed, or do I specify the values as I have done in the question? – Mischiefmaker

set.seed is only there so that this example is reproducible. It'll make it so that the calls to sample always return the same results. In practice you probably don't want that, though as the analysis is finalized it comes in handy for making sure you always have identical results (i.e., for publication). – Estrade

This is a tad bit confusing. without setting seed I don't think the computer understands what "eval" is doing. Here is my code getting stuck rextester.com/ZJRHQ60262 – Mischiefmaker

In your linked code, you create a variable in the main environment called r but do not copy this into the list version. So the lapply uses the unsampled copy of r, which has too many values. This makes results a list where the results of each equation has a different number of values, which cannot be index-replaced (with [ ]). My advice would be to use seq directly in the list, i.e., list(a = seq(1, 10, 1).... – Estrade

Answer 2

If every independent value can take on the same value (e.g. seq(from = 0.0001, to = 1000, by = 0.1)), we can approach this with much greater rigor and avoid the possibility of generating duplicates. First we create a masterFun that is essentially a wrapper for all of the functions you want to define:

masterFun <- function(y) {
    ## y is a vector with 6 values
    ## y[1] -->> x
    ## y[2] -->> t
    ## y[3] -->> v
    ## y[4] -->> w
    ## y[5] -->> n
    ## y[6] -->> f

    fA <- function(x, t) {x * t - 2*x}
    fB <- function(v, x) {v - x^2}
    fC <- function(x, w, t) {x - w*t - t*t}
    fD <- function(n, f, t) {(n - f)/t}

    ## one can easily filter out negative
    ## results as @jdobres has done.

    c(a = fA(y[1], y[2]), b = fB(y[3], y[1]), 
      c = fC(y[1], y[4], y[2]), d = fD(y[5], y[6], y[2]))
}

Now, using permuteSample, which is capable of generating random permutations of a vector and subsequently applying any given user defined function to each permutation, from RcppAlgos (I am the author), we have:

## Not technically the domain, but this variable name
## is concise and very descriptive
domain <- seq(from = 0.0001, to = 1000, by = 0.1)

library(RcppAlgos)

          ## number of variables ... x, t, v, w, n, f
          ##           ||
          ##           \/
permuteSample(domain, m = 6, repetition = TRUE,
                 n = 3, seed = 123, FUN = masterFun)
[[1]]
            a              b              c              d 
218830.316100 -608541.146040 -310624.596670      -1.415869 

[[2]]
            a              b              c              d 
371023.322880 -482662.278860 -731052.643620       1.132836 

[[3]]
             a               b               c               d 
18512.60761001 -12521.71284001 -39722.27696002     -0.09118721

In short, the underlying algorithm is capable of generating the n^th lexicographical result, which allows us to apply a mapping from 1 to "# of total permutations" to the permutations themselves. For example, given the permutations of the vector 1:3:

permuteGeneral(3, 3)
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    1    3    2
[3,]    2    1    3
[4,]    2    3    1
[5,]    3    1    2
[6,]    3    2    1

We can easily generate the 2^nd and the 5^th permutation above without generating the first permutation or the first four permutations:

permuteSample(3, 3, sampleVec = c(2, 5))
     [,1] [,2] [,3]
[1,]    1    3    2
[2,]    3    1    2

This allows us to have a more controlled and tangible grasp of our random samples as we can now think of them in a more familiar way (i.e. a random sample of numbers).

If you actually want to see which variables were used in the above calculation, we simply drop the FUN argument:

permuteSample(domain, m = 6, repetition = TRUE, n = 3, seed = 123)
         [,1]     [,2]     [,3]     [,4]     [,5]     [,6]
[1,] 780.7001 282.3001 951.5001 820.8001 289.1001 688.8001
[2,] 694.8001 536.0001  84.9001 829.2001 757.3001 150.1001
[3,] 114.7001 163.4001 634.4001  80.4001 327.2001 342.1001

Australorp