# Probit Regression

#1

Is possible to make a Probit Regression with Math.NET?

(Christoph Rüegg) #2

Not out of the box.

For reference: https://en.wikipedia.org/wiki/Probit_model

(Peter Vanderwaart) #3

This aroused my curiosity. I’ve never had any exposure to probit models in any context, but I thought I’d give it a try. I based the program on this paper: On Binary Choice Models: Logit and Probit (PDF). by Thomas B. Fomby at SMU. It has the advantage of simplicity. It’s contents are all I know about Probit Models.

• The sample data was constructed to be strongly polar so I could easily see if the program was giving plausible answers. It uses a random number generator for noise, so no two runs will be exactly the same.

• The dot product function in function F produces the same calculation as in the paper, just different object types and notation.

• The use of the ? operator in function ll provides the same “either/or” as the multiplication by Y[i] in the paper. I originally wrote the function as in the paper, but the CDF in function F returned 1.0 for extreme values leading to a Math.Log(0) error in the function ll.

• I have not tried any other data set, nor have I compared the results to output from any other program.

private void btnProbit_Click(object sender, EventArgs e)
{
// data for Probit Model
// two independent variable of hard data plus one indepenent variable of noise
Random RanGen = new Random();

``````  int nObs = 21; // number of observations
int nIv = 3; // number of independent variables

Vector<double>[] X = new Vector<double>[nObs];
Vector<double> Y = Vector<double>.Build.Dense(nObs);

// create test data

for (int i = 0; i < nObs; i++)
{
Vector<double> v = Vector<double>.Build.Dense(nIv);
v[0] = Convert.ToDouble(i)+ 5.0*RanGen.NextDouble();
v[1] = Convert.ToDouble(20 - i)+ 5.0*RanGen.NextDouble();
v[2] = RanGen.NextDouble();
X[i] = v;
Y[i] = (i <= 10 ? 1.0 : 0.0);
}

// define functions
var F = new Func<Vector<double>, int, double>((b, i) => Normal.CDF(0.0, 1.0, b.DotProduct(X[i])));
var ll = new Func<Vector<double>, int, double>((b, i) => (Y[i] > 0.5 ? Math.Log(F(b, i)) : Math.Log(1.0 - F(b, i))));
var LL = new Func<Vector<double>, double>((b) => -1.0 * Y.MapIndexed((i, y) => ll(b, i)).Sum()); //Y used only as a source of indexes

// solve
var obj = ObjectiveFunction.Value(LL);
var solver = new NelderMeadSimplex(1e-3, maximumIterations: 1000);
var initialGuess = new DenseVector(new[] { 0.0, 0.0, 0.0 });

var result = solver.FindMinimum(obj, initialGuess);

try
{

Console.WriteLine("Iterations: " + result.Iterations.ToString());
Console.WriteLine("Minimum: " + result.MinimizingPoint.ToString());
Console.WriteLine("Reason: " + result.ReasonForExit.ToString());
Console.WriteLine("");
Console.WriteLine("x[1],  x[2],  x[3],   y,   F(x)");

for (int i = 0; i < nObs; i++)
{
for (int j = 0; j < nIv; j++)
{
Console.Write(X[i][j].ToString("#0.0##") + ", ");
}
double rr = F(result.MinimizingPoint, i);
Console.WriteLine(Y[i].ToString("#0.0##") + ", " + rr.ToString("#0.0####"));
}

}
catch (Exception eSolve)
{
Console.WriteLine(eSolve.Message);
}
``````

}