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.

Comments:

  • 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);
      }
    
      ////////// - NelderMeadSimplex - /////////
    
      // 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);
      }
    

    }