Как сохранить выходные данные цикла при пересчете значений параметров для одновременных уравнений в R? - PullRequest
Новая
       24

Как сохранить выходные данные цикла при пересчете значений параметров для одновременных уравнений в R?

3 голосов
/ 07 октября 2019

У меня есть код, который оценивает значения параметров для 4 одновременных уравнений. Мне особенно интересно хранить все комбинации параметров, когда a + b (хранится в результатах $ ab) больше 3000. Если оно больше 3000, то оно кодируется как «Да». Я хочу написать цикл for, который будет перебирать код, чтобы проверить, есть ли a + b> 3000, и сохранить соответствующие значения. Затем я хочу, чтобы программа зациклилась 1000 раз и сохранила значения параметров для соответствующего «Да». Я пытаюсь сохранить вывод, но он не дает мне никаких результатов. 

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)

Ответы [ 3 ]

4 голосов
/ 17 октября 2019

a+b равно (x * t - 2 * x) + (v - x^2), что является просто квадратом в x, поэтому вы можете аналитически решить a+b>3000 для x, v и t.

Неравенство равно x^2 + (2-t)x + (3000-v) < 0.

Подставим T = 2-t и V = 3000-v, затем x^2 + Tx + V < 0.

Для того, чтобы любые значения были меньше нуля, он должен иметь два действительных корня,что означает, что T^2 - 4V > 0 - то есть V < (T^2)/4. (https://en.wikipedia.org/wiki/Quadratic_formula)

Учитывая T и V, удовлетворяющие этому неравенству, значения x, для которых a+b>3000 - это значения между корнями квадратичного, то есть |2x+T| < sqrt(T^2 - 4V).

Так что, если вы хотите выбрать значения, которые удовлетворяют условию, будет просто зациклить диапазон значений t, выбрать значение v, которое удовлетворяет V < (T^2)/4, а затем выбрать x, который попадает в соответствующий диапазон.

Вот один из способов ... 

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)
2 голосов
/ 19 октября 2019

Основано, в основном, на вашем описании, вот решение по тидиверу: 

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" ...

Создано в 2019-10-19 пакетом представлением (v0.3.0)

Дайте мне знать, если я неправильно понимаю цель.

2 голосов
/ 17 октября 2019

Вот решение , основанное на грубой силе - оцениваются все комбинации x и t, а также комбинации x и v: 

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
...