I’m trying to solve a simple parameter estimation problem in Julia similar to how fminsearch works in Octave/Matlab. The RAFF.jl package seemed like a nice approach, but perhaps there’s a better way?
For context, I have a photo of several lampposts in a public square that I know are all the same height (H), but I don’t know what the height is, and I know nothing about the camera, but I’m assuming that each pixel extends over a small angle such that there are r pixels per radian. I know the distance (Dst) from the camera to each lamppost so the problem is find r and H that minimize the difference between the tangent of the angle and H/Dst,
In Octave (n,Dst defined below):
f = @(x) sum((tan(n/x(1)) - x(2)./Dst).^2);
C = fminsearch(f,[1500,12]) = [1503.360 12.900]
The solutions generated by RAFF don’t match the Octave solution, so I’m guessing that either I’m using the wrong package, or not setting RAFF up properly. Here are the results:
# Data inputs to the problem
n = [765.0, 332.0, 217.0, 165.0, 283.0, 297.0, 814.0, 434.0, 289.0];
Dst = [22.71666667, 45.64, 71.58333333, 90.22666667, 59.63, 59.33666667, 23.82333333, 46.88333333, 72.88];
# Write data in format useful for RAFF.jl
x = [n Dst zeros(9,1)];
# Function to be minimized over parameters p
f(x,p) = (tan(x[1]/p[1]) - p[2]/x[2])^2
# Jacobian
begin
gmodel!(g,x,p) = begin
g[1] = 2*x[1]/p[1]^2*sec(x[1]/p[1])^2*(tan(x[1]/p[1]) - p[2]/x[2])
g[2] = 2/p[2]*(tan(x[1]/p[1]) - p[2]/x[2])
end
end
# Fit using RAFF
fit = raff(f,gmodel!,x,2; initguess = [1000.0, 10.0], ε=1.0e-10)
# Output (solution only)
Solution (.solution) = [1000.0, 10.0]
# Without Jacobian
fit = raff(f,x,2; initguess = [1000.0, 10.0], ε=1.0e-10)
Solution (.solution) = [56821.055712178044, 0.31544255532138893]
# No initial guess
fit = raff(f,x,2)
Solution (.solution) = [-0.9778814405815958, -3.5093908435870933]
# Function re-written to use Dst as first column in x, n as second column (no zeros)
f1(x,p) = p[1]*atan(p[2]/x[1])
Solution (.solution) = [477.50199698046083, 39.851464869409654]
# Function re-written to use n as first column in x, Dst as second column
f2(x,p) = p[2]/(tan(x[1]/p[1]))
Solution (.solution) = [-0.8529985482107982, 5.582472558856189]