0

我正在尝试从面板数据观察中构建一个转换矩阵，以获得加权转换矩阵的 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))

IE，

日志（exp（Y_transition）/（1 + exp（Y_transition）））

其中 Y_transition = xFeature1 + yFeature2 表示该转换，x 和 y 是未知数。

例如，对于个人 120421006，这将归结为具有三个元素的表达式，因为他转换了三次，并且函数将返回

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

日志（exp（12000x + 0y）/ 1 + exp（12000x + 0y））+

对数（exp（10000x + 1y）/ 1 + exp（10000x + 1y））

这里有一个问题：我需要 x 和 y 作为未知数返回，因为目标是获得所有个体可能性的总和，以便将其传递给 ML 估计器。您将如何自动化为所有 ID 返回此输出的函数？

提前谢谢了

4

2 回答 2

1

首先，您必须决定您的功能必须有多灵活。我让它相当僵硬，但你可以根据自己的口味改变它。

首先，您必须输入将在优化器中提供的初始猜测参数。然后，声明要在估计中使用的数据和变量。

假设您总是只有 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)
于 2018-02-09T18:27:38.737 回答
1

创建函数：

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))) )"

但这没有意义，为什么您要将它们分组然后将它们全部相加。这与仅将它们相加而不使用 ID 对它们进行分组相同，这将更简单、更快

于 2018-02-09T18:18:11.873 回答