我有一个代码可以评估4个联立方程的参数值。当a + b（存储在results $ 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)
[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值，然后选择一个[C0 ]的范围内。
这里是一种方法...
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)
解决方案-评估data.table和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
由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" ...
（v0.3.0）在2019-10-19创建
让我知道我是否误解了目标。