<< Statistics`ContinuousDistributions`

Off[General::spell1]

Intervals[m_, s_, P_]:= 
  {CDF[NormalDistribution[m, s], P[[1]]], 
   CDF[NormalDistribution[m, s], P[[3]]] - 
     CDF[NormalDistribution[m, s], P[[1]]], 
   CDF[NormalDistribution[m, s], P[[2]]] - 
     CDF[NormalDistribution[m, s], P[[3]]], 
   1 - CDF[NormalDistribution[m, s], P[[2]]]}

alpha0[m_, s_]:= Random[NormalDistribution[m, s]]

iChoices[m_, s_, si_, P_]:= 
   Module[{a0, ak, a},
          a0 = alpha0[m, s];
          ak = Table[Random[NormalDistribution[0, si]], {i, 1, 3}];
          a = a0 + ak;
          Which[a[[1]] < P[[1]] && a[[2]] < P[[2]] && a[[3]] < P[[3]], 1,
                a[[1]] < P[[1]] && a[[2]] > P[[2]] && a[[3]] < P[[3]], 2,
                a[[1]] > P[[1]] && a[[2]] < P[[2]] && a[[3]] < P[[3]], 3,
                a[[1]] > P[[1]] && a[[2]] < P[[2]] && a[[3]] > P[[3]], 4,
                a[[1]] < P[[1]] && a[[2]] > P[[2]] && a[[3]] > P[[3]], 5,
                a[[1]] > P[[1]] && a[[2]] > P[[2]] && a[[3]] > P[[3]], 6,
                a[[1]] < P[[1]] && a[[2]] < P[[2]] && a[[3]] > P[[3]], 7,
                a[[1]] > P[[1]] && a[[2]] > P[[2]] && a[[3]] < P[[3]], 8]]

ChoiceDist[m_, s_, si_, P_]:= 
   Module[{c},
          c = Table[iChoices[m, s, si, P], {i, 1, 100}];
          Map[Length, Table[Flatten[Position[c, k]], {k, 1, 8}]]]

a0i[m_, s_]:= 
  Table[Quantile[NormalDistribution[m, s], 0.005 + 0.01 * i], {i, 0, 99}]

iProbs[m_, s_, si_, P_]:= 
   Module[{a0},
          a0 = a0i[m, s];
          Table[{N[CDF[NormalDistribution[a0[[i]], si], P[[1]]] * 
                   CDF[NormalDistribution[a0[[i]], si], P[[2]]] * 
                   CDF[NormalDistribution[a0[[i]], si], P[[3]]]], 
                 N[CDF[NormalDistribution[a0[[i]], si], P[[1]]] * 
                   (1 - CDF[NormalDistribution[a0[[i]], si], P[[2]]]) * 
                   CDF[NormalDistribution[a0[[i]], si], P[[3]]]], 
                 N[(1 - CDF[NormalDistribution[a0[[i]], si], P[[1]]]) * 
                   CDF[NormalDistribution[a0[[i]], si], P[[2]]] * 
                   CDF[NormalDistribution[a0[[i]], si], P[[3]]]], 
                 N[(1 - CDF[NormalDistribution[a0[[i]], si], P[[1]]]) * 
                   CDF[NormalDistribution[a0[[i]], si], P[[2]]] * 
                   (1 - CDF[NormalDistribution[a0[[i]], si], P[[3]]])], 
                 N[CDF[NormalDistribution[a0[[i]], si], P[[1]]] * 
                   (1 - CDF[NormalDistribution[a0[[i]], si], P[[2]]]) * 
                   (1 - CDF[NormalDistribution[a0[[i]], si], P[[3]]])], 
                 N[(1 - CDF[NormalDistribution[a0[[i]], si], P[[1]]]) * 
                   (1 - CDF[NormalDistribution[a0[[i]], si], P[[2]]]) * 
                   (1 - CDF[NormalDistribution[a0[[i]], si], P[[3]]])], 
                 N[CDF[NormalDistribution[a0[[i]], si], P[[1]]] * 
                   CDF[NormalDistribution[a0[[i]], si], P[[2]]] * 
                   (1 - CDF[NormalDistribution[a0[[i]], si], P[[3]]])], 
                 N[(1 - CDF[NormalDistribution[a0[[i]], si], P[[1]]]) * 
                   (1 - CDF[NormalDistribution[a0[[i]], si], P[[2]]]) * 
                   CDF[NormalDistribution[a0[[i]], si], P[[3]]]]}, 
                 {i, 1, 100}]]

Probs[m_, s_, si_, P_] := Map[Mean, Transpose[iProbs[m, s, si, P]]]

LnL1[m_, s_, si_]:= 
  Module[{P1, LOGp, d, LogL}, 
         P1 = {0.145, 0.263, 0.200};
         LOGp = Log[Probs[m, s, si, P1]];
         d1 = {6, 5, 16, 9, 4, 42, 2, 16};
         LogL = d1.LOGp]

LnL4[m_, s_, si_]:= 
  Module[{P4, LOGp, d, LogL}, 
         P4 = {0.205, -0.564, 0.075};
         LOGp = Log[Probs[m, s, si, P1]];
         d4 = {1, 8,  5, 9, 4, 50, 0, 23};
         LogL = d4.LOGp]

LnL6[m_, s_, si_]:= 
  Module[{P6, LOGp, d, LogL}, 
         P6 = {0.295, 0.370, 0.339};
         LOGp = Log[Probs[m, s, si, P6]];
         d6 = {12, 3, 16, 3, 5, 44, 2, 15};
         LogL = d6.LOGp]

LnL[m_, s_, si_]:= LnL1[m, s, si] + LnL4[m, s, si]

Table[LnL1[0.341, 0.160, 0.222], {i, -2, 2}]