fminsearch in R is worse than in Matlab
fminsearch, matlab, optimization, r
Solution
The problem here is that the likelihood surface is very badly behaved -- there are both multiple minima and discontinuous jumps -- which will make the results you get with different optimizers almost arbitrary. I will admit that MATLAB's optimizers are remarkably robust, but I would say that it's pretty much a matter of chance (and where you start) whether an optimizer will get to the global minimum for this case, unless you use some form of stochastic global optimization such as simulated annealing.
I chose to use R's built-in optimizer (which uses Nelder-Mead by default) rather than `fminsearch` from the `pracma` package.
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
This is better than the `fminsearch` results, but still different from the MATLAB results, and worse than them:
## 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)
The `bbmle` package offers an experimental/not very well documented set of tools for exploring optimization surfaces:
library(bbmle)
ss <- slice2D(fun=spline_fit,opt1$par,nt=51)
library(lattice)
A 2D "slice" around the `optim`-estimated parameters. The circles show the optim fit (solid) and the minimum value within each slice (open).
png("splom1.png")
print(splom(ss))
dev.off()
A 'slice' between the matlab and optim fits shows that the surface is quite rugged:
ss2 <- bbmle:::slicetrans(matlab_fit,opt1$par,spline_fit)
png("slice1.png")
print(plot(ss2))
dev.off()
Problem
There is my data (x and y columns are relevant): https://www.dropbox.com/s/b61a7enhoa0p57p/Simple1.csv What I need is to fit the data with the polyline. Matlab code that does this is: ``` 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) ``` And this does work fine, producing following optimization: [102.9842, 191.0006, 421.9912] I've implemented the same idea in 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 } ``` But I get very bad results: ``` > 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 ``` As you can see, it stops after 2 iterations and doesn't produce good points. I'll be very glad for any help in resolving this issue. Also, if anyone knows how to implement this in C# using any free library, it will be even better. I know whereto get polyfit, but not fminsearch.