How can I estimate unknown parameters in nonlinear forward model?

Question

I have a table in $Z=0$. Above the table in $(YZ,ZC)$ I have an instrument, C. Think of the instrument as a monocular rotating around the X axis. With the instrument I can measure at which angle, $\theta$, I observe a target (red dot) on the table. I have approximate measurements of $(YZ,ZC)$ and want to estimate these more accurately. Also the instrument has an initial rotation of r that I want to estimate more closely.

From the figure we see that: $$ ZC*tan(\theta + r) = (YL - YC) $$

I can place as many targets on the table, N, as I want and measure their YL. However measuring YL and identifying the target in the monocular is a bit of work so I prefer to have N as small as possible.

Also the angle $\theta$ is restricted to: $-10*\pi/180 <= \theta <= 10*\pi/180$.

Let us assume that $$ YC = -0.25, ZC = 1, r = 4*\pi/180 $$

Here is some code to generate N pairs of $(\theta,YL)$ which constitute the "measurements"/data:

import math
import random
import numpy as np

YC = -0.25
ZC = 1
r = 4
theta = []
YL = []
def printThetaYL(N,mu, sigma):
    thetaVec = np.linspace(-10,10,N)

    for theta_i in thetaVec:
          YL_i = ZC*math.tan(math.radians(r + theta_i)) + YC
          noise = np.random.normal(mu, sigma)
          YL_i = YL_i + noise
          print("(theta={:.3f},YL={:.3f})".format(theta_i,YL_i))
          theta.append(theta_i)
          YL.append(YL_i)

printThetaYL(5, 0, 0.0005)

For initial estimate I want to use: $$ YC \approx -0.2, ZC \approx 0.95, r \approx 3*\pi/180 $$

How can I estimate the parameters $(YZ,ZC,r)$ accurately using Python with as few target points, N, as possible?

I tried by introducing the error/loss function: $$ SE = \frac{1}{N}\sum_{i=1}^{N} [YC + ZC*tan(\theta_i + r) - YLi]^2 $$

# Loss function SE=f(YC,ZC,r):
def SE(x):
    r  = x[0]
    YC = x[1]
    ZC = x[2]
    N = len(theta)
    E = 0
    for i in range(0,N):
        theta_i = theta[i]
        YL_i = YL[i]
        Ei = YC + ZC * math.tan(math.radians(r + theta_i)) - YL_i
        E = E + Ei**2
    return E/N

And minimize this:

res = scipy.optimize.minimize(SE, x0=[3.0, -0.2, 0.95])

However the results were not so great. For N=20 I got the estimates: $$ r=3.000,YC=-0.232,ZC=1.002 $$ The results for ZC is good (2 mm off), but the results for $\theta$ (1 degree off) and YC (16.8 mm off) are not good. Interestingly I got the same results without noise ($\sigma=0$). Is there some inherit problem with the way my "experiment" is defined?

Answer 1

import numpy as np
from scipy.optimize import minimize


YC = -0.25
ZC = 1
r  = 4*np.pi/180
N  = 20;
thetamin =-10*np.pi/180
thetamax = 10*np.pi/180


def printThetaYL(mu, sigma):
    theta = np.linspace(thetamin,thetamax,N)
    data = []
    for theta_i in theta:
          YL_i = ZC*np.tan(r + theta_i) + YC
          noise = np.random.normal(mu, sigma)
          YL_i = YL_i + noise
          data.append([theta_i,YL_i])
    return data

data = printThetaYL(0, 0.0005)

def obj(X):
    (ZC,YC,r) = X
    s = 0
    for sample in data:
        (theta,YL) = sample
        res = ZC*np.tan(theta + r) + YC - YL
        s += res*res
    return s

x0 = [0,0,0]
bds = ((0,None),(None,None),(-np.pi/2,np.pi/2))
sol = minimize(obj, x0, bounds = bds)

print("obj = {:.5f}".format(sol.fun))
(ZC,YC,r) = sol.x
print("ZC = {:.4f}, YC = {:.4f}, r = {:.4f}".format(ZC,YC,r))