使用 optim 找到最小化，同时也强制参数总和为 1

Question

我正在尝试使用 R 和 optim 来计算混合比例。因此，例如，假设我有一个岩石成分，60% SiO2 和 40% CaO。我想知道我必须混合多少两个不同的阶段才能生产出我的岩石。假设第 2 阶段是 35% SiO2 和 65% CaO，第 2 阶段是 80% SiO2 和 20% CaO。

*** 编辑：我已更新代码以包含第三阶段，并尝试使用组合包。我还尝试设置优化搜索范围的界限。

#Telling R the composition of both phases and the target rock

library(compositions)

Phase1 <- c(35, 65)
Phase2 <- c(80, 20)
Phase3 <- c(10, 90)
Target_Composition <- c(60, 40)

#My function to minimize
my.function <- function(n){
    n <- clo(n) #I though this would make 1 = n[1] + n[2] + n[3]
    (Target_Composition[1] - (n[1]*Phase1[1] + n[2]*Phase2[1] + n[3]*Phase3[1]))^2 + 
    (Target_Composition[2] - (n[1]*Phase1[2] + n[2]*Phase2[2] + n[3]*Phase3[2]))^2 
}

然后我在优化中运行它：

optim(c(.54,.46), my.function)

当我运行它时，我得到 0.536 和 0.487，这确实是最小值，但是，我需要设置一个额外的参数 n[1] + n[2] = 1 。有没有办法使用优化来做到这一点？这就是让我困惑的地方。作为旁注，我想要解决的实际问题有更多的阶段，每个阶段都更复杂，但是，一旦我让这部分工作，我就会扩大规模。

我现在使用：

optim(c(.5, .4, .1), my.function, lower=0, upper=1, method="L-BFGS-B")

我现在得到 0.4434595、0.4986680 和 0.2371141 作为结果，它们的总和不为 1。我该如何解决这个问题？另外，我在 0 和 1 之间设置搜索范围的新方法是否有效？

谢谢你的帮助。

Answer 1

对于一个参数估计 -optimize()函数用于最小化函数。

\n\n

对于两个或多个参数估计，optim()函数用于最小化函数。基 R 中还有另一个函数constrOptim()，可用于执行不等式约束的参数估计。在您的问题中，您打算应用框约束。该L-BFGS-B求解器是合适的。

\n\n

注意： 并非所有求解器都会收敛，但它们都会根据求解器中实现的某些条件而终止（您可以使用control参数控制它们），因此您必须检查从求解器获得的最小值是否通过 KKT（Karush、Kuhn 和Tucker）条件，以确保您的求解器实际上已收敛，同时最佳地估计参数。

\n\n

KKT条件

\n\n

\n
1. 对于最小值，梯度（一阶导数）必须为 \xe2\x80\x9czero\xe2\x80\x9d
\n
2. 对于最小值，Hessian（二阶偏导数）必须是正定对称矩阵。
\n

\n\n

您还可以使用 hessian 矩阵的行列式来查找最小值是局部的、全局的还是鞍点。

\n\n

在这里，我展示了为给定表达式创建目标函数的两种方法，例如符号代数形式和矩阵形式。我还比较了这两个函数的输出，但是，我使用矩阵形式的目标函数来显示优化步骤。请注意，优化步骤满足约束条件您在问题中提到的

\n\n

n1 + n2 + n3 = 1\nn1 = (-2, 2)\nn2 = (-1, 1)\nn3 = (-1, 1)\n

\n\n

optimx您可能想尝试使用使工作变得更轻松的库。

\n\n

