如何从参数集递归重新采样，直到R中满足所需的条件

Question

我有一个数据集，其中A1、A2、A3、B1、B2、B3作为状态变量和存储在值中的参数。我试图找到一个条件，当A1，A2，A3，B1，B2，B3都是阳性。我想循环这个代码，这样它就可以遍历所有可能的参数值组合，如果A1到B3为正，则将其存储在单独的输出中(尽管在R中可能需要几个小时(

library(Ryacas)
fromval <- 0.1
toval <- 500
byval <- 0.1
repval <- 5000
f <- seq(from = fromval, to = toval, by = byval)
f1 <- seq(from = fromval, to = toval, by = byval)
ep <- seq(from = fromval, to = toval, by = byval)
a1 <- seq(from = fromval, to = toval, by = byval)
h <- seq(from = fromval, to = toval, by = byval)
e1 <- seq(from = fromval, to = toval, by = byval)
m1 <- seq(from = fromval, to = toval, by = byval)
mp <- seq(from = fromval, to = toval, by = byval)
ab <- seq(from = fromval, to = toval, by = byval)
r <- seq(from = fromval, to = toval, by = byval)
fp <- seq(from = fromval, to = toval, by = byval)
R <- seq(from = fromval, to = toval, by = byval)
values <- list(f1=f1, ep=ep, a1=a1, h=h, e1=e1, m1=m1, mp=mp,f=f,ab=ab,r=r,fp=fp,R=R)
eqs <- list(A1= expression((a1*ep*f*f1*fp*R - e1*ep*fp^2*m1*(h + R) + ep^2*f1*(ab*fp - a1*mp*r)*(h + R) - sqrt(ep^2*(-4*a1*fp*mp*(e1*fp*m1 + ep*f1*r)*(h + R)*(a1*f*f1*R - e1*m1*(fp - a1*mp)*(h + R)) +  (ab*ep*f1*fp*(h + R) - e1*fp*m1*(fp - 2*a1*mp)*(h + R) + a1*f1*(f*fp*R + ep*mp*r*(h + R))) ^2)))/(2*a1*ep*f1*(e1*fp*m1 + ep*f1*r)*(h + R))),

A2= expression((ep*fp*(a1*f^2*f1^2*R^2 - e1*f*f1*m1*(fp - a1*mp)*R*(h + R) + e1^2*fp*m1^2*mp*(h + R)^2) + ep^2*f1*(h + R)*(ab*fp*(f*f1*R + e1*m1*mp*(h + R)) + mp*r*(-(a1*f*f1*R) + e1*m1*(2*fp - a1*mp)*(h + R))) + e1*h*m1*mp*sqrt(ep**2*(-4*a1*fp*mp*(e1*fp*m1 + ep*f1*r)*(h + R)* (a1*f*f1*R - e1*m1*(fp - a1*mp)*(h + R)) +  (ab*ep*f1*fp*(h + R) - e1*fp*m1*(fp - 2*a1*mp)*(h + R) + a1*f1*(f*fp*R + ep*mp*r*(h + R)))^2)) + f*f1*R*sqrt(ep^2*(-4*a1*fp*mp*(e1*fp*m1 + ep*f1*r)*(h + R)* (a1*f*f1*R - e1*m1*(fp - a1*mp)*(h + R)) +  (ab*ep*f1*fp*(h + R) - e1*fp*m1*(fp - 2*a1*mp)*(h + R) + a1*f1*(f*fp*R + ep*mp*r*(h + R)))**2)) + e1*m1*mp*R*sqrt(ep^2*(-4*a1*fp*mp*(e1*fp*m1 + ep*f1*r)*(h + R)* (a1*f*f1*R - e1*m1*(fp - a1*mp)*(h + R)) + (ab*ep*f1*fp*(h + R) - e1*fp*m1*(fp - 2*a1*mp)*(h + R) + a1*f1*(f*fp*R + ep*mp*r*(h + R)))^2)))/(2*ep^2*f1*fp*mp*(h + R)*(a1*f*f1*R - e1*m1*(fp - a1*mp)*(h + R)))),
A3= expression((ep^2*f1*(ab*fp + a1*mp*r)*(h + R) + ep*fp*(a1*f*f1*R - e1*m1*(fp - 2*a1*mp)*(h + R)) - sqrt(ep^2*(-4*a1*fp*mp*(e1*fp*m1 + ep*f1*r)*(h + R)*(a1*f*f1*R - e1*m1*(fp - a1*mp)*(h + R)) + (ab*ep*f1*fp*(h + R) - e1*fp*m1*(fp - 2*a1*mp)*(h + R) + a1*f1*(f*fp*R + ep*mp*r*(h + R))) ^2)))/(2*a1*fp*(e1*fp*m1 + ep*f1*r)*(h + R))),

B1 = expression((1/(2* a1* ep* f1* (e1* fp* m1 + ep* f1* r)* (h + R)))*(a1* ep* f* f1* fp* R - e1* ep* fp^2* m1* (h + R) + ep^2* f1* (ab* fp - a1* mp* r) *(h + R) +sqrt(ep^2* (-4* a1* fp* mp* (e1* fp* m1 + ep* f1* r)* (h + R)* (a1* f* f1* R - e1* m1* (fp - a1* mp)* (h + R)) + (ab* ep* f1* fp* (h + R) - e1* fp* m1* (fp - 2* a1* mp)* (h + R) + a1* f1* (f* fp* R + ep* mp* r* (h + R)))^2)))),

B2= expression(-(-(ep*fp*(a1*f^2*f1^2*R^2 - e1*f*f1*m1*(fp - a1*mp)*R*(h + R) + e1^2*fp*m1^2*mp*(h + R)^2)) - ep^2*f1*(h + R)*(ab*fp*(f*f1*R + e1*m1*mp*(h + R)) + mp*r*(-(a1*f*f1*R) + e1*m1*(2*fp - a1*mp)*(h + R))) + e1*h*m1*mp*sqrt(ep^2*(-4*a1*fp*mp*(e1*fp*m1 + ep*f1*r)*(h + R)*(a1*f*f1*R - e1*m1*(fp - a1*mp)*(h + R)) + (ab*ep*f1*fp*(h + R) - e1*fp*m1*(fp - 2*a1*mp)*(h + R) + a1*f1*(f*fp*R + ep*mp*r*(h + R)))^2)) + f*f1*R*sqrt(ep^2*(-4*a1*fp*mp*(e1*fp*m1 + ep*f1*r)*(h + R)*(a1*f*f1*R - e1*m1*(fp - a1*mp)*(h + R)) + (ab*ep*f1*fp*(h + R) - e1*fp*m1*(fp - 2*a1*mp)*(h + R) + a1*f1*(f*fp*R + ep*mp*r*(h + R)))^2)) + e1*m1*mp*R*sqrt(ep^2*(-4*a1*fp*mp*(e1*fp*m1 + ep*f1*r)*(h + R)*(a1*f*f1*R - e1*m1*(fp - a1*mp)*(h + R)) + (ab*ep*f1*fp*(h + R) - e1*fp*m1*(fp - 2*a1*mp)*(h + R) + a1*f1*(f*fp*R + ep*mp*r*(h + R)))^2)))/(2*ep^2*f1*fp*mp*(h + R)*(a1*f*f1*R - e1*m1*(fp - a1*mp)*(h + R)))),

B3= expression((ep^2*f1*(ab*fp + a1*mp*r)*(h + R) + ep*fp*(a1*f*f1*R - e1*m1*(fp - 2*a1*mp)*(h + R)) + sqrt(ep^2*(-4*a1*fp*mp*(e1*fp*m1 + ep*f1*r)*(h + R)*(a1*f*f1*R - e1*m1*(fp - a1*mp)*(h + R)) + (ab*ep*f1*fp*(h + R) - e1*fp*m1*(fp - 2*a1*mp)*(h + R) + a1*f1*(f*fp*R + ep*mp*r*(h + R)))^2)))/(2*a1*fp*(e1*fp*m1 + ep*f1*r)*(h + R)))
)

samples <- 5000
values.sampled <- lapply(values, sample, samples)
results <- sapply(eqs, eval, envir = values.sampled)
values.sampled.df <- data.frame(values.sampled)
results <- data.frame(results)
values.sampled.df$sl <- c(1:samples)
results$sl <- c(1:samples)
df <- merge(results,values.sampled.df,by="sl")
dfraw <- df
df$A1[df$A1 <= 0] <- NA
df$A2[df$A2 <= 0] <- NA
df$A3[df$A3 <= 0] <- NA
df$B1[df$B1 <= 0] <- NA
df$B2[df$B2 <= 0] <- NA
df$B3[df$B3 <= 0] <- NA
df <- df[complete.cases(df), ]
df

Answer 1

这将批量进行100k个随机参数组合。如果找到结果，它会将其写入CSV。

library(tidyverse)
batch_size <- 100000
batch_num <- 1
repeat {
print(paste0("Processing batch (size = ", batch_size, "): ", batch_num))

search_batch <- 
matrix(
runif(batch_size * 12, min = 0.1, max = 500),
ncol = 12
) %>%
as_tibble() %>%
set_names(c("f1", "ep", "a1", "h", "e1", "m1", "mp", "f", "ab", "r", "fp", "R"))
result <-
search_batch %>%
mutate(
A1 = eval(eqs$A1, envir = .),
A2 = eval(eqs$A2, envir = .),
B1 = eval(eqs$B1, envir = .),
B2 = eval(eqs$B2, envir = .),
B3 = eval(eqs$B3, envir = .)
) %>%
filter(
A1 > 0,
A2 > 0,
B1 > 0,
B2 > 0,
B3 > 0
)

if (nrow(result) > 0) {
print("Found result!")
write_csv(result, "results.csv", append = TRUE)
break
}
batch_num <- batch_num + 1
}