这个,完整代码,需要数据文件哪~没办法发啊~数据,就是数据导入不进来,您看看吧!多谢
##Read parameters set data
ParamData <- read.csv("Coeffi.csv", header = TRUE)
##count the number of row of ParamData
rr <- nrow(ParamData)
print(rr)
##Read choice set data
Data<-read.csv("Sequence.csv",header=TRUE)
##count the number of row
hh<-nrow(Data)
print(hh)
##number of choice alternative
##ch <- 5
##Initial value of variables
Sd <- sum(Data[,2]== 1); Se <- sum(Data[,2]== 2);Sf <- sum(Data[,2]== 3); Sg <- sum(Data[,2]== 4); Sh <- sum(Data[,2]== 5)
cat("d:",Sd," e:",Se," f:",Sf," g:",Sg," h:",Sh,"\n")
##Transfer data for easily calculation
TData <- t(Data)
ParamData1 <- matrix (ParamData[,2], nrow = 10000, ncol = 1)
ParamData2 <- matrix (ParamData[,3], nrow = 10000, ncol = 1)
TData3 <- matrix(TData[3,], nrow = 1, ncol = 283)
TData4 <- matrix(TData[3,], nrow = 1, ncol = 283)
TData5 <- matrix(TData[5,], nrow = 1, ncol = 283)
TData6 <- matrix(TData[6,], nrow = 1, ncol = 283)
TData7 <- matrix(TData[7,], nrow = 1, ncol = 283)
TData8 <- matrix(TData[8,], nrow = 1, ncol = 283)
TData9 <- matrix(TData[9,], nrow = 1, ncol = 283)
TData10 <- matrix(TData[10,], nrow = 1, ncol = 283)
TData11 <- matrix(TData[11,], nrow = 1, ncol = 283)
TData12 <- matrix(TData[12,], nrow = 1, ncol = 283)
##calculate utility of each alternative
d <- ParamData1 %*% TData3 + ParamData2 %*% TData8
e <- ParamData1 %*% TData4 + ParamData2 %*% TData9
f <- ParamData1 %*% TData5 + ParamData2 %*% TData10
g <- ParamData1 %*% TData6 + ParamData2 %*% TData11
h <- ParamData1 %*% TData7 + ParamData2 %*% TData12
##calculate the expected utility of each alternative
for (i in 1:10000){
for (j in 1:283){
Eh <- 0
## Eg <- h
Ef[j] <- g[i,j] * integrate(dnorm,-Inf, as.numeric(g[i,j] - h[i,j]), mean = 0, sd = 2)$value + h[i,j] * (1-integrate (dnorm,-Inf, as.numeric(g[i,j] - h[i,j]), mean = 0, sd = 2))$value
Ee[j] <- f[i,j] * integrate(dnorm, -Inf, as.numeric(f[i,j] - Ef[i,j]), mean = 0, sd = 3)$value + Ef[i,j] * (1 - integrate(dnorm, -Inf, as.numeric(f[i,j] - Ef[i,j]), mean = 0, sd = 3))$value
Ed[j] <- e[i,j] * integrate(dnorm, -Inf, as.nuemric(e[i,j] - Ee[i,j]), mean = 0, sd = 4)$value + Ee[i,j] * (1 - integrate(dnorm, -Inf, as.numeric(e[i,j] - Ee[i,j]), mean = 0, sd = 4))$value
##Calculate the probability of each alternative at each junction
PPPd[j] = integrate(dnorm, -Inf, d[i,j] - Ed[i,j], mean = 0, sd = 5)$value
PPPe[j] = integrate(dnorm, -Inf, e[i,j] - Ee[i,j], mean = 0, sd = 4)$value
PPPf[j] = integrate(dnorm, -Inf, f[i,j] - Ef[i,j], mean = 0, sd = 3)$value
PPPg[j] = integrate(dnorm, -Inf, g[i,j] - h[i,j], mean = 0, sd = 2)$value
PPPh[j] <- 1
##Calculate the final probability of each alternative
PPd <- integrate(dnorm, -Inf, d[i,j] - Ed[i,j], mean = 0, sd = 5)$value
PPe <- PPPe[j] * (1 - PPPd[j])
PPf <- PPPf[j] * (1 - PPPd[j]) * (1 - PPPe[j])
PPg <- PPPg[j] * (1 - PPPd[j]) * (1 - PPPe[j]) * (1 - PPPf[j])
PPh <- PPPh[j] * (1 - PPPd[j]) * (1 - PPPe[j]) * (1 - PPPf[j]) * (1 - PPPg[j])
##
Pd[j] <- (PPd !=0)*PPd + (PPd==0)
Pe[j] <- (PPe !=0)*PPe + (PPe==0)
Pf[j] <- (PPf !=0)*PPf + (PPf==0)
Pg[j] <- (PPg !=0)*PPg + (PPg==0)
Ph[j] <- (PPh !=0)*PPh + (PPh==0)
##selection result
Cd <- Data[,2]== 1
Ce <- Data[,2]== 2
Cf <- Data[,2]== 3
Cg <- Data[,2]== 4
Ch <- Data[,2]== 5
##log-likelihood function
LL[i] <- colSums( Cd*log(Pd) + Ce*log(Pe) + Cf*log(Pf) + Cg*log(Pg) + Ch*log(Ph) )
j = j+1
}
i = i+1
}