(Local) In[]:=
β={0.25,0.5,-0.75};
(Local) In[]:=
samples_x=Flatten[Table[{1,x1,x2},{x1,-2,2,0.05},{x2,-2,2,0.05}],1];​​Length[samples_x]
(Local) Out[]=
6561
(Local) In[]:=
randomVariables=RandomVariate[LogisticDistribution[0,1],Length[samples_x]];​​Show[Histogram[randomVariables,50,"ProbabilityDensity"],Plot[PDF[LogisticDistribution[0,1],x],{x,-6,6}]]
(Local) Out[]=
-5
0
5
0.00
0.05
0.10
0.15
0.20
0.25
(Local) In[]:=
samples_y=Table[If[randomVariables[[i]]<samples_x[[i]].β,1,0],{i,Length[samples_x]}];
(Local) In[]:=
samples_xy=Table[Flatten[{samples_x[[i]],samples_y[[i]]}],{i,Length[samples_x]}];​​Length[samples_xy]
(Local) Out[]=
6561
(Local) In[]:=
SmoothHistogram3D[Map[#[[2;;3]]&,Select[samples_xy,#[[4]]==1&]]]
(Local) Out[]=
(Local) In[]:=
Φ[x_]=CDF[LogisticDistribution[0,1],x]
(Local) Out[]=
1
1+
-x

(Local) In[]:=
logLikelihood[β_]:=Sum[samples_y[[i]]*Log[Φ[samples_x[[i]].β]]+(1-samples_y[[i]])*Log[1-Φ[samples_x[[i]].β]],{i,Length[samples_x]}]
(Local) In[]:=
NMaximize[logLikelihood[{β0,β1,β2}],{β0,β1,β2}]
(Local) Out[]=
{-3811.08,{β00.219385,β10.507949,β2-0.76393}}