3

有我的数据（x 和 y 列是相关的）： https ://www.dropbox.com/s/b61a7enhoa0p57p/Simple1.csv

我需要的是用折线拟合数据。执行此操作的 Matlab 代码是：

spline_fit.m:
function [score, params] = spline_fit (points, x, y)

min_f = min(x)-1;
max_f = max(x);

points = [min_f points max_f];
params = zeros(length(points)-1, 2);

score = 0;
for i = 1:length(points)-1
    in = (x > points(i)) & (x <= points(i+1));
    if sum(in) > 2
        p = polyfit(x(in), y(in), 1);
        pred = p(1)*x(in) + p(2);
        score = score + norm(pred - y(in));
        params(i, :) = p;
    else
       params(i, :) = nan;
    end
end


test.m:
%Find the parameters
r = [100,250,400];
p = fminsearch('spline_fit', r, [], x, y)
[score, param] = spline_fit(p, x, y)

%Plot the result
y1 = zeros(size(x));
p1 = [-inf, p, inf];
for i = 1:size(param, 1)
    in = (x > p1(i)) & (x <= p1(i+1));
    y1(in) = x(in)*param(i,1) + param(i,2);
end

[x1, I] = sort(x);
y1 = y1(I);

plot(x,y,'x',x1,y1,'k','LineWidth', 2)

这确实可以正常工作，产生以下优化：[102.9842, 191.0006, 421.9912]

我在 R 中实现了同样的想法：

library(pracma);
spline_fit <- function(x, xx, yy) {

  min_f = min(xx)-1;
  max_f = max(xx);

  points = c(min_f, x, max_f)
  params = array(0, c(length(points)-1, 2));

  score = 0;
  for( i in 1:length(points)-1)
  {
    inn <- (xx > points[i]) & (xx <= points[i+1]);
    if (sum(inn) > 2)
    {
      p <- polyfit(xx[inn], yy[inn], 1);
      pred <- p[1]*xx[inn] + p[2];
      score <- score + norm(as.matrix(pred - yy[inn]),"F");
      params[i,] <- p;
    }
    else
      params[i,] <- NA;
  }  
  score
}

但我得到了非常糟糕的结果：

> fminsearch(spline_fit,c(100,250,400), xx = Simple1$x, yy = Simple1$y)
$xval
[1] 100.1667 250.0000 400.0000

$fval
[1] 4452.761

$niter
[1] 2

如您所见，它在 2 次迭代后停止并且不会产生好点。

对于解决此问题的任何帮助，我将非常高兴。

此外，如果有人知道如何使用任何免费库在 C# 中实现这一点，那就更好了。我知道从哪里获得 polyfit，但不知道 fminsearch。

4

1 回答 1

17

这里的问题是似然面的表现非常糟糕——同时存在多个最小值和不连续的跳跃——这将使您使用不同的优化器获得的结果几乎是任意的。我承认 MATLAB 的优化器非常健壮，但我想说，对于这种情况，优化器是否会达到全局最小值几乎是一个机会问题（以及你从哪里开始），除非你使用某种形式的随机全局优化比如模拟退火。

我选择使用 R 的内置优化器（默认使用 Nelder-Mead）而不是包中fminsearch的。pracma

spline_fit <- function(x, xx = Simple1$x, yy=Simple1$y) {

    min_f = min(xx)-1
    max_f = max(xx)

    points = c(min_f, x, max_f)
    params = array(0, c(length(points)-1, 2))

    score = 0
    for( i in 1:(length(points)-1))
    {
        inn <- (xx > points[i]) & (xx <= points[i+1]);
        if (sum(inn) > 2)
        {
            p <- polyfit(xx[inn], yy[inn], 1);
            pred <- p[1]*xx[inn] + p[2];
            score <- score + norm(as.matrix(pred - yy[inn]),"F");
            params[i,] <- p;
        }
        else
            params[i,] <- NA;
    }  
    score
}

library(pracma) ## for polyfit
Simple1 <- read.csv("Simple1.csv")
opt1 <- optim(fn=spline_fit,c(100,250,400), xx = Simple1$x, yy = Simple1$y)
## [1] 102.4365 201.5835 422.2503

这比fminsearch结果好，但仍然与 MATLAB 结果不同，并且比它们差：

## Matlab results:
matlab_fit <- c(102.9842, 191.0006, 421.9912)
spline_fit(matlab_fit, xx = Simple1$x, yy = Simple1$y)
## 3724.3
opt1$val
## 3755.5  (worse)

该bbmle软件包提供了一组用于探索优化表面的实验性/未充分记录的工具：

library(bbmle)
ss <- slice2D(fun=spline_fit,opt1$par,nt=51)
library(lattice)

optim估计参数周围的二维“切片” 。圆圈显示最佳拟合（实心）和每个切片内的最小值（打开）。

png("splom1.png")
print(splom(ss))
dev.off()

在此处输入图像描述

matlab 和 optim 拟合之间的“切片”表明表面非常粗糙：

ss2 <- bbmle:::slicetrans(matlab_fit,opt1$par,spline_fit)
png("slice1.png")
print(plot(ss2))
dev.off()

在此处输入图像描述

于 2013-11-17T19:38:01.957 回答