# load library\nlibrary(\'compositions\')\n\n## function for getting algebraic expressions \nget_fexpr <- function( phase_len, n_ingredients ) # len = length of n or phase1 or phase2 or target composition\n{\n  ## phase_len = number of phases (eg: phase1, phase2, phase3: = 3)\n  ## n_ingredients = number of ingredients that form a composition (Eg: sio2 and cao: = 2)\n  ## y = target composition\n  ## x = phase data as c(phase1[1], phase2[1], phase1[2], phase2[2])\n  ## n = parameters to be estimated\n\n  p_n <- paste( rep("n", phase_len), 1:phase_len, sep = \'\')  # n\n  p_x <- paste( rep("x", phase_len), paste( "[", 1:( phase_len * n_ingredients ), "]", sep = \'\'), sep = \'\' ) # x\n  p_y <- paste( "y[", 1:n_ingredients, "]", sep = ""  )    # y\n\n  # combine n, x, and y\n  p_n_x <- paste( "(", paste( p_n, p_x, sep = "*" ), ")", sep = \'\')\n  p_n_x <- lapply( split( p_n_x, ceiling( seq_along( p_n_x ) / phase_len ) ), paste, collapse = " + " )\n  p_n_x <- lapply( p_n_x, function( x ) paste( "(", x, ")", sep = "" ) )\n\n  # get deviations and sum of squares\n  dev    <- mapply( paste, p_y, p_n_x, sep = " - " )  # deviations\n  dev_sq <- paste( "(", dev, ")**2", sep = \'\', collapse = " + ")  # sum of squares of deviations\n\n  return( dev_sq )\n}\n\nget_fexpr( phase_len = 1, n_ingredients = 1 )\n# [1] "(y[1] - ((n1*x[1])))**2"\nget_fexpr( phase_len = 1, n_ingredients = 2 )\n# [1] "(y[1] - ((n1*x[1])))**2 + (y[2] - ((n1*x[2])))**2"\nget_fexpr( phase_len = 2, n_ingredients = 2 )\n# [1] "(y[1] - ((n1*x[1]) + (n2*x[2])))**2 + (y[2] - ((n1*x[3]) + (n2*x[4])))**2"\nget_fexpr( phase_len = 3, n_ingredients = 2 )\n# [1] "(y[1] - ((n1*x[1]) + (n2*x[2]) + (n3*x[3])))**2 + (y[2] - ((n1*x[4]) + (n2*x[5]) + (n3*x[6])))**2"\nget_fexpr( phase_len = 4, n_ingredients = 2 )\n# [1] "(y[1] - ((n1*x[1]) + (n2*x[2]) + (n3*x[3]) + (n4*x[4])))**2 + (y[2] - ((n1*x[5]) + (n2*x[6]) + (n3*x[7]) + (n4*x[8])))**2"\n\n## objective functions for max/min\n## matrix form\nmyfun1 <- function( n, phase_data, Target_Composition )\n{\n  print(n)\n  Target_Composition <- matrix(Target_Composition, ncol = length(Target_Composition), byrow = TRUE)\n  dot_product <- n %*% phase_data   # get dot product\n  sum((Target_Composition - dot_product)^2)\n}\n\n## algebraic expression form\nmyfun2 <- function( y, x, n, fexpr )\n{\n  names(n) <- paste( "n", 1:length( n ), sep = "" ) # assign names to n\n  list2env( as.list(n), envir = environment() )     # create new variables with a list of n\n  return( eval( parse( text = fexpr ) ) )\n} \n\n## Comparison of functions of matrix and algebriac forms\n## 1. raw data for two parameters estimation\nTarget_Composition <- clo( c(60, 40), total = 1 )\nPhase1 <- clo( c(35, 65), total = 1 )\nPhase2 <- clo( c(80, 20), total = 1 )\nstopifnot( length( Phase1 ) == length( Phase2 ))  \nn      <- clo( c(0.54, 0.46) , total = 1 )\n\n## data for matrix form\nphase_data_concat <- c(Phase1, Phase2)\nn_row <- length(phase_data_concat) / length(Phase1)\nn_col <- length(phase_data_concat) / n_row\nphase_data <- matrix( data = phase_data_concat, \n                      nrow = n_row, \n                      ncol = n_col, \n                      byrow = TRUE,\n                      dimnames = list(c(\'phase1\', \'phase2\'), \n                                      c(\'sio2\', \'cao\')))\n\n## data for algebraic form\ny <- Target_Composition\nx <- c(matrix( data = phase_data_concat, nrow = nrow( phase_data ), ncol = ncol( phase_data ), byrow = TRUE ))\nn <- n\nfexpr <- get_fexpr( phase_len = length( n ), n_ingredients = 2 )\n\n\n# compare algebraic form and matrix form functions\nmyfun1(n, phase_data, Target_Composition)\n# [1] 0.0036979999999999969\nmyfun2( y = y, x = x, n = n, fexpr = fexpr )\n# [1] 0.0036979999999999969\n\n## 2. raw data for three parameters estimation\nTarget_Composition <- clo( c(60, 40), total = 1)\nPhase1 <- clo( c(35, 65), total = 1)\nPhase2 <- clo( c(80, 20), total = 1)\nPhase3 <- clo( c(10, 90), total = 1)\nstopifnot( length( Phase1 ) == length( Phase2 ) && length( Phase1 ) == length( Phase3 ))  \nn <- clo( c(0.5, 0.4, 0.1) )   # starting guess for n1, n2, and n3\n\n## data for matrix form\nphase_data_concat <- c(Phase1, Phase2, Phase3)\nn_row <- length(phase_data_concat) / length(Phase1)\nn_col <- length(phase_data_concat) / n_row\nphase_data <- matrix( data = phase_data_concat, \n                      nrow = n_row, \n                      ncol = n_col, \n                      byrow = TRUE,\n                      dimnames = list(c(\'phase1\', \'phase2\', \'phase3\'), \n                                      c(\'sio2\', \'cao\')))\n## data for algebraic form\ny <- Target_Composition\nx <- c( matrix( phase_data_concat, nrow = nrow( phase_data), ncol = ncol(phase_data), byrow = TRUE ) )\nn <- n\nfexpr <- get_fexpr( phase_len = length( n ), n_ingredients = 2 )\n\n# compare algebraic form and matrix form functions\nmyfun1(n, phase_data, Target_Composition)\n# [1] 0.01805\nmyfun2( y = y, x = x, n = n, fexpr = fexpr )\n# [1] 0.01805\n\n\n## Optimization using matrix form objective function (myfun1)\n# three parameter estimation\nphase_data\n#                       sio2                 cao\n# phase1 0.34999999999999998 0.65000000000000002\n# phase2 0.80000000000000004 0.20000000000000001\n# phase3 0.10000000000000001 0.90000000000000002\n\n# target data\nTarget_Composition <- clo( c(60, 40), total = 1 )\n# [1] 0.59999999999999998 0.40000000000000002\n\nn <- c(0.5, 0.4, 0.1)\n# [1] 0.50000000000000000 0.40000000000000002 0.10000000000000001\n\n## Harker diagram; also called scatterplot of two componenets without any transformation\nplot( phase_data, type = "b", main = "Harker Diagram" )   \n\noptim_model <- optim(par = n, \n                     fn = myfun1,\n                     method = "L-BFGS-B", \n                     lower = c(-2, -1, -1),  # lower bounds: n1 = -2; n2 = -1; n3 = -1 \n                     upper = c( 2, 1, 1 ),   # upper bounds: n1 = 2;  n2 = 1   \n                     phase_data = phase_data,\n                     Target_Composition = Target_Composition)\n\n# [1]  0.50000000000000000   0.40000000000000002   0.10000000000000001\n# [1]  0.50100000000000000   0.40000000000000002   0.10000000000000001\n# [1]  0.49900000000000000   0.40000000000000002   0.10000000000000001\n# [1]  0.50000000000000000   0.40100000000000002   0.10000000000000001\n# [1]  0.50000000000000000   0.39900000000000002   0.10000000000000001\n# [1]  0.50000000000000000   0.40000000000000002   0.10100000000000001\n# [1]  0.500000000000000000  0.400000000000000022  0.099000000000000005\n# [1]  0.443000000000007166  0.514000000000010004 -0.051999999999988944\n# [1]  0.444000000000007167  0.514000000000010004 -0.051999999999988944\n# [1]  0.442000000000007165  0.514000000000010004 -0.051999999999988944\n# [1]  0.443000000000007166  0.515000000000010005 -0.051999999999988944\n# [1]  0.443000000000007166  0.513000000000010004 -0.051999999999988944\n# [1]  0.443000000000007166  0.514000000000010004 -0.050999999999988943\n# [1]  0.443000000000007166  0.514000000000010004 -0.052999999999988945\n# [1]  0.479721654922847740  0.560497432130581008 -0.020709332414779191\n# [1]  0.480721654922847741  0.560497432130581008 -0.020709332414779191\n# [1]  0.478721654922847739  0.560497432130581008 -0.020709332414779191\n# [1]  0.479721654922847740  0.561497432130581009 -0.020709332414779191\n# [1]  0.479721654922847740  0.559497432130581007 -0.020709332414779191\n# [1]  0.47972165492284774   0.56049743213058101  -0.01970933241477919\n# [1]  0.479721654922847740  0.560497432130581008 -0.021709332414779191\n# [1]  0.474768384608834304  0.545177697419992557 -0.019903455841806163\n# [1]  0.475768384608834305  0.545177697419992557 -0.019903455841806163\n# [1]  0.473768384608834303  0.545177697419992557 -0.019903455841806163\n# [1]  0.474768384608834304  0.546177697419992558 -0.019903455841806163\n# [1]  0.474768384608834304  0.544177697419992557 -0.019903455841806163\n# [1]  0.474768384608834304  0.545177697419992557 -0.018903455841806163\n# [1]  0.474768384608834304  0.545177697419992557 -0.020903455841806164\n# [1]  0.474833910147636595  0.544703104470840138 -0.019537864476362268\n# [1]  0.475833910147636596  0.544703104470840138 -0.019537864476362268\n# [1]  0.473833910147636594  0.544703104470840138 -0.019537864476362268\n# [1]  0.474833910147636595  0.545703104470840139 -0.019537864476362268\n# [1]  0.474833910147636595  0.543703104470840137 -0.019537864476362268\n# [1]  0.474833910147636595  0.544703104470840138 -0.018537864476362267\n# [1]  0.474833910147636595  0.544703104470840138 -0.020537864476362269\n# [1]  0.474834452107517901  0.544702005703077585 -0.019536411001123268\n# [1]  0.475834452107517902  0.544702005703077585 -0.019536411001123268\n# [1]  0.473834452107517901  0.544702005703077585 -0.019536411001123268\n# [1]  0.474834452107517901  0.545702005703077586 -0.019536411001123268\n# [1]  0.474834452107517901  0.543702005703077584 -0.019536411001123268\n# [1]  0.474834452107517901  0.544702005703077585 -0.018536411001123267\n# [1]  0.474834452107517901  0.544702005703077585 -0.020536411001123269\n\noptim_model$par   # values of n after minimization of function my.function using starting guess of n1 = 0.54, n2 = 0.46\n# [1]  0.474834452107517901  0.544702005703077585 -0.019536411001123268\n\nsum(optim_model$par)\n# [1] 1.0000000468094723\n\n# different starting guess - n\nn <- c(0.2, 0.2, 0.6)\n\noptim_model <- optim(par = n, \n                     fn = myfun1,\n                     method = "L-BFGS-B", \n                     lower = c(-2, -1, -1),  # lower bounds: n1 = -2; n2 = -1; n3 = -1 \n                     upper = c( 2, 1, 1 ),   # upper bounds: n1 = 2;  n2 = 1   \n                     phase_data = phase_data,\n                     Target_Composition = Target_Composition)\n\n# [1] 0.20000000000000001 0.20000000000000001 0.59999999999999998\n# [1] 0.20100000000000001 0.20000000000000001 0.59999999999999998\n# [1] 0.19900000000000001 0.20000000000000001 0.59999999999999998\n# [1] 0.20000000000000001 0.20100000000000001 0.59999999999999998\n# [1] 0.20000000000000001 0.19900000000000001 0.59999999999999998\n# [1] 0.20000000000000001 0.20000000000000001 0.60099999999999998\n# [1] 0.20000000000000001 0.20000000000000001 0.59899999999999998\n# [1] 0.014000000000008284 0.571999999999969644 0.103999999999989656\n# [1] 0.015000000000008284 0.571999999999969644 0.103999999999989656\n# [1] 0.013000000000008283 0.571999999999969644 0.103999999999989656\n# [1] 0.014000000000008284 0.572999999999969645 0.103999999999989656\n# [1] 0.014000000000008284 0.570999999999969643 0.103999999999989656\n# [1] 0.014000000000008284 0.571999999999969644 0.104999999999989657\n# [1] 0.014000000000008284 0.571999999999969644 0.102999999999989655\n# [1] 0.13382855816928069 0.72372846274181657 0.20610638896226857\n# [1] 0.13482855816928069 0.72372846274181657 0.20610638896226857\n# [1] 0.13282855816928069 0.72372846274181657 0.20610638896226857\n# [1] 0.13382855816928069 0.72472846274181657 0.20610638896226857\n# [1] 0.13382855816928069 0.72272846274181657 0.20610638896226857\n# [1] 0.13382855816928069 0.72372846274181657 0.20710638896226857\n# [1] 0.13382855816928069 0.72372846274181657 0.20510638896226857\n# [1] 0.11766525503937687 0.67373774947575715 0.20873609146356592\n# [1] 0.11866525503937687 0.67373774947575715 0.20873609146356592\n# [1] 0.11666525503937687 0.67373774947575715 0.20873609146356592\n# [1] 0.11766525503937687 0.67473774947575715 0.20873609146356592\n# [1] 0.11766525503937687 0.67273774947575715 0.20873609146356592\n# [1] 0.11766525503937687 0.67373774947575715 0.20973609146356592\n# [1] 0.11766525503937687 0.67373774947575715 0.20773609146356592\n# [1] 0.11787907521862623 0.67218907774694359 0.20992907381396153\n# [1] 0.11887907521862623 0.67218907774694359 0.20992907381396153\n# [1] 0.11687907521862623 0.67218907774694359 0.20992907381396153\n# [1] 0.11787907521862623 0.67318907774694359 0.20992907381396153\n# [1] 0.11787907521862623 0.67118907774694359 0.20992907381396153\n# [1] 0.11787907521862623 0.67218907774694359 0.21092907381396153\n# [1] 0.11787907521862623 0.67218907774694359 0.20892907381396153\n# [1] 0.11788084371929158 0.67218549229424496 0.20993381673316230\n# [1] 0.11888084371929158 0.67218549229424496 0.20993381673316230\n# [1] 0.11688084371929158 0.67218549229424496 0.20993381673316230\n# [1] 0.11788084371929158 0.67318549229424496 0.20993381673316230\n# [1] 0.11788084371929158 0.67118549229424496 0.20993381673316230\n# [1] 0.11788084371929158 0.67218549229424496 0.21093381673316230\n# [1] 0.11788084371929158 0.67218549229424496 0.20893381673316230\n\noptim_model$par   # values of n after minimization of function my.function using starting guess of n1 = 0.54, n2 = 0.46\n# [1] 0.11788084371929158 0.67218549229424496 0.20993381673316230\n\nsum(optim_model$par)\n# [1] 1.0000001527466988\n

