我正在尝试在R中编写一段优化代码来计算一组关于果蝇蝇的生物学问题的未知值。
数据帧由13列组成（为清楚起见，下面的代码只显示了9列），行数不同。前三列包含收集的数据，其余列使用各种公式计算。其中两列Missing_C和Missing_D最初填充了空数据，优化问题中的两列代表初始值。
Time.min. Prob_C Prob_D Miss_C Miss_D Event_C Event_D Risk_C Risk_D
1 0 1.00 1.00 0 0 0.00 0.00 86.00 78.00
2 5 0.98 0.97 0 0 1.93 1.98 84.07 76.02
3 16 0.84 0.95 0 0 10.67 1.90 73.40 74.12
4 17 0.50 0.75 0 0 21.02 12.85 52.38 61.27
5 20 0.30 0.50 0 0 14.97 15.32 37.42 45.95
作为使用的一些公式的示例，Event_C和Risk_C使用for循环计算如下：
#define values for events_c and risk_c with for loops`
temp_events_c <-vector()
temp_risk_c <-vector()
for (i in 2:no_rows) {
temp_events_c <- ((prob_c[i] * risk_c[i-1]) - (prob_c[i] * miss_c[i-1]) - (prob_c[i-1] * risk_c[i-1]) + (prob_c[i-1] * miss_c[i-1])) / (prob_c[i] - (2 * prob_c[i-1]))
events_c[i] <- temp_events_c
temp_risk_c <- risk_c[i-1] - miss_c[i-1] - events_c[i]
risk_c[i] <- temp_risk_c
}
根据这些数据，我还有一个收集的值（本例中为9.1），它与表中的值有关。以下函数定义了此值与列Event_C，Event_D和上面未显示的两列（Expected_C和Expected_D）之间的关系，其中这些列的总和由x [1]，x [2]，x [3]表示， X [4]：
fn <- function(x) ((x[1]-x[2])^2)/x[2] + ((x[3]-x[4])^2)/x[4]
然后，我想使用最小化算法，即来自slsqp的nloptr来计算Miss_C和Miss_D中最终满足该单个值的值。优化的额外代码将是这样的：
x0 <- c(Miss_C,Miss_D)
heq <- function(x) (((x[1]-x[2])^2)/x[2] + ((x[3]-x[4])^2)/x[4]) - 9.1 # heq == 0
slsqp(x0, fn, gr = NULL,
hin = NULL, heq = heq)
显然，这不起作用，因为初始值不直接包含在需要解决的函数中，这就是我坚持的观点！我不确定这是否是一个优化问题或更多的一般R编码问题 - 无论哪种方式，任何帮助将非常感激。
干杯，安德鲁
*编辑 - 根据要求添加完整代码*
#input variables
time_vector <- c(0,5,16,17,20)
prob_c <- c(1,0.977,0.835,0.5,0.30)
prob_d <- c(1,0.974,0.949,0.75,0.50)
miss_c <- c(0,0,0,0,0)
miss_d <- c(0,0,0,0,0)
#get number of rows
no_rows <- length(time_vector)
#fill events columns with dummy data
events_c <- c(0:(no_rows - 1))
events_d <- c(0:(no_rows - 1))
#define starting number at risk
risk_c_t0 <- 86
risk_d_t0 <- 78
#add t0 risk to each column
risk_c <- risk_c_t0
risk_d <-risk_d_t0
#fill risk columns with dummy data
risk_c[2:no_rows] <- c(2:no_rows)
risk_d[2:no_rows] <- c(2:no_rows)
#re-define values for events_c and risk_c with for loops
temp_events_c <-vector()
temp_risk_c <-vector()
for (i in 2:no_rows) {
temp_events_c <- ((prob_c[i] * risk_c[i-1]) - (prob_c[i] * miss_c[i-1]) - (prob_c[i-1] * risk_c[i-1]) + (prob_c[i-1] * miss_c[i-1])) / (prob_c[i] - (2 * prob_c[i-1]))
events_c[i] <- temp_events_c
temp_risk_c <- risk_c[i-1] - miss_c[i-1] - events_c[i]
risk_c[i] <- temp_risk_c
}
#re-define values for events_t with for loops
temp_events_d <-vector()
temp_risk_d <-vector()
for (j in 2:no_rows) {
temp_events_d <- ((prob_d[j] * risk_d[j-1]) - (prob_d[j] * miss_d[j-1]) - (prob_d[j-1] * risk_d[j-1]) + (prob_d[j-1] * miss_d[j-1])) / (prob_d[j] - (2 * prob_d[j-1]))
events_d[j] <- temp_events_d
temp_risk_d <- risk_d[j-1] - miss_d[j-1] - events_d[j]
risk_d[j] <- temp_risk_d
}
#calculate total risk, events and expected
total_risk <- risk_c + risk_d
total_events <- events_c + events_d
expected_c <- (risk_c * (total_events/total_risk))
expected_d <- (risk_d * (total_events/total_risk))
#place values into dataframe
df1 <- data.frame(time_vector,prob_c,prob_d, miss_c, miss_d, events_c, events_d, risk_c, risk_d, total_risk, total_events, expected_c, expected_d)
#sum of values
sum_events_C <- sum(events_c)
sum_events_d <- sum(events_d)
sum_expected_c <- sum(expected_c)
sum_expected_d <- sum(expected_d)
#chi_sq formula
chi_sq_combo <- (((sum_events_C - sum_expected_c)^2)/sum_expected_c) + (((sum_events_d - sum_expected_d)^2)/sum_expected_d)
#### end of table calculations before sim
#x <- c(sum_events_C, sum_expected_c, sum_events_d, sum_expected_d)
#x0 <- c(miss_c,miss_d) #inital values
#fn <- function(x) ((x[1]-x[2])^2)/x[2] + ((x[3]-x[4])^2)/x[4]
#heq <- function(x) (((x[1]-x[2])^2)/x[2] + ((x[3]-x[4])^2)/x[4]) - 6.5 # heq == 0
#slsqp(x0, fn, gr = NULL,
# hin = NULL, heq = heq)
重述上面的评论，我认为问题是使用优化来找到两个产生目标卡方值的值。可能导致问题的复杂因素是可能存在许多产生目标的解决方案，因此可能需要添加一些其他要求以使答案唯一。
要做到这一点，你需要一个两个变量的函数来计算使用这些变量的卡方值与目标值之间的差值的平方，然后你最小化它。
例如，
fn2 <- function(x) {
c <- x[1]
d <- x[2]
chisq <- (((c - sum_expected_c)^2)/sum_expected_c) +
(((d - sum_expected_d)^2)/sum_expected_d)
(chisq - 6.5)^2
}
for (i in 1:no_rows) {
x0 <- c(df1$miss_c[i],df1$miss_d[i]) #initial values
res <- nloptr::slsqp(x0, fn2)
miss_c[i] <- res$par[1]
miss_d[i] <- res$par[2]
}
这给了所有5次相同的值，所以我可能没有完全理解你。