Автоматизация функции для возврата выражения с математическими константами и неизвестными

Я пытаюсь построить матрицу переходов из наблюдений за данными панели, чтобы получить оценки ML взвешенной матрицы переходов. Ключевым шагом является получение индивидуальной функции правдоподобия для отдельных лиц. Скажем, у вас есть следующий фрейм данных:

ID          Feature1  Feature2  Transition
120421006   10000        1         ab
120421006   12000        0         ba
120421006   10000        1         ab
123884392    3000        1         ab
123884392    2000        0         ba
908747738    1000        1         ab

Идея состоит в том, чтобы вернуть для каждого агента логарифмическую вероятность его пути. Например, для агента 120421006 это сводится к (игнорируя начальный термин)

LL = log(exp(Yab)/1 + exp(Yab)) + log(exp(Yba)/(1 + exp(Yba))) + log(exp(Yab)/1 + exp(Yab))

i.e,

журнал (exp (Y_transition)/(1 + exp (Y_transition)))

где Y_transition = xFeature1 + yFeature2 для этого перехода, а x и y неизвестны.

Например, для индивидуума 120421006 это сведется к выражению с тремя элементами, так как он делает переход трижды, и функция вернет

LL = log(exp(10000x + 1y)/ 1 + exp(10000x + 1y)) +

log(exp(12000x + 0y)/ 1 + exp(12000x + 0y)) +

log(exp(10000x + 1y)/ 1 + exp(10000x + 1y))

И вот загвоздка: мне нужно, чтобы x и y возвращались как неизвестные, поскольку цель состоит в том, чтобы получить сумму вероятностей всех людей, чтобы передать ее оценщику ML. Как бы вы автоматизировали функцию, которая возвращает этот вывод для всех идентификаторов?

Спасибо заранее


r log-likelihood
person Arrebimbomalho    schedule 09.02.2018    source источник

Ответы (2)


arrow_upward
1
arrow_downward

Сначала вы должны решить, насколько гибкой должна быть ваша функция. Я оставляю его довольно жестким, но вы можете изменить его по своему вкусу.

Во-первых, вы должны ввести начальные параметры предположения, которые вы предоставите оптимизатору. Затем объявите свои данные и переменные, которые будут использоваться в вашей оценке.

Предполагая, что у вас всегда будет только 2 переменные (вы можете изменить их позже)

y <- function(initial_param, data, features){

  x = initial_param[1]
  y = initial_param[2]

  F1 = data[, features[1]]
  F2 = data[, features[2]]

  LL = log(exp(F1 * x + F2 * y) / (1 + exp(F1 * x + F2 * y)))

  return(-sum(LL))
}

Эта функция возвращает сумму минус логарифмическая вероятность, учитывая, что большинство оптимизаторов пытаются найти параметры, при которых ваша функция достигает минимума по умолчанию.

Чтобы найти ваши параметры, просто введите в приведенную ниже функцию функцию правдоподобия y, начальные параметры, набор данных и вектор с именами ваших переменных:

nlm(f = y,  initial_param = your_starting_guess, data = your_data,
                  features = c("name_of_first_feature", "name_of_second_feature"), iterlim=1000, hessian=F)
person Felipe Alvarenga    schedule 09.02.2018
comment
Спасибо, Фелипе! Мне просто нужна логарифмическая вероятность, а не решение проблемы, так как реальный набор данных имеет огромное количество переменных, индивидуумов и переходов, и, возможно, было бы лучше перенести это на C и использовать правильный оптимизатор... если я попросите return(LL) в последней части вашего кода, он вернет список логарифмических правдоподобий для агентов, верно? (Я не могу проверить код прямо сейчас :)) Тогда, если бы я мог получить сумму по элементам списка, я был бы установлен... - person Arrebimbomalho; 09.02.2018
comment
sum(LL) вернет список вероятностей для каждого наблюдения. Затем вы должны свернуть этот результат на уровень агента, если хотите получить сумму вероятности для каждого агента. - person Felipe Alvarenga; 09.02.2018
comment
Ok. Но я имею в виду, что если есть какой-то способ получить просто выражение вероятности для каждого агента без фактического решения... вывод должен возвращать LL = log(exp(F1 * x + F2 * y) / (1 + exp(F1 * x + F2 * y))) со значениями F1 и F2 для каждого агента, но x и y как неизвестные... шаг решения отложен на потом, когда этот шаг масштабируется :) - person Arrebimbomalho; 09.02.2018

arrow_upward
1
arrow_downward

Создайте функцию:

fun=function(x){
a=paste0("exp(",x[1],"*x","+",x[2],"*y)")
parse(text=paste("sum(",paste0("log(",a,"/(1+",a,"))"),")"))
}