\n\n

可视化： \n我修改myfun1为myfun3使用框约束可视化目标函数。我使用了两个参数估计模型。

\n\n

# modified function for visualization\nmyfun3 <- function( n1, n2, phase_data, Target_Composition )\n{\n  Target_Composition <- matrix(Target_Composition, ncol = length(Target_Composition), byrow = TRUE)\n  dot_product <- c(n1, n2) %*% phase_data   # get dot product\n  sum((Target_Composition - dot_product)^2)\n}\nmyfun3 <- Vectorize(FUN = myfun3, vectorize.args = list( "n1", "n2"))\n\nTarget_Composition <- clo( c(60, 40), total = 1 )\nPhase1 <- clo( c(35, 65), total = 1 )\nPhase2 <- clo( c(80, 20), total = 1 )\nstopifnot( length( Phase1 ) == length( Phase2 ))  \nn      <- clo( c(0.54, 0.46) , total = 1 )\n\n## data for matrix form\nphase_data_concat <- c(Phase1, Phase2)\nn_row <- length(phase_data_concat) / length(Phase1)\nn_col <- length(phase_data_concat) / n_row\nphase_data <- matrix( data = phase_data_concat, \n                      nrow = n_row, \n                      ncol = n_col, \n                      byrow = TRUE,\n                      dimnames = list(c(\'phase1\', \'phase2\'), \n                                      c(\'sio2\', \'cao\')))\n\n\nn1 = seq(-2, 2, 0.1)  # bounds for n1\nn2 = seq(-1, 1, 0.1)  # bounds for n2\nz <- outer( n1, n2, phase_data, Target_Composition, FUN = myfun3)\npersp(n1, n2, z,theta=150)\n

\n\n

\n