R中矩阵的优化

我对R中的模型的优化/校准是个新手,但我很想学习并真正需要一些帮助。我的问题与人口模型有关。

我做了一些调查,找到了帮助 here 和 here 但他们都没有完全回答我的问题。

我有一个标量矩阵(倾向性),其中每列必须总计为1。这些倾向用于估计特定人口(按年龄划分的人)可能产生的家庭数量。倾向性模型倾向于高估历史上的家庭数量(我知道真实的家庭数量)。

我想通过调整属性来校准模型,以最小化家庭数量中的错误,这样列仍会增加到1,初始值为零的属性必须保持为零。

简单示例:

  # Propensities matrix
  mtx <- matrix(c(0.00,  0.00,  0.85,  0.00,  0.15, 0.35,  0.45,  0.00,  
               0.20,  0.00, 0.65,  0.15,  0.00,  0.20,  0.00), ncol = 3)

  # Population by age cohort
  pop <- c(2600, 16200, 13400)

  # True number of households
  target <- c(7000, 4500, 5500)

  # Function to optimise
  hh <- function(mtx, pop, target) {
    # Estimate living arrangements
    x <- mtx %*% pop 
    # Estimate number of households using parent cohorts (1,2 and 4)
    x <- c(x[1,1]/2, x[2,1]/2, x[4,1]) - target
    return(x)
  }

我没有包含任何优化/校准步骤的代码,因为这会让我很尴尬,而且我也没有任何东西可以工作!

理想情况下,在这个过程的最后,我将有一组对许多不同区域都很好地概括的倾向。关于如何实现它有什么建议吗?有用的链接?

更新

下面的代码片段执行enrico建议的本地搜索方法。

library(tidyverse)
library(NMOF)

data <- list(mtx = matrix(c(0.00,  0.00,  0.90,  0.00,  0.10, 0.25,  0.50,  0.00,   
                            0.25,  0.00, 0.60,  0.20,  0.00,  0.20,  0.00), ncol = 3),
             pop = c(2600, 16200, 13400),
             target = c(7190, 4650, 5920))

# True mtx
mtx.true <- matrix(c(0.00,  0.00,  0.75,  0.00,  0.25, 0.35,  0.45,  0.00,
                     0.20,  0.00, 0.65,  0.15,  0.00,  0.20,  0.00), ncol = 3)

# Function to optimise
households <- function(x, data) {

  # Estimate living arrangements
  z <- x %*% data$pop 

  # Estimate number of households using parent cohorts (1,2 and 4)
  z <- c(z[1,1]/2, z[2,1]/2, z[4,1]) - data$target
  sum(abs(z))

}

# Local search function to perturb propensities
neighbour <- function(x, data) {

  # Choose random column from mtx
  i <- sample(1:ncol(x), 1)
  # Select two non-zero propensities from mtx column
  j <- which(x[, i] != 0) %>% sample(2, replace = FALSE)

  # Randomnly select one to perturb positively 
  x[j[1], i] <- 0.1 * (1 - x[j[1], i]) + x[j[1], i]
  # Perturb second propensity to ensure mtx column adds to 1
  x[j[2], i] <- x[j[2], i] + (1 - sum(x[,i]))

  x

}

# Local search algorithm inputs 
localsearch <- list(x0 = data$mtx,             
                    neighbour = neighbour,
                    nS = 50000,                
                    printBar = FALSE)

# Execute 
now <- Sys.time()
solution <- LSopt(OF = households, algo = localsearch, data)
#> 
#> Local Search.
#> Initial solution:  2695 
#> Finished.
#> Best solution overall: 425.25
Sys.time() - now
#> Time difference of 6.33272 secs

# Inspect propensity matrices
print(solution$xbest)
#>           [,1]   [,2] [,3]
#> [1,] 0.0000000 0.3925  0.6
#> [2,] 0.0000000 0.4250  0.2
#> [3,] 0.2937976 0.0000  0.0
#> [4,] 0.0000000 0.1825  0.2
#> [5,] 0.7062024 0.0000  0.0
print(mtx.true)
#>      [,1] [,2] [,3]
#> [1,] 0.00 0.35 0.65
#> [2,] 0.00 0.45 0.15
#> [3,] 0.75 0.00 0.00
#> [4,] 0.00 0.20 0.20
#> [5,] 0.25 0.00 0.00

谢谢!