Cómo aplicar el método de los mínimos cuadrados no lineales

Question

Estoy tratando de implementar el ejemplo de estimación de un ángulo entre un objetivo $\textbf{x}$ y un sensor $x_{p}$. Estoy utilizando el ejemplo en este libro. Hay tres mediciones disponibles del ángulo en radianes de la siguiente manera

$$ z = \begin{bmatrix} z(1) \\ z(2) \\ z(3) \end{bmatrix} = \textbf{h}(\textbf{x}, x_{p}) + w $$

donde $\textbf{x} = [x \ y]$ es la posición del objetivo, $x_{p} = [x_{px} \ x_{py}]$ es la posición del sensor y $w$ es el ruido. La función real de este sensor es

$$ \textbf{h}(\textbf{x}, x_{p}) = tan^{-1} \frac{y - x_{py}}{x - x_{px}} $$

La versión linealizada de $\textbf{h}(.)$ se proporciona a través de la Jacobiana $J$ de la siguiente manera

$$ J = \frac{\partial h(\textbf{x},x_{p})}{\partial \textbf{x}} = \begin{bmatrix} \frac{\partial h(\textbf{x},x_{p_{1}})}{\partial x} & \frac{\partial h(\textbf{x},x_{p_{1}})}{\partial y} \\ \frac{\partial h(\textbf{x},x_{p_{2}})}{\partial x} & \frac{\partial h(\textbf{x},x_{p_{2}})}{\partial y} \\ \frac{\partial h(\textbf{x},x_{p_{3}})}{\partial x} & \frac{\partial h(\textbf{x},x_{p_{3}})}{\partial y} \\ \end{bmatrix} $$

suponiendo que la varianza del ruido es la misma, la fórmula de Mínimos Cuadrados es

$$ \hat{x} = (J^{T}J)^{-1} J^{T} z $$

En mi simulación,

z =
    0,7548
    0,9402
    2,2379
J =
   -0,0621    0,0552
   -0,0552    0,0621
    0,2500    0,2500
xbar =  <--- incorrecto
   -2,7493
   11,7078

Como puedes ver, $\hat{x}$ no solo es incorrecto, sino que su dimensión es extraña. Espero que sea 3x1 pero obtengo 2x1

---

Editado: La siguiente imagen ilustra el problema.

Answer 1

En primer lugar, vamos a repasar algunos conceptos de estimación de LS (Mínimos Cuadrados). Deje que $K$ sea un entero positivo y que $\mathbf{1}_K = [1 \ 1 \ \cdots \ 1]^T \in \mathbb{R}^{K \times 1}$ sea un vector de unos ($K$-dimensional), donde $^T$ denota la transposición. Considere un único sensor haciendo $K$ observaciones (ruidosas) del mismo ángulo, $\phi$. Luego, denotando con $\mathbf{y} \in \mathbb{R}^{K \times 1}$ las observaciones del verdadero ángulo $\phi$, se tiene que $$ \mathbf{y} = \phi \mathbf{1}_K + \mathbf{w}, \quad \mathbf{w} \sim (\mathbf{0}, \ \sigma^2 \mathbf{I}_K) $$ donde $\mathbf{I}_K$ es la matriz identidad de $K \times K$ y $\mathbf{w} \in \mathbb{R}^{K \times 1}$ es un término ruidoso (del cual no nos importa su distribución específica), pero suponemos que tiene media cero, no correlacionado y con la misma varianza. Con la única excepción de la hipótesis de media cero, las otras hipótesis, para LS, no son cruciales (solo por simplicidad). La estimación "mejor" (en el sentido de LS) del ángulo $\phi$ es dada por $$ (*) \quad \hat{\phi}_{LS} = (\mathbf{1}_K^T \mathbf{1}_K)^{-1} \mathbf{1}_K^T \mathbf{y} $$ Bajo la hipótesis adicional de que $\mathbf{w} \sim \mathcal{N}(\mathbf{0}, \sigma^2 \mathbf{I}_K)$, lo anterior $(*)$ es también la estimación de ML (Máxima Verosimilitud) de $\phi$. La ecuación $(*)$ es el famoso estimador de media muestral, ya que puede ser reescrito como $$ \hat{\phi} = \frac{1}{K} \sum_{k=1}^K y_k $$ donde $y_k$ es el $k$-ésimo componente del vector $\mathbf{y}$. La media muestral se puede utilizar incluso si se tiene una matriz de covarianza desconocida $\mathbf{C$: solo la hipótesis de media cero (que es el escenario típico con AOA, Ángulo de Arribo) tiene relevancia (si la media no es cero, pero se conoce, se puede restar y reescribir el problema con el ruido de media cero; si la media no es cero y es desconocida, no se puede hacer mucho, los resultados estarán sesgados). Si se tienen $N$ sensores, entonces la posición del objetivo se puede estimar usando triangulación.

Ahora, este no es el caso. La fórmula que utiliza la matriz Jacobiana está estimando directamente la posición del objetivo: sin embargo, algunas cosas deben ser claras.

En el vector $\mathbf{z} \in \mathbb{R}^{N \times 1}$, para que lo que escribió tenga sentido, realmente debería recolectar los $N$ ángulos estimados $\hat{\phi}_n$, con $n=1,\dots, N$. Cada estimación de ángulo es una media muestral realizada en $K$ observaciones del mismo ángulo, es decir, $$ \mathbf{y}_n = \phi_n \mathbf{1}_K + \mathbf{w}_n, \quad n=1,\dots,N $$ y luego, definiendo $\mathbf{z} \triangleq [\hat{\phi}_1 \ \hat{\phi_2} \cdots \hat{\phi}_N]^T \in \mathbb{R}^{N \times 1}$, tiene sentido decir que $$ \hat{\mathbf{x}} = (\mathbf{J}^T \mathbf{J})^{-1} \mathbf{J}^T \mathbf{z} $$ donde $\hat{\mathbf{x}} \in \mathbb{R}^{2 \times 1}$ es la posición estimada del objetivo (y $\mathbf{J} \in \mathbb{R}^{N \times 2}$ es la matriz Jacobiana, obviamente). Una comprobación dimensional simple demuestra que $\hat{\mathbf{x}}$ es un vector de $2 \times 1$ (ya que el problema está en 2D).

¿Y los errores? Los errores provienen de una combinación de diversos factores. En primer lugar, la varianza del ruido, $\sigma^2$, es realmente un factor limitante en el rendimiento: el impacto de $\sigma^2$ puede reducirse aumentando $K$ (si es posible, siempre es una buena idea). En segundo lugar, el procedimiento de LS utiliza sólo la primera derivada de la función (de toda la serie de Taylor) y por lo tanto es una aproximación. Ciertamente, aumentando $N$ (es decir, más sensores), el rendimiento será (estadísticamente) mejor, pero las aproximaciones siguen siendo aproximaciones: se puede verificar "cuánto" su aproximación es problemática no generando ruido en absoluto, dejando que el algoritmo haga su trabajo y luego graficando los errores (o la ECDF, Función de Distribución Acumulada Empírica, u otra métrica que se adapte a su caso, por ejemplo, error RMS, con RMS siendo Root Mean Square).

El límite absoluto de los rendimientos de cualquier estimador está determinado por el Límite Inferior de Cramér-Rao (CRLB), que se puede calcular bajo suposiciones estándar. Se puede encontrar más información sobre el CRLB en la literatura.