##多项logit模型计算

##读入数据文件

data<-read.csv("E:\\毕设\\R\\suda4.csv",header=T)

hh<-nrow(data)

print(hh)

ch<-7

b0<-c(0,0,0,0,0,0)

Sfoot<-sum(data[,2]==1);Sbike<-sum(data[,2]==2);Sebike<-sum(data[,2]==3);Sbus<-sum(data[,2]==4);Scar<-sum(data[,2]==5);Staxi<-sum(data[,2]==6);Sfbus<-sum(data[,2]==7)

cat("foot:",Sfoot,"bike:",Sbike,"ebike:",Sebike,"bus:",Sbus,"taxi:",Staxi,"car:",Scar,"fbus:",Sfbus)

##定义Logit模型的对数似然函数

fr<-function(x){

LL=0

##效用计算

foot<-x[40]*data[,9]+x[1]*data[,22]+x[2]*data[,23]+x[3]*data[,24]+x[4]*data[,25]+x[5]*data[,26]+x[6]*data[,27]+x[7]*data[,28]+x[8]*data[,29]+x[9]*data[,30]+x[10]*data[,31]+x[11]*data[,32]+x[12]*data[,33]+x[13]*data[,34]+x[15]*data[,35]+x[16]*data[,36]+x[53]*matrix(1,nrow=hh,ncol=1)

bike<-x[41]*data[,10]+x[1]*data[,22]+x[2]*data[,23]+x[3]*data[,24]+x[4]*data[,25]+x[5]*data[,26]+x[6]*data[,27]+x[7]*data[,28]+x[8]*data[,29]+x[9]*data[,30]+x[10]*data[,31]+x[11]*data[,32]+x[12]*data[,33]+x[13]*data[,34]+x[15]*data[,35]+x[16]*data[,36]+x[54]*matrix(1,nrow=hh,ncol=1)

ebike<-x[42]*data[,11]+x[1]*data[,22]+x[2]*data[,23]+x[3]*data[,24]+x[4]*data[,25]+x[5]*data[,26]+x[6]*data[,27]+x[7]*data[,28]+x[8]*data[,29]+x[9]*data[,30]+x[10]*data[,31]+x[11]*data[,32]+x[12]*data[,33]+x[13]*data[,34]+x[15]*data[,35]+x[16]*data[,36]+x[55]*matrix(1,nrow=hh,ncol=1)

bus<-x[43]*data[,12]+x[44]*data[,13]+x[45]*data[,14]+x[1]*data[,22]+x[2]*data[,23]+x[3]*data[,24]+x[4]*data[,25]+x[5]*data[,26]+x[6]*data[,27]+x[7]*data[,28]+x[8]*data[,29]+x[9]*data[,30]+x[10]*data[,31]+x[11]*data[,32]+x[12]*data[,33]+x[13]*data[,34]+x[15]*data[,35]+x[16]*data[,36] +x[56]*matrix(1,nrow=hh,ncol=1)

taxi<-x[46]*data[,15]+x[47]*data[,16]+x[48]*data[,17]+x[1]*data[,22]+x[2]*data[,23]+x[3]*data[,24]+x[4]*data[,25]+x[5]*data[,26]+x[6]*data[,27]+x[7]*data[,28]+x[8]*data[,29]+x[9]*data[,30]+x[10]*data[,31]+x[11]*data[,32]+x[12]*data[,33]+x[13]*data[,34]+x[15]*data[,35]+x[16]*data[,36]+x[57]*matrix(1,nrow=hh,ncol=1)

car<-x[49]*data[,18]+x[50]*data[,19]+x[1]*data[,22]+x[2]*data[,23]+x[3]*data[,24]+x[4]*data[,25]+x[5]*data[,26]+x[6]*data[,27]+x[7]*data[,28]+x[8]*data[,29]+x[9]*data[,30]+x[10]*data[,31]+x[11]*data[,32]+x[12]*data[,33]+x[13]*data[,34]+x[15]*data[,35]+x[16]*data[,36]+x[58]*matrix(1,nrow=hh,ncol=1)

fbus<-x[51]*data[,20]+x[52]*data[,21]+x[1]*data[,22]+x[2]*data[,23]+x[3]*data[,24]+x[4]*data[,25]+x[5]*data[,26]+x[6]*data[,27]+x[7]*data[,28]+x[8]*data[,29]+x[9]*data[,30]+x[10]*data[,31]+x[11]*data[,32]+x[12]*data[,33]+x[13]*data[,34]+x[15]*data[,35]+x[16]*data[,36]

##效用函数的指数化

Efoot<-exp(foot)*data[,3]

Ebike<-exp(bike)*data[,4]

Eebike<-exp(ebike)*data[,5]

Ebus<-exp(bus)*data[,6]

Etaxi<-exp(taxi)*data[,7]

Ecar<-exp(car)*data[,8]

Efbus<-exp(fbus)*data[,9]

PPfoot<-Efoot/(Efoot+Ebike+Eebike+Ebus+Etaxi+Ecar+Efbus)

PPbike<-Ebike/(Efoot+Ebike+Eebike+Ebus+Etaxi+Ecar+Efbus)

PPebike<-Eebike/(Efoot+Ebike+Eebike+Ebus+Etaxi+Ecar+Efbus)

PPbus<-Ebus/(Efoot+Ebike+Eebike+Ebus+Etaxi+Ecar+Efbus)

PPtaxi<-Etaxi/(Efoot+Ebike+Eebike+Ebus+Etaxi+Ecar+Efbus)

PPcar<-Ecar/(Efoot+Ebike+Eebike+Ebus+Etaxi+Ecar+Efbus)

PPfbus<-Ebus/(Efoot+Ebike+Eebike+Ebus+Etaxi+Ecar+Efbus)

##仅将选择结果进行确实率的有效化

Pfoot<-(PPfoot!=0)*PPfoot+(PPfoot==0)

Pbike<-(PPbike!=0)*PPbus+(PPbus==0)

Pebike<-(PPebike!=0)*PPebike+(PPebike==0)

Pbus<-(PPbus!=0)*PPbus+(PPbus==0)

Ptaxi<-(PPtaxi!=0)*PPtaxi+(PPtaxi==0)

Pcar<-(PPcar!=0)*PPcar+(PPcar==0)

Pfbus<-(PPfbus!=0)*PPfbus+(PPfbus==0)

##选择结果

Cfoot<-data[,2]==1

Cbike<-data[,2]==2

Cebike<-data[,2]==3

Cbus<-data[,2]==4

Ctaxi<-data[,2]==5

Ccar<-data[,2]==6

Cfbus<-data[,2]==7

##对数似然值的计算

LL<-colSums(Cfoot*log(Pfoot)+Cbike*log(Pbike)+Cebike*log(Pebike)+Cbus*log(Pbus)+Ctaxi*log(Ptaxi)+Ccar*log(Pcar)+Cfbus*log(Pfbus))

return(LL)

}

##最大似然值的计算

res<-optim(b0,fr,method="BFGS",hessian=TRUE,control=list(fnscale=-1))

##参数估计

b<-res$par

hhh<-res$hessian

##t值的计算

tval<-b/sqrt(-diag(solve(hhh)))

##初步似然值

L0<-Sfoot*log(Sfoot/hh)+Sbike*log(Sbike/hh)+Sebike*log(Sebike/hh)+Sbus*log(Sbus/hh)+Staxi*log(Staxi/hh)+Scar*log(Scar/hh)+Sfbus*log(Sfbus/hh)

##最终似然值

LL<-res$value

##拟合度计算

##结果输出

##p^2的值

cat("roh = ",(L0-LL)/L0)

##修正的p^2值

cat("rohbar=",(L0-(LL-length(b)))/L0)

print(res)

print(tval)

错误于optim(b0, fr, method = "BFGS", hessian = TRUE, control = list(fnscale = -1)) :

'vmmin'的初始值不能为无穷大

>

这个代码是仿照教授给我的代码编的，用原来的就不会有这个问题，求大神解决！

小弟先谢过各位