﻿### Multinomial Logit model estimation

### データファイルの読み込み
Data <- read.csv("all.csv",header=T)

## データ数:Dataの行数を数える
hh <- nrow(Data)
print(hh)

## パラメータの初期値の設定
b0 <- numeric(11) 

##### Logit model の対数尤度関数の定義 #####

fr <- function(x) {
  ### パラメータの宣言：
  ## 定数項
  b1 <- x[1]
  b2 <- x[2]
  b3 <- x[3]
  b4 <- x[4]
  
  ## 目的地までの所要時間
  d1 <- x[5]
  

 
  ## 料金

  f1 <- x[6]
  

  ## 朝ダミー

  m1 <- x[7]
  

  ## 夜ダミー

  n1 <- x[8]
  

  ## 出勤ダミー

  j1 <- x[9]
  

  ## 散歩回遊ダミー

  w1 <- x[10]
  
  
  ## 買い物ダミー

  s1 <- x[11]
 


  
  ## 対数尤度のための変数を宣言
  LL = 0

  ### 今回用いる交通手段は以下の3つ．
  ## 鉄道(train)
  ## バス(bus)
  ## 自動車(car)
  ## 自転車(bike)
  ## 徒歩(walk)
    
  ## 効用の計算:説明変数にしたい列を入れる．
                                              # 時間                       # 料金              # 定数項
  train  <- Data$代替手段生成可否train*exp(d1*Data$総所要時間train +f1*Data$費用train +m1*Data$朝ダミー +s1*Data$買い物ダミー +n1*Data$夜ダミー +j1*Data$出勤ダミー + b1*matrix(1,nrow =hh,ncol=1))
  bus    <- Data$代替手段生成可否bus  *exp(d1*Data$総所要時間bus   +f1*Data$費用bus   + b2*matrix(1,nrow =hh,ncol=1))
  car    <- Data$代替手段生成可否car  *exp(d1*Data$所要時間car +f1*Data$費用car                      + b3*matrix(1,nrow =hh,ncol=1))
  bike   <- Data$代替手段生成可否bike *exp(d1*Data$所要時間bike                     + b4*matrix(1,nrow =hh,ncol=1))
  walk   <- Data$代替手段生成可否walk *exp(d1*Data$所要時間walk +w1*Data$散歩回遊ダミー 　)
 

  ### 選択確率の計算
  ## 分母となる，各々のexp(V)の和をつくる
  deno <- (car + train + bus + bike + walk)

  ## それぞれ計算する
  Ptrain <- Data$代替手段生成可否train*(train / deno)
  Pbus   <- Data$代替手段生成可否bus  *(bus   / deno)
  Pcar   <- Data$代替手段生成可否car  *(car   / deno)
  Pbike  <- Data$代替手段生成可否bike *(bike  / deno)
  Pwalk  <- Data$代替手段生成可否walk *(walk  / deno)

  ## 選択確率が0になってしまった場合に起こる問題の回避
   
　 Ptrain <- (Ptrain!=0)*Ptrain + (Ptrain==0)
    
　 Pbus   <- (Pbus!=0)*Pbus   + (Pbus==0)
    
   Pcar   <- (Pcar!=0)*Pcar   + (Pcar==0)
   
   Pbike  <- (Pbike!=0)*Pbike + (Pbike==0)
   Pwalk  <- (Pwalk!=0)*Pwalk + (Pwalk==0)
  


  ## 選択結果
  Ctrain   <- Data$代表交通手段コード =="200"
  Cbus     <- Data$代表交通手段コード =="240"
  Ccar     <- Data$代表交通手段コード =="100"
  Cbike    <- Data$代表交通手段コード =="410"
  Cwalk    <- Data$代表交通手段コード =="420"
  

  ## 対数尤度の計算
  LL <- colSums(Ctrain*log(Ptrain) + Cbus*log(Pbus) + Ccar*log(Pcar) + Cbike*log(Pbike) + Cwalk*log(Pwalk))

}

##### 対数尤度関数frの最大化#####

##パラメータ値の最適化 
res <- optim(b0,fr, method = "Nelder-Mead", hessian = TRUE, control=list(fnscale=-1))

## パラメータ推定値、ヘッセ行列 
b   <- res$par
hhh <- res$hessian

## t値の計算
tval <- b/sqrt(-diag(solve(hhh)))

## 初期尤度
L0 <- fr(b0)

## 最終尤度
LL <- res$value

##### 結果の出力 #####
print(res)

## 初期尤度
print(L0)

## 最終尤度
print(LL)

##ρ^2値
print((L0-LL)/L0)

## 修正済ρ^2値
print((L0-(LL-length(b)))/L0)

##パラメータ推定値 
print(b)

## t値 
print(tval)

