У меня проблемы с оценкой модели полиномиального логита, мой код предоставляет коэффициенты, отличные от коэффициентов, оцененных с помощью "mlogit".
Мой код:
ll.mnl <- function(beta,x,y){
b1 <- beta[1:3]
b2 <- beta[4:6]
x1 <- x[,1:3]
x2 <- x[,4:6]
u1 <- exp(x1%*%b1)
u2 <- exp(x2%*%b2)
u_total <- 1+u1+u2
p1 <- u1/u_total
p2 <- u2/u_total
p3 <- 1-p1-p2
ll <- prod(p1^y[,1])*prod(p2^y[,2])*prod(p3^y[,3])
log_ll <- -log(ll)
return(log_ll)
}
x_auto <- cbind(1,basep$t.auto,basep$c.auto)
x_bus <- cbind(1,basep$t.bus,basep$c.bus)
x_rail <- cbind(1,basep$t.rail,basep$c.rail)
x <- cbind(x_auto,x_bus)
y <- cbind(basep$y.auto,basep$y.bus,basep$y.rail)
bb<-c(2,5,0.1,1,-3,-5)
optim_ll.mnl<-
optim(bb,x=x,y=y,fn=ll.mnl,method="BFGS",control=list(trace=FALSE,
fnscale=1,maxit=10000))
optim_ll.mnl
и результаты:
optim_ll.mnl
par
[1] 0.0198305210 0.0054421463 -0.0118098673 -0.4719502006
-0.0007723349 0.0121851407
$value
[1] 63.12357
$counts
function gradient
120 15
$convergence
[1] 0
$message
NULL
код, используемый в "mlogit":
H3 <- mlogit.data(heat, shape = "wide", choice = "depvar", varying = c(3:8))
m6 <- mlogit(depvar ~ 0 | oc+ic,H3,reflevel = "er",method="bfgs")
summary(m6)
и результаты "mlogit":
Call:
mlogit(formula = depvar ~ 0 | oc + ic, data = H3, reflevel = "er",
method = "bfgs", print.level = 0)
Frequencies of alternatives:
er gc gr
0.100 0.725 0.175
bfgs method
8 iterations, 0h:0m:0s
g'(-H)^-1g = 4.54E-07
gradient close to zero
Coefficients :
Estimate Std. Error z-value Pr(>|z|)
gc:(intercept) 0.5786927 2.3728927 0.2439 0.80733
gr:(intercept) -1.7413422 2.6891081 -0.6476 0.51727
gc:oc 0.0031484 0.0141568 0.2224 0.82401
gr:oc -0.0366012 0.0288349 -1.2693 0.20432
gc:ic 0.0011817 0.0038685 0.3055 0.76001
gr:ic 0.0085216 0.0040121 2.1240 0.03367 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Log-Likelihood: -58.473
McFadden R^2: 0.04882
Likelihood ratio test : chisq = 6.0023 (p.value = 0.19898)</raw>
my database is: (first four individuals)
<raw>id depvar ic.gc ic.gr ic.er oc.gc oc.gr oc.er income age rooms region y.gc y.gr y.er
1 gc 866.00 962.64 995.76 199.69 151.72 505.60 7 25 6 ncostl 1 0 0
2 gc 727.93 758.89 894.69 168.66 168.66 486.49 5 60 5 scostl 1 0 0
3 gc 599.48 783.05 900.11 165.58 137.80 404.74 4 65 2 ncostl 1 0 0
4 er 835.17 793.06 831.04 180.88 147.14 425.22 2 50 4 scostl 0 0 1
Я должен создать код вручную, мне нужно оценивать коэффициенты по переменным и по режимам в дополнение к их перехватам, поэтому я использую вашу помощь. Заранее большое спасибо.