Error estándar de covarianza estimada

Question

Dejar $X_1,...,X_n$ y $Y_1,...,Y_n$ser dos muestras aleatorias independientes de$\mathcal{N}(\mu, \sigma^2)$ donde ambos $\mu$ y $\sigma$ son parámetros desconocidos.

Calculo su covarianza usando: $$\hat{\operatorname{cov}}(X, Y) = \operatorname{E}{\big[(X_i - \operatorname{E}[X])(Y_i - \operatorname{E}[Y])\big]} $$

con reemplazar $\operatorname{E}[X]$ y $\operatorname{E}[Y]$ por la media muestral correspondiente.

¿Cómo calculo el error estándar de $\hat{\operatorname{cov}}(X, Y)$?

---

Editar: El comentario de Adam L. Taylor es válido. Si se sabe que ambas variables son independientes, no sería necesario estimar su covarianza; por lo que esta suposición debe descartarse.

Answer 1

La independencia de $X$ y $Y$hace que este problema sea sencillo. Para facilitar la notación, suponga$\mu=0$. Luego

$Cov (X,Y)=S_{XY}=E[(X-\mu_X)(Y-\mu_Y)]=E[XY]$,

y el estimador $\hat S_{XY} = \frac{1}{n}\sum_{i=1}^n x_iy_i$ tiene expectativa cero, entonces $Var (\hat S_{XY}) = E[\hat S_{XY}^2]$.

\begin{align} Var(\hat S_{XY}) &= E[\hat S_{XY}^2] \\ &= E\left[\left(\frac{1}{n}\sum_{i=1}^n x_iy_i\right)^2\right] \\ &= \frac{1}{n^2}E\left[\sum_{i=1}^n x_i^2y_i^2+2\sum_{i<j}x_iy_ix_jy_j\right] \\ &= \frac{1}{n^2}nE\left[X^2Y^2\right] + 0 \\ &= \frac{1}{n}E[X^2]E[Y^2] \\ &= \frac{1}{n} \sigma^4 \end{align}

Entonces, el error estándar de $\hat S_{XY} = \sqrt{Var(\hat S_{XY})}=\sigma^2/\sqrt{n}$.

Answer 2

Eso es interesante :-)

Ahora, si entendí bien tu pregunta, entonces el truco es pensar en términos de funciones y no enfocarte solo en COV. Sin embargo, las siguientes cosas son inicialmente importantes:

• Es una distribución normal.
• La covarianza es (solo) una función como cualquier otra.

Si necesita la desviación estándar para el resultado de la función COV, automáticamente asume que:

• Los elementos $X_i, Y_i$ puede tener desviaciones estándar asociadas a estos valores de medición.
• La desviación estándar $E[X], E[Y]$ ya se conoce y se puede calcular mediante la media ($X$),media($Y$). ¿Tiene alguna razón para no confiar en ella asumiendo un valor diferente?

Ahora, dada esta información, puede usar el método estándar para calcular la estimación del error resultante de una función con Propagación de incertidumbre gaussiana . Importante es la limitación; Funciona solo para variables de distribución normal. La varianza$\sigma_y^2$ de una variable $y$ que constan de otras variables inciertas $x$ y sus correspondientes variaciones $\sigma_x^2$, como

$y = x_1 + x_2 + ...+ x_n$

$\sigma_{y}^2 = \sigma_{x_1}^2+\sigma_{x_2}^2+...+\sigma_{x_n}^2$

se puede calcular como (notación matricial):

$\sigma_y^2 = \mathbf{A\Sigma}_{xx}\mathbf{A}^\mathrm{T}$.

dónde $\mathbf{A}$es la matriz jacobiana y$\mathbf{\Sigma}_{xx}$es la matriz de varianza-covarianza para los valores$X_i,Y_i$correspondiente a la función. En la diagonal, debe colocar las variaciones para$X_i,Y_i$los valores fuera de la diagonal son covarianzas entre ellos (es posible que desee asumir 0 para ellos). Tenga en cuenta que esta es una solución general para funciones no lineales y utiliza solo un (el primer) término de linealización. Es rápido y generalmente es el camino a seguir en aplicaciones productivas, pero puede tener errores de aproximación en comparación con una solución analítica pura.

Otra opción es hacer una pequeña simulación de Montecarlo. Para lograr esto, puede probar$X_i,Y_i$con su incertidumbre esperada y calcular su covarianza. Ahora bien, si lo hace varias (miles) de veces, obtiene una estimación justa del error resultante. Aquí hay un pseudocódigo para OCTAVE / MATLAB:

% Clean stuff before start to avoid variable conflicts
clc
clear all

% These are the values
X = [ 1 2 3 4 5 ].';
Y = [ 5 4 3 2 1 ].';

% How many tries do you want to have
n_samples = 10000;

% prepare the resulting error
cov_res = zeros( n_samples , 1 );

% loop the computation through n_samples
for i = 1 : n_samples
   
    % generate random distributed noise, 1 sigma [-0.1:0.1]
    x_error_sample = 0.1 * randn( size( X , 1 ) , 1 );
    y_error_sample = 0.1 * randn( size( Y , 1 ) , 1 );
    
    % Compute the covariance matrix for X and Y
    cov_i = cov( X + x_error_sample , Y + y_error_sample );
    
    % Pick only the covariance
    cov_res( i ) = cov_i( 1 , 2 );
end

% covariance estimator can be chosen by your own metric (e.g. mean,median,...)
mean( cov_res )

% The error of this estimation can be chosen by your own metric (e.g. std,rms,var,...)
std( cov_res )

Este enfoque también se puede utilizar para cualquier distribución de $X$ y $Y$, simplemente reemplace el término randn con su elección.

Saludos

Answer 3

Además de la solución analítica de @ abstrusiosity, puede aplicar el bootstrap aquí. Esto tiene la ventaja de funcionar incluso cuando las dos muestras no son independientes (la verdadera covarianza no es$0$)

library(tidyverse)

# Simulate data
m = 0
s = 2
n = 100
X = rnorm(n, m, s)
Y = rnorm(n, m, s)

(expected_se = (s**2) / sqrt(n)) # @abstrusiosity's solution
# [1] 0.2828427

nboot = 200
bootstrap_cov = map_dbl(1:nboot, function(i){
  rx = sample(X, n, replace = T)
  ry = sample(X, n, replace = T)
  cov(rx, ry)
})
sd(bootstrap_cov)
# [1] 0.270266

También puede probar que la estimación de arranque converge con la solución analítica:

sim_bootstrap = function(i){
  # Simulate fresh data
  X = rnorm(n, m, s)
  Y = rnorm(n, m, s)
  bootstrap_cov = map_dbl(1:nboot, function(i){
    # Do bootstrap sample
    rx = sample(X, n, replace = T)
    ry = sample(X, n, replace = T)
    cov(rx, ry)
  })
  sd(bootstrap_cov)
}
results = map_dbl(1:200, sim_bootstrap)
mean(results)
# [1] 0.2817101

qplot(results) + geom_vline(xintercept=expected_se, color='red') +
  labs(x='Standard Error', caption='Red line shows analytic SE')