How to store output of a loop when resampling of parameter values for simultaneous equations in R?

user_6820344 0
Viewed 191

I have a code that evaluates parameter values for 4 simultaneous equations. I'm particularly interested to store all the parameter combinations when a + b (stored in results$ab) is greater than 3000. If it is greater than 3000, then it's coded as "Yes". I want to write a for-loop that will loop through the code to check for if a + b > 3000 and store the corresponding values. Then, I want the program to loop 1000 times, and store the parameter values for corresponding "Yes". I am trying to store the output, but it's not yielding me any results. 

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)


  values <- list(x = x, t = t, v = v, w = w, n = n, f = f)
  for(i in 1:1000){
    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)
    )

    for(i in 1:1000){
    samples <- 10000
    values.sampled <- lapply(values, sample, samples)
    results[i] <- sapply(eqs, eval, envir = values.sampled)
    results[i] <- data.frame(results)
    results$ab[i] <- results$a[i] + results$b[i]
    results$binary[i] <- ifelse(results$ab[i] > 3000, "Yes","No")
    output[i] <- results[results$binary=="Yes",]

  }

what <- as.list(output)
differential-equations r for-loop loops
3 Answers

a+b is equal to (x * t - 2 * x) + (v - x^2), which is just a quadratic in x, so you can solve a+b>3000 analytically for x, v and t.

The inequality is x^2 + (2-t)x + (3000-v) < 0.

Substitute T = 2-t and V = 3000-v, then x^2 + Tx + V < 0.

For this to have any values less than zero, it must have two real roots, which implies that T^2 - 4V > 0 - i.e V < (T^2)/4. (https://en.wikipedia.org/wiki/Quadratic_formula)

Given T and V satisfying this inequality, the values of x for which a+b>3000 are those between the roots of the quadratic, i.e. |2x+T| < sqrt(T^2 - 4V).

So if you want a selection of values that satisfy the condition, it should be straightforward to loop through a range of values of t, pick a value of v that satisfies V < (T^2)/4, and then pick an x that falls in the appropriate range.

Here is one way...

t <- 1:1000
T <- 2 - t
V <- sapply((T ^ 2) / 4, function(z) runif(min = 0, max = z, n = 1)) #assumes V>0 (???)
v <- 3000 - V
x <- (sapply(sqrt(T ^ 2 - 4 * V), function(z) runif(min = -z, max = z, n = 1)) - T) / 2
ab <- (x * t - 2 * x) + (v - x ^ 2) #all >3000 (except for t=2, where ab=3000 exactly)
user_77274290

Here's a data.table solution that relies on brute force - all combinations of x and t as well as combinations of x and v are evaluated:

library(data.table)
# create data
dt <- as.data.table(replicate(6, seq(from = 0.001, to = 1000, by = 1), simplify = F))
names(dt) <- c('x', 't','v','w','n','f')

# if your data.frame has all the combinations you want:
dt[, lapply(eqs, eval, envir = .SD)][, a_b := a + b][]

# This does all of the combinations for eqs a and b - takes 22 seconds for 10,000 rows
dt[, {x_t = CJ(x,t)
      x_v = CJ(x,v)
      a_b = eval(eqs$a, envir = x_t) + eval(eqs$b, envir = x_v)
      output = a_b > 3000L
      list(x = x_t[[1]], t = x_t[[2]], v = x_v[[2]], ab = a_b, output = output)
      }
   ][output == T,]

              x       t       v       ab output
     1:   3.001 754.001 754.001 3001.750   TRUE
     2:   3.001 755.001 755.001 3005.751   TRUE
     3:   3.001 756.001 756.001 3009.752   TRUE
     4:   3.001 757.001 757.001 3013.753   TRUE
     5:   3.001 758.001 758.001 3017.754   TRUE
    ---                                        
479045: 992.001 998.001 998.001 4966.005   TRUE
479046: 992.001 999.001 999.001 5959.006   TRUE
479047: 993.001 998.001 998.001 3977.004   TRUE
479048: 993.001 999.001 999.001 4971.005   TRUE
479049: 994.001 999.001 999.001 3981.004   TRUE
user_113550660

Based mostly on your description, here's a tidyverse solution:

library(tidyverse)

n_samples_per_trial <- 1e4
n_trials            <- 1e3
variable_ranges     <- seq(from = 0.0001, to = 1000, by = 0.1)
cutoff              <- 3e3

result_list <- rerun(.n = n_trials, {

  rerun(.n = 6, {
    sample(variable_ranges, n_samples_per_trial, replace = TRUE)
  }) %>% 
    as_tibble(.name_repair = ~ c("x", "t", "v", "w", "n", "f")) %>% 
    mutate(a      = x * t - 2 * x,
           b      = v - x^2,
           c      = x - w*t - t*t,
           d      = (n - f)/t,
           ab     = a + b,
           binary = if_else(ab > cutoff, "Yes", "No")) %>% 
    filter(binary == "Yes")

})

result_list %>% head(2) %>% glimpse()
#> List of 2
#>  $ :Classes 'tbl_df', 'tbl' and 'data.frame':    4840 obs. of  12 variables:
#>   ..$ x     : num [1:4840] 720.7 22.9 44.7 20.7 10.6 ...
#>   ..$ t     : num [1:4840] 794 325 768 400 531 ...
#>   ..$ v     : num [1:4840] 172 962 853 901 264 ...
#>   ..$ w     : num [1:4840] 842.8 274.3 18.9 957.2 321.8 ...
#>   ..$ n     : num [1:4840] 605 322 671 890 965 ...
#>   ..$ f     : num [1:4840] 397 557 299 691 825 ...
#>   ..$ a     : num [1:4840] 570795 7397 34222 8245 5604 ...
#>   ..$ b     : num [1:4840] -519236 437 -1145 473 152 ...
#>   ..$ c     : num [1:4840] -1298899 -194750 -603673 -543387 -452411 ...
#>   ..$ d     : num [1:4840] 0.261 -0.723 0.486 0.497 0.264 ...
#>   ..$ ab    : num [1:4840] 51558 7834 33077 8718 5756 ...
#>   ..$ binary: chr [1:4840] "Yes" "Yes" "Yes" "Yes" ...
#>  $ :Classes 'tbl_df', 'tbl' and 'data.frame':    4815 obs. of  12 variables:
#>   ..$ x     : num [1:4815] 400 806 410 477 169 ...
#>   ..$ t     : num [1:4815] 992 938 893 573 257 ...
#>   ..$ v     : num [1:4815] 453 773 456 279 933 ...
#>   ..$ w     : num [1:4815] 778 254 189 183 784 ...
#>   ..$ n     : num [1:4815] 629 811 411 529 667 ...
#>   ..$ f     : num [1:4815] 903 330 995 908 340 ...
#>   ..$ a     : num [1:4815] 395543 754765 365563 271843 43214 ...
#>   ..$ b     : num [1:4815] -159307 -649186 -167726 -226773 -27764 ...
#>   ..$ c     : num [1:4815] -1754435 -1117341 -966696 -431818 -267523 ...
#>   ..$ d     : num [1:4815] -0.277 0.513 -0.654 -0.662 1.275 ...
#>   ..$ ab    : num [1:4815] 236236 105579 197837 45070 15450 ...
#>   ..$ binary: chr [1:4815] "Yes" "Yes" "Yes" "Yes" ...

Created on 2019-10-19 by the reprex package (v0.3.0)

Let me know if I misunderstand the goal.

user_56934870
Related