Python scipy.optimize.leastsq 到 Java org.apache.commons.math3.fitting.leastsquares
Python scipy.optimize.leastsq to Java org.apache.commons.math3.fitting.leastsquares
我试着模仿这个algorithm, developed in Python, that calculates geolocation based on seen Wifi stations positions, itself based on this idea。
该算法首先使用 Numpy 函数来计算观测到的纬度和经度的基本加权平均值。为了尽量减少可能的 Wifi 位置错误的影响，它还使用“scipy.optimize.leastsq”方法以统计方式计算，如果可能的话，更精确的位置。
我想在 Java Android 平台上实现相同的行为。
对于所有其他计算，我成功地依赖 org.apache.commons.math3。
所以对于最小二乘问题，我逻辑上尝试依赖 https://commons.apache.org/proper/commons-math/userguide/leastsquares.html.
我的问题，如果我很好理解的话，是 Scipy 为我管理了雅可比函数定义的复杂性，而我糟糕的数学技能不允许我正确定义最小二乘问题的模型。我尝试了一些基于此 example 的实验，这似乎接近我需要的东西，但结果并不好，因为我不知道如何处理 "jacobian" 部分。
有人为thispost做了同样的事情，有人可以为我做同样的事情并尝试以简单的方式解释吗？
有关 Python 部分工作原理的更多详细信息：
使用的“scipy.optimize.leastsq”语句是：
(lat, lon), cov_x, info, mesg, ier =
scipy.optimize.leastsq(func, initial, args=data, full_output=True)
其中数据为：latitude/longitude/age中的milliseconds/signal强度，例如：data = numpy.array([(43.48932915, 1.66561772, 1000, -20), (43.48849093, 1.6648176, 2000, -10), (43.48818612, 1.66615113, 3000, -50)])
初始计算的是加权平均值latitude/longitude，本例中：initial = 43.48864654, 1.66550075
函数是
def func(initial, data):
return numpy.array([
geographic distance((float(point[latitude]), float(point[longitude])), (initial[latitude], initial[longitude])).meters * min(math.sqrt(2000.0 / float(point[age])), 1.0) / math.pow(float(point[signal strength]), 2)
结果是：43.4885401095, 1.6648660983
我在 Java 中的实验，我替换了数据值并更改了 "modelI" 的计算方式。我简化了信号强度和年龄值。但事实证明，这还不够。
double modelI = calculateVincentyDistance(o.getY(), o.getX(), center.getY(), center.getX())* Math.min(Math.sqrt(2000.0/1000.0), 1.0) / Math.pow(-10, 2);
我也打算试试https://github.com/odinsbane/least-squares-in-java，但我不确定是否正确使用它，因为我不掌握它的工作方式。
仅供参考，我使用 Vincenty 距离计算，例如可以用 Haversine 或 Euclidean 代替。
感谢您的帮助！
代码不易移植，因为 SciPy 提供了更通用的 Least-squares minimization interface while Apache Commons Math provides curve fitting。还有许多优化问题可以重述为曲线拟合。在 Python 代码中最小化
F(current_point) = Sum{ (distance(known_point[i], current_point) * weight[i])^2 } -> min
Java曲线拟合问题有点不同：
F(current_point) = Sum{ (target_value[i] - model[i](current_point))^2 } -> min
因此，可以通过将所有 target_value 分配给 0 并使 model[i] 计算从 current_point 到 known_point[i] 的加权距离来创建等效拟合问题。
在一般情况下，此类问题无法使用公式进行精确求解，而会使用一些数值优化方法。这里还有另一个区别：Java 实现明确要求您为优化器提供方法来计算被优化函数的导数。如果未提供 Dfun，Python 代码似乎使用了某种差异区分符。您可以在 Java 中手动或使用 FiniteDifferencesDifferentiator but for simple formulas it might be easier to code them explicitly using DerivativeStructure
执行类似的操作
static class PositionInfo {
public final double latitude;
public final double longitude;
public final int ageMs;
public final int strength;
public PositionInfo(double latitude, double longitude, int ageMs, int strength) {
this.latitude = latitude;
this.longitude = longitude;
this.ageMs = ageMs;
this.strength = strength;
}
public double getWeight() {
return Math.min(1.0, Math.sqrt(2000.0 / ageMs)) / (strength * strength);
}
}
static DerivativeStructure getWeightedEuclideanDistance(double tgtLat, double tgtLong, PositionInfo knownPos) {
DerivativeStructure varLat = new DerivativeStructure(2, 1, 0, tgtLat); // latitude is 0-th variable of 2 for derivatives up to 1
DerivativeStructure varLong = new DerivativeStructure(2, 1, 1, tgtLong); // longitude is 1-st variable of 2 for derivatives up to 1
DerivativeStructure latDif = varLat.subtract(knownPos.latitude);
DerivativeStructure longDif = varLong.subtract(knownPos.longitude);
DerivativeStructure latDif2 = latDif.pow(2);
DerivativeStructure longDif2 = longDif.pow(2);
DerivativeStructure dist2 = latDif2.add(longDif2);
DerivativeStructure dist = dist2.sqrt();
return dist.multiply(knownPos.getWeight());
}
// as in https://en.wikipedia.org/wiki/Haversine_formula
static DerivativeStructure getWeightedHaversineDistance(double tgtLat, double tgtLong, PositionInfo knownPos) {
DerivativeStructure varLat = new DerivativeStructure(2, 1, 0, tgtLat);
DerivativeStructure varLong = new DerivativeStructure(2, 1, 1, tgtLong);
DerivativeStructure varLatRad = varLat.toRadians();
DerivativeStructure varLongRad = varLong.toRadians();
DerivativeStructure latDifRad2 = varLat.subtract(knownPos.latitude).toRadians().divide(2);
DerivativeStructure longDifRad2 = varLong.subtract(knownPos.longitude).toRadians().divide(2);
DerivativeStructure sinLat2 = latDifRad2.sin().pow(2);
DerivativeStructure sinLong2 = longDifRad2.sin().pow(2);
DerivativeStructure summand2 = varLatRad.cos().multiply(varLongRad.cos()).multiply(sinLong2);
DerivativeStructure sum = sinLat2.add(summand2);
DerivativeStructure dist = sum.sqrt().asin();
return dist.multiply(knownPos.getWeight());
}
使用这样的准备你可以做这样的事情：
public static void main(String[] args) {
// latitude/longitude/age in milliseconds/signal strength
final PositionInfo[] data = new PositionInfo[]{
new PositionInfo(43.48932915, 1.66561772, 1000, -20),
new PositionInfo(43.48849093, 1.6648176, 2000, -10),
new PositionInfo(43.48818612, 1.66615113, 3000, -50)
};
double[] target = new double[data.length];
Arrays.fill(target, 0.0);
double[] start = new double[2];
for (PositionInfo row : data) {
start[0] += row.latitude;
start[1] += row.longitude;
}
start[0] /= data.length;
start[1] /= data.length;
MultivariateJacobianFunction distancesModel = new MultivariateJacobianFunction() {
@Override
public Pair<RealVector, RealMatrix> value(final RealVector point) {
double tgtLat = point.getEntry(0);
double tgtLong = point.getEntry(1);
RealVector value = new ArrayRealVector(data.length);
RealMatrix jacobian = new Array2DRowRealMatrix(data.length, 2);
for (int i = 0; i < data.length; i++) {
DerivativeStructure distance = getWeightedEuclideanDistance(tgtLat, tgtLong, data[i]);
//DerivativeStructure distance = getWeightedHaversineDistance(tgtLat, tgtLong, data[i]);
value.setEntry(i, distance.getValue());
jacobian.setEntry(i, 0, distance.getPartialDerivative(1, 0));
jacobian.setEntry(i, 1, distance.getPartialDerivative(0, 1));
}
return new Pair<RealVector, RealMatrix>(value, jacobian);
}
};
LeastSquaresProblem problem = new LeastSquaresBuilder()
.start(start)
.model(distancesModel)
.target(target)
.lazyEvaluation(false)
.maxEvaluations(1000)
.maxIterations(1000)
.build();
LeastSquaresOptimizer optimizer = new LevenbergMarquardtOptimizer().
withCostRelativeTolerance(1.0e-12).
withParameterRelativeTolerance(1.0e-12);
LeastSquaresOptimizer.Optimum optimum = optimizer.optimize(problem);
RealVector point = optimum.getPoint();
System.out.println("Start = " + Arrays.toString(start));
System.out.println("Solve = " + point);
}
P.S。重量的逻辑对我来说似乎很可疑。在您引用的问题中，OP 对半径有一些估计，然后它是一个明显的重量。使用以 logarithmic dBm 测量的信号强度的反向平方对我来说似乎很奇怪。
我试着模仿这个algorithm, developed in Python, that calculates geolocation based on seen Wifi stations positions, itself based on this idea。
该算法首先使用 Numpy 函数来计算观测到的纬度和经度的基本加权平均值。为了尽量减少可能的 Wifi 位置错误的影响，它还使用“scipy.optimize.leastsq”方法以统计方式计算，如果可能的话，更精确的位置。
我想在 Java Android 平台上实现相同的行为。
对于所有其他计算，我成功地依赖 org.apache.commons.math3。 所以对于最小二乘问题，我逻辑上尝试依赖 https://commons.apache.org/proper/commons-math/userguide/leastsquares.html.
我的问题，如果我很好理解的话，是 Scipy 为我管理了雅可比函数定义的复杂性，而我糟糕的数学技能不允许我正确定义最小二乘问题的模型。我尝试了一些基于此 example 的实验，这似乎接近我需要的东西，但结果并不好，因为我不知道如何处理 "jacobian" 部分。
有人为thispost做了同样的事情，有人可以为我做同样的事情并尝试以简单的方式解释吗？
有关 Python 部分工作原理的更多详细信息：
使用的“scipy.optimize.leastsq”语句是：
(lat, lon), cov_x, info, mesg, ier =
scipy.optimize.leastsq(func, initial, args=data, full_output=True)
其中数据为：latitude/longitude/age中的milliseconds/signal强度，例如：data = numpy.array([(43.48932915, 1.66561772, 1000, -20), (43.48849093, 1.6648176, 2000, -10), (43.48818612, 1.66615113, 3000, -50)])
初始计算的是加权平均值latitude/longitude，本例中：initial = 43.48864654, 1.66550075
函数是
def func(initial, data):
return numpy.array([
geographic distance((float(point[latitude]), float(point[longitude])), (initial[latitude], initial[longitude])).meters * min(math.sqrt(2000.0 / float(point[age])), 1.0) / math.pow(float(point[signal strength]), 2)
结果是：43.4885401095, 1.6648660983
我在 Java 中的实验，我替换了数据值并更改了 "modelI" 的计算方式。我简化了信号强度和年龄值。但事实证明，这还不够。
double modelI = calculateVincentyDistance(o.getY(), o.getX(), center.getY(), center.getX())* Math.min(Math.sqrt(2000.0/1000.0), 1.0) / Math.pow(-10, 2);
我也打算试试https://github.com/odinsbane/least-squares-in-java，但我不确定是否正确使用它，因为我不掌握它的工作方式。
仅供参考，我使用 Vincenty 距离计算，例如可以用 Haversine 或 Euclidean 代替。
感谢您的帮助！
代码不易移植，因为 SciPy 提供了更通用的 Least-squares minimization interface while Apache Commons Math provides curve fitting。还有许多优化问题可以重述为曲线拟合。在 Python 代码中最小化
F(current_point) = Sum{ (distance(known_point[i], current_point) * weight[i])^2 } -> min
Java曲线拟合问题有点不同：
F(current_point) = Sum{ (target_value[i] - model[i](current_point))^2 } -> min
因此，可以通过将所有 target_value 分配给 0 并使 model[i] 计算从 current_point 到 known_point[i] 的加权距离来创建等效拟合问题。
在一般情况下，此类问题无法使用公式进行精确求解，而会使用一些数值优化方法。这里还有另一个区别：Java 实现明确要求您为优化器提供方法来计算被优化函数的导数。如果未提供 Dfun，Python 代码似乎使用了某种差异区分符。您可以在 Java 中手动或使用 FiniteDifferencesDifferentiator but for simple formulas it might be easier to code them explicitly using DerivativeStructure
static class PositionInfo {
public final double latitude;
public final double longitude;
public final int ageMs;
public final int strength;
public PositionInfo(double latitude, double longitude, int ageMs, int strength) {
this.latitude = latitude;
this.longitude = longitude;
this.ageMs = ageMs;
this.strength = strength;
}
public double getWeight() {
return Math.min(1.0, Math.sqrt(2000.0 / ageMs)) / (strength * strength);
}
}
static DerivativeStructure getWeightedEuclideanDistance(double tgtLat, double tgtLong, PositionInfo knownPos) {
DerivativeStructure varLat = new DerivativeStructure(2, 1, 0, tgtLat); // latitude is 0-th variable of 2 for derivatives up to 1
DerivativeStructure varLong = new DerivativeStructure(2, 1, 1, tgtLong); // longitude is 1-st variable of 2 for derivatives up to 1
DerivativeStructure latDif = varLat.subtract(knownPos.latitude);
DerivativeStructure longDif = varLong.subtract(knownPos.longitude);
DerivativeStructure latDif2 = latDif.pow(2);
DerivativeStructure longDif2 = longDif.pow(2);
DerivativeStructure dist2 = latDif2.add(longDif2);
DerivativeStructure dist = dist2.sqrt();
return dist.multiply(knownPos.getWeight());
}
// as in https://en.wikipedia.org/wiki/Haversine_formula
static DerivativeStructure getWeightedHaversineDistance(double tgtLat, double tgtLong, PositionInfo knownPos) {
DerivativeStructure varLat = new DerivativeStructure(2, 1, 0, tgtLat);
DerivativeStructure varLong = new DerivativeStructure(2, 1, 1, tgtLong);
DerivativeStructure varLatRad = varLat.toRadians();
DerivativeStructure varLongRad = varLong.toRadians();
DerivativeStructure latDifRad2 = varLat.subtract(knownPos.latitude).toRadians().divide(2);
DerivativeStructure longDifRad2 = varLong.subtract(knownPos.longitude).toRadians().divide(2);
DerivativeStructure sinLat2 = latDifRad2.sin().pow(2);
DerivativeStructure sinLong2 = longDifRad2.sin().pow(2);
DerivativeStructure summand2 = varLatRad.cos().multiply(varLongRad.cos()).multiply(sinLong2);
DerivativeStructure sum = sinLat2.add(summand2);
DerivativeStructure dist = sum.sqrt().asin();
return dist.multiply(knownPos.getWeight());
}
使用这样的准备你可以做这样的事情：
public static void main(String[] args) {
// latitude/longitude/age in milliseconds/signal strength
final PositionInfo[] data = new PositionInfo[]{
new PositionInfo(43.48932915, 1.66561772, 1000, -20),
new PositionInfo(43.48849093, 1.6648176, 2000, -10),
new PositionInfo(43.48818612, 1.66615113, 3000, -50)
};
double[] target = new double[data.length];
Arrays.fill(target, 0.0);
double[] start = new double[2];
for (PositionInfo row : data) {
start[0] += row.latitude;
start[1] += row.longitude;
}
start[0] /= data.length;
start[1] /= data.length;
MultivariateJacobianFunction distancesModel = new MultivariateJacobianFunction() {
@Override
public Pair<RealVector, RealMatrix> value(final RealVector point) {
double tgtLat = point.getEntry(0);
double tgtLong = point.getEntry(1);
RealVector value = new ArrayRealVector(data.length);
RealMatrix jacobian = new Array2DRowRealMatrix(data.length, 2);
for (int i = 0; i < data.length; i++) {
DerivativeStructure distance = getWeightedEuclideanDistance(tgtLat, tgtLong, data[i]);
//DerivativeStructure distance = getWeightedHaversineDistance(tgtLat, tgtLong, data[i]);
value.setEntry(i, distance.getValue());
jacobian.setEntry(i, 0, distance.getPartialDerivative(1, 0));
jacobian.setEntry(i, 1, distance.getPartialDerivative(0, 1));
}
return new Pair<RealVector, RealMatrix>(value, jacobian);
}
};
LeastSquaresProblem problem = new LeastSquaresBuilder()
.start(start)
.model(distancesModel)
.target(target)
.lazyEvaluation(false)
.maxEvaluations(1000)
.maxIterations(1000)
.build();
LeastSquaresOptimizer optimizer = new LevenbergMarquardtOptimizer().
withCostRelativeTolerance(1.0e-12).
withParameterRelativeTolerance(1.0e-12);
LeastSquaresOptimizer.Optimum optimum = optimizer.optimize(problem);
RealVector point = optimum.getPoint();
System.out.println("Start = " + Arrays.toString(start));
System.out.println("Solve = " + point);
}
P.S。重量的逻辑对我来说似乎很可疑。在您引用的问题中，OP 对半径有一些估计，然后它是一个明显的重量。使用以 logarithmic dBm 测量的信号强度的反向平方对我来说似乎很奇怪。