by(test[2:3],test[,1],fun)

sum(log(exp(c(10000, 12000, 10000) * x + c(1, 0, 1) * y)/(1 + 
    exp(c(10000, 12000, 10000) * x + c(1, 0, 1) * y))))
-------------------------------------------------------------------- 
sum(log(exp(c(3000, 2000) * x + c(1, 0) * y)/(1 + exp(c(3000, 
    2000) * x + c(1, 0) * y))))
-------------------------------------------------------------------- 
sum(log(exp(1000 * x + 1 * y)/(1 + exp(1000 * x + 1 * y))))

на примере x=0 и y=3 мы можем решить это:

x=0
y=3
sapply(by(test[2:3],test[,1],fun),eval)
[1] -0.79032188 -0.74173453 -0.04858735

в вашем примере выше:

x=0
y=3
 log(exp(10000*x + 1*y)/ (1 + exp(10000*x + 1*y))) +#There should be paranthesis
  log(exp(12000*x + 0*y)/ (1 + exp(12000*x + 0*y))) + 
  log(exp(10000*x + 1*y)/( 1 + exp(10000*x + 1*y)))
[1] -0.7903219

Чтобы получить то, что вам нужно в комментариях:

fun1=function(x){
    a=paste0("exp(",x[1],"*x","+",x[2],"*y)")
    paste("sum(",paste0("log(",a,"/(1+",a,"))"),")")
    }

paste(by(test[2:3],test[,1],fun1),collapse = "+")
1] "sum( log(exp(c(10000, 12000, 10000)*x+c(1, 0, 1)*y)/(1+exp(c(10000, 12000, 10000)*x+c(1, 0, 1)*y))) )+sum( log(exp(c(3000, 2000)*x+c(1, 0)*y)/(1+exp(c(3000, 2000)*x+c(1, 0)*y))) )+sum( log(exp(1000*x+1*y)/(1+exp(1000*x+1*y))) )"

Но это не имеет смысла, почему вы группируете их, а затем суммируете их все. Это то же самое, что просто суммировать их без группировки по идентификатору, что было бы проще и быстрее.

person Onyambu    schedule 09.02.2018
comment
Супер! Я, вероятно, экспортирую функцию, состоящую из сумм вероятностей, либо в Java, либо в C++, поскольку реальный набор данных душит R, хотя попробовать функцию optim.R, которая имеет некоторые итерационные возможности и алгоритмы, вероятно, является хорошей идеей... Кстати, знаешь ли ты, как получить единственную функцию, состоящую из суммы вероятностей с константами, из твоего замечательного кода - person Arrebimbomalho; 09.02.2018
comment
Извините, я не понимаю, что вы имеете в виду - person Onyambu; 09.02.2018
comment
Что ж, ваш код в первом прямоугольнике возвращает выражение, которое я искал (вероятность агентов с x и y, оставленных в качестве переменных), и мне не было интересно, есть ли способ собрать все это в большое выражение, например. ................................................. .............сумма(log(exp(c(10000, 12000, 10000) * x + c(1, 0, 1) * y)/(1 + exp(c(10000) , 12000, 10000) * x + c(1, 0, 1) * y)))) сумма(log(exp(c(3000, 2000) * x + c(1, 0) * y)/(1 + exp(c(3000, 2000) * x + c(1, 0) * y)))) , т. е. суммированные вероятности агентов, но со свободными переменными x и y. - person Arrebimbomalho; 09.02.2018
comment
так вы хотите суммировать все вероятности? - person Onyambu; 09.02.2018
comment
Что-то вроде. Я хочу получить выражение, состоящее из суммы отдельных вероятностей, но с еще свободными и у. Таким образом, если есть способ автоматического получения выражения sum(log(exp(c(10000, 12000, 10000) * x + c(1, 0, 1) * y)/(1 + exp(c(10000, 12000, 10000) * x + c(1, 0, 1) * y)))) + sum(log(exp(c(3000, 2000)) * x + c(1, 0) * y)/(1 + exp( c(3000, 2000) * x + c(1, 0) * y)))) это было бы здорово (это действительно объединяет вывод вашей функции выше) - person Arrebimbomalho; 09.02.2018
comment
Идея состоит в том, чтобы затем передать функцию правильному решателю, используя, скажем, BFGS или сопряженный градиент на языке, который может обрабатывать большие данные и большие входные данные... - person Arrebimbomalho; 09.02.2018
comment
Я могу это сделать! Только не разбирайте уравнения в конце: - person Onyambu; 09.02.2018
comment
в функции удалите часть parse(text= и верните только вставку. Затем при выполнении просто сверните их с помощью знака + - person Onyambu; 09.02.2018