I'm writing a code to have a travel demand forecasting using the Multinomial Logit Model without using packages. However, I'm having an error saying: Error in fn(par, ...) :
dims [product 121] do not match the length of object [0]
I have in 121 numbers of sample size.
I'd appreciate if you take a look at my code.
My data is shown below:
FModeFr Age Drlicense Work Income PosAuto CostCar CostBus CostTaxi CostSharedtaxi Distance TimeCar TimeBus TimeTaxi TimeSharedtaxi TimeWalk
1 1 2 1 3 3 1 1000 500 5000 1000 5 30 30 30 30 60
2 5 2 1 5 2 4 1400 500 2000 1000 2 15 12 15 15 24
My code is shown below:
hh <-nrow(data)
fr <-function(x) {
Car <-x[1]
Bus <-x[2]
Taxi <-x[3]
Sharedtaxi <-x[4]
Age <-x[5]
Drlicense <-x[6]
Work <-x[7]
Income <-x[8]
Cost <-x[9]
Distance <-x[10]
Time <-x[11]
PosAuto <-x[12]
LL=0
cCar <-Carmatrix(1,nrow=hh,ncol=1)
cBus <-Busmatrix(1,nrow=hh,ncol=1)
cTaxi <-Taximatrix(1,nrow=hh,ncol=1)
cSharedtaxi <-Sharedtaximatrix(1,nrow=hh,ncol=1)
cAge <-Age*(data$X.Age)/10
cPosAuto <-PosAuto*(data$X.PosAuto)
cDrlicense <-Drlicense*(data$X.Drlicense)
cWork <-Work*(data$X.Work)
cIncome <-Income*(data$X.Income)/100000
cCostCar <-Car*(data$X.CostCar)/100
cCostBus <-Bus*(data$X.CostBus)/100
cCostTaxi <-Taxi*(data$X.CostTaxi)/100
cCostSharedtaxi <-Sharedtaxi*(data$X.CostSharedtaxi)/10
cDistance <-Distance*(data$X.Distance)/10
cTimeCar <-Time*(data$X.CarMin)/10
cTimeBus <-Time*(data$X.BusMin)/10
cTimeTaxi <-Time*(data$X.TaxiMin)/10
cTimeSharedtaxi <-Time*(data$X.SharedtaxiMin)/10
cTimeWalk <-Time*(data$X.WalkMin)/100
vCar <-cCar+cAge+cWork+cIncome+cDistance+cTimeCar
vBus <-cBus+cPosAuto+cAge+cIncome+cCostBus+cDistance+cTimeBus
vTaxi <-cTaxi+cIncome+cCostTaxi+cDistance
vSharedtaxi <-cSharedtaxi+cIncome+cCostSharedtaxi+cDistance
vWalk <-cIncome+cDistance+cTimeWalk
PosCar <-data$X.PosAuto*exp(vCar)
PosBus <-exp(vBus)
PosTaxi <-exp(vTaxi)
PosSharedtaxi <-exp(vSharedtaxi)
PosWalk <-exp(vWalk)
deno <-PosCar+PosBus+PosTaxi+PosSharedtaxi+PosWalk
PCar <-PosCar/deno
PBus <-PosBus/deno
PTaxi <-PosTaxi/deno
PSharedtaxi <-PosSharedtaxi/deno
PWalk <-PosWalk/deno
PCar <-(PCar!=0)*PCar+(PCar==0)
PBus <-(PBus!=0)*PBus+(PBus==0)
PTaxi <-(PTaxi!=0)*PTaxi+(PTaxi==0)
PSharedtaxi <-(PSharedtaxi!=0)*PSharedtaxi+(PSharedtaxi==0)
PWalk <-(PWalk!=0)*PWalk+(PWalk==0)
SCar <-(data$X.FModeFr==1)
SBus <-(data$X.FModeFr==2)
STaxi <-(data$X.FModeFr.==3)
SSharedtaxi <-(data$X.FModeFr.==4)
SWalk <-(data$X.FModeFr.==5)
OutPut <-data.frame(PCar,PBus,PTaxi,PSharedtaxi,PWalk,SCar,SBus,STaxi,SSharedtaxi,SWalk)
write.table(OutPut,"Output.csv",quote=F,col.names=T,append=T,sep=",")
LL <-colSums(SCarlog(PCar)+SBuslog(PBus)+STaxilog(PTaxi)+SSharedtaxilog(PSharedtaxi)+SWalk*log(PWalk))
}
b0 <-numeric(12)
res_BFGS <-optim(b0,fr, method="BFGS", hessian=TRUE, control=list(fnscale=-1))
show <-function(res,b0){
b <-res$par
hhh <-res$hessian
tval <-b/sqrt(-diag(solve(hhh)))
L0 <-hh*log(1/5)
LL <-res$value
print(b)
print(tval)
print(hhh)
print(L0)
print(LL)
print((L0-LL)/L0)
print((L0-(LL-length(b)))/L0)
}
show(res_BFGS,b0)