(Local) In[]:=
β={-0.25,0.5,-0.5};
(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[NormalDistribution[0,1],Length[samples_x]];
(Local) In[]:=
Show[Histogram[randomVariables,30,"ProbabilityDensity"],Plot[PDF[NormalDistribution[0,1],x],{x,-9,9},PlotRange{0,Full}]]
(Local) Out[]=
-2
0
2
4
0.0
0.1
0.2
0.3
0.4
(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[NormalDistribution[0,1],x]
(Local) Out[]=
1
2
Erfc-
x
2

(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[]=
{-3492.18,{β0-0.236155,β10.486374,β2-0.494407}}