Параметры программы наименьших квадратов не меняются - PullRequest
0 голосов
/ 22 мая 2019

Я пытаюсь создать программу, которая может аппроксимировать параметры для экспоненциальной функции вида y = ae^bx, которая соответствует любому заданному набору данных, я посмотрел вверх и попытался понять и использовать алгоритм Марквардта-Левенберга, хотя в В моей реализации параметры на каждой итерации не меняются. Может кто-нибудь помочь указать, что не так с моей реализацией, или я просто не до конца понимаю алгоритм?

public class regression {

  boolean exponential;
  matrix data, estimate, gradient, hessian, diagHessian;
  double scalar;

  public regression(boolean exponential, double[][] data) {
    this.exponential = exponential;
    this.data = new matrix(data.length, data[0].length);
    this.data.copy(data);
  }

  public void approximate() {

    this.findDerivative(this.exponential);
    for (int i = 0; i < 10000; i++) {
      double s = this.evaluateError();
      double x = this.currentError();

      //if new error is smaller than old error 
      if (s < x) {
        this.estimate = this.evaluate();
        reload();
        this.scalar = this.scalar / 10;
      } else {
        this.scalar = this.scalar * 10.0;
      }
    }
  }

  public matrix evaluate() {
    matrix a = new matrix(2, 2);
    matrix b = new matrix(2, 2);

    b = this.diagHessian;
    b.times(this.scalar);

    a = this.hessian.subtract(b);

    matrix r = a.inverse();
    matrix x = this.gradient;
    matrix d = r.multiply(x);
    matrix y = this.estimate.subtract(d);

    return y;
  }

  public double currentError() {
    double error = 0.0;
    double w = this.estimate.arrayM[0][0];
    double z = this.estimate.arrayM[1][0];

    for (int i = 0; i < this.data.size; i++) {
      error += (this.data.arrayM[i][1] -
        (w * java.lang.Math.exp(z * this.data.arrayM[i][0]))) * (this.data.arrayM[i]
        [1] - (w * java.lang.Math.exp(z * this.data.arrayM[i][0])));
    }

    return error * 0.5;
  }

  //evaluate the estimate matrix at the next iteration and 
  //return a new matrix not effecting the current iteration
  public double evaluateError() {
    matrix a = new matrix(2, 2);
    matrix b = new matrix(2, 2);

    b = this.diagHessian;
    b.times(this.scalar);

    a = this.hessian.subtract(b);

    matrix r = a.inverse();
    matrix x = this.gradient;
    matrix d = r.multiply(x);
    matrix y = this.estimate.subtract(d);

    double error = 0.0;
    double w = y.arrayM[0][0];
    double z = y.arrayM[1][0];

    for (int i = 0; i < this.data.size; i++) {
      error += (this.data.arrayM[i][1] -
        (w * java.lang.Math.exp(z * this.data.arrayM[i][0]))) * (this.data.arrayM[i]
        [1] - (w * java.lang.Math.exp(z * this.data.arrayM[i][0])));
    }

    return error * 0.5;
  }

  public static matrix transpose(matrix b) {
    matrix a = new matrix(b.sizeX, b.size);

    for (int i = 0; i < b.size; i++) {
      for (int y = 0; y < b.sizeX; y++) {
        a.arrayM[y][i] = b.arrayM[i][y];
      }
    }
    return a;
  }

  public void reload() {
    this.calcGradient();
    this.calcHessian();
    this.calcDiag();
  }


  public void findDerivative(boolean expo) {
    if (expo) {
      //create estimate
      this.estimate = new matrix(2, 1);
      double a = 40000;
      double b = 0.00012;
      this.estimate.arrayM[0][0] = a;
      this.estimate.arrayM[1][0] = b;
      this.gradient = new matrix(2, 1);
      this.calcGradient();
      this.hessian = new matrix(2, 2);
      this.calcHessian();
      this.diagHessian = new matrix(2, 2);
      this.calcDiag();
      this.scalar = 0.5;
    }
  }

  public void calcDiag() {
    for (int i = 0; i < this.hessian.size; i++) {
      this.diagHessian.arrayM[i][i] = 1.0;
    }
  }

  public void calcGradient() {
    if (this.exponential) {
      double a = this.estimate.arrayM[0][0];
      double b = this.estimate.arrayM[1][0];
      double s = 0.0;

      for (int i = 0; i < this.data.size; i++) {
        s += -(this.data.arrayM[i][1] -
            a * java.lang.Math.exp(b * data.arrayM[i][0])) *
          (java.lang.Math.exp(b * data.arrayM[i][0]));
      }

      this.gradient.arrayM[0][0] = s;
      s = 0.0;

      for (int i = 0; i < this.data.size; i++) {
        s +=
          -(this.data.arrayM[i][1] - a * java.lang.Math.exp(b * data.arrayM[i][0])) *
          (this.data.arrayM[i][0] * a * java.lang.Math.exp(b * data.arrayM[i][0]));
      }
      this.gradient.arrayM[1][0] = s;
    }
  }

  public void calcHessian() {
    if (this.exponential) {
      double a = this.estimate.arrayM[0][0];
      double b = this.estimate.arrayM[1][0];
      double s = 0.0;
      matrix jacobian = new matrix(this.data.size, 2);

      for (int i = 0; i < jacobian.size; i++) {

        jacobian.arrayM[i][0] = Math.exp(b * this.data.arrayM[i][0]);
        jacobian.arrayM[i][1] = a * this.data.arrayM[i[0] * Math.exp(b * this.data.arrayM[i][0]);
      }
      matrix c = transpose(jacobian);
      matrix d = c.multiply(jacobian);
      this.hessian = d;
    }
  }
}
...