Оптимизация матрицы в R

Question

Я новичок в оптимизации / калибровке моделей в R, но я хочу учиться и действительно нуждаюсь в некоторой помощи. Мой вопрос касается демографического моделирования.

Я провел некоторые исследования и нашел здесь и здесь помощь, но ни один из них не ответил на мой вопрос.

У меня есть матрица скаляров (склонностей), где каждый столбец должен составлять до 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)
  }

Я не включил свой код для этапа оптимизации / калибровки, так как это было бы неловко, и я не смог заставить что-либо работать!

В идеале, в конце этого процесса у меня будет один набор склонностей, который будет хорошо обобщаться для множества разных регионов. Любой совет о том, как я должен идти о достижении этого? Полезные ссылки?

Обновить

Приведенный ниже фрагмент кода выполняет метод локального поиска, предложенный Энрико.

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

Спасибо!

Answer 1

Я могу только прокомментировать часть оптимизации.

Код, который вы предоставили, достаточно; только ваша целевая функция оценивает вектор. Вам нужно будет преобразовать этот вектор в одно число, которое нужно минимизировать, например, сумму квадратов или абсолютных значений.

Когда дело доходит до методов, я бы попробовал эвристику; на самом деле, я бы попробовал метод локального поиска. Эти методы работают с решением через функции, которые вы определяете; таким образом, вы можете закодировать свое решение в виде матрицы. Более конкретно, вам понадобятся две функции: целевая функция (которая у вас по сути есть) и функция соседства, которая принимает в качестве входных данных решение и модифицирует его. В вашем конкретном случае это может занять матрицу, выбрать два ненулевых элемента из одного столбца, увеличить один и уменьшить другой. Таким образом, сумма столбца останется неизменной.

Возможно, учебник http://enricoschumann.net/files/NMOF_Rmetrics2012.pdf представляет интерес, с кодом R http://enricoschumann.net/files/NMOF_Rmetrics2012.R.