Estadística con múltiples variables e introducción al lenguaje R

Algunos Conceptos Básicos de Estadística

Introducción

En este capítulo se presentarán algunos conceptos de estadística básica de una variable que serán necesarios para el desarrollo de métodos en varias variables. Se presentarán algunos conceptos de inferencia estadística, junto con algunas distribuciones elementales. Suponemos que el lector ha tenido algún contacto con los conceptos de estadística en una variable como parte de algún curso a nivel universitario.

Variables Aleatorias

Los modelos matemáticos se dividen en tres clases generales: (1) Puramente deterministas, (2) estáticos o deterministas con algunos componentes aleatorios sencillos, y (3) estocásticos. Cualquier observación de un modelo determinista es estrictamente una función de sus parámetros y variables tales como el tiempo, el espacio, o energías o estímulos varios. Los modelos puramente deterministas son típicos de la física clásica, donde por ejemplo, la distancia recorrida por un cuerpo en caída libre es una función de su altura y velocidad inicial. Si se ignoran efectos transitorios como la turbulencia atmosférica, error observacional entre otros, el desplazamiento se puede calcular exactamente dados un intervalo de tiempo y la constante gravitacional.

En el segundo tipo de modelo, cada observación es una función de un componente estrictamente determinista y un término aleatorio atribuible a errores en las medidas o a las variaciones en el muestreo de las observaciones o de las variables de entrada. Los modelos examinados aquí son principalmente de este tipo con la restricción de que la componente aleatoria será simplemente agregada a la parte determinista.

Los modelos estocásticos se construyen a partir de eventos o componentes fundamentalmente aleatorios para explicar fenómenos dinámicos o evolutivos: pueden ser tan simples como una secuencia de lanzamientos de una moneda al aire o tan complejos como los procesos de nacimiento y muerte que describen la dinámica de una población biológica.

Definamos ahora con más precisión que significa el concepto de variación aleatoria o los componentes aleatorios en los modelos tipo 2 y 3 mencionados anteriormente. Empezamos definiendo una variable aleatoria discreta o una que pueda contener un número contable de valores.

Supóngase que un experimento sólo puede resultar en uno y sólo uno de k resultados llamados $E_{1},...,E_{k}$. Estos resultados son mutuamente exclusivos, en el sentido de que la ocurrencia de uno niega la existencia de cualquier otro. A cada evento se le asigna un valor numérico $p_{i}$ entre cero y uno. Este valor numérico se le conoce como la probabilidad $P(E_{i})$ de ese evento. Entonces $p_{i}$ es la probabilidad de que ocurra el evento $E_{i}$ en una corrida del experimento. Dentro del marco del experimento se le asigna la probabilidad de cero a los eventos imposibles y la probabilidad de uno a los eventos que ocurrirán con entera certidumbre. De modo que, por exclusividad mutua de los eventos, en una corrida del experimento

${\begin{matrix}P(E_{i}\cap E_{j})&=&0\\P(E_{i}\cup E_{j})&=&p_{i}+p_{j}\end{matrix}}$ donde el símbolo de intersección $\cap$ denota el evento $E_{i}$

y $E_{j}$ y el símbolo de unión $\cup$ indica el evento $E_{i}$ o $E_{j}$.

Usando la propiedad aditiva de las probabilidades de eventos mutuamente exclusivos la probabilidad total del conjunto de eventos es

$P(E_{1}\cup E_{2}\cdots \cup E_{k})=p_{1}+\cdots +p_{k}=1$

Ahora se le asigna el valor numérico $x_{i}$ al i-ésimo resultado, donde, por conveniencia, los resultados se ordenan ascendentemente según sus valores $x_{i}$. Se define entonces la variable aleatoria discreta $X$ como la cantidad que toma el valor $x_{i}$ con probabilidad de $p_{i}$ en cada corrida del experimento. Por ejemplo, si el experimento consiste en lanzamientos de una moneda, el resultado "cara" puede tener una cierta probabilidad $p_{1}$ y el resultado "cruz" puede tener una probabilidad $p_{2}$. Esta variable aleatoria puede ser descrita entonces por una función de probabilidad $f(x_{i})$ la cual especifica las probabilidades con las que la variable aleatoria $X$ adquiere los valores de cara o cruz.

$X$	$f(x_{i})$
cara	$p_{1}$
cruz	$p_{2}$

Se hace notar que siempre la suma de todas las probabilidades debe ser igual a uno y que implícitamente se le ha asignado una probabilidad de cero a eventos irrelevantes tales como la moneda cayendo al canto.

Las variables aleatorias que se verán más adelante tomarán valores sobre alguna región continua en lugar de hacerlo sobre un conjunto de eventos discretos y se les llamará variables aleatorias continuas, variantes continuas o simplemente variantes. Todos estos términos son sinónimos. La variable aleatoria continua $X$ definida sobre el dominio de los números reales se caracteriza por su función de distribución

$F(x)=P(X\leq x)\;\;-\infty <x<\infty$

la cual proporciona la probabilidad de que $X$ sea menor o igual que algun valor de $x$ en su dominio. Ya que $X$ es continua, $P(X=x)=0$.

Si $F(x)$ es una función absolutamente continua, la versión continua de la función discreta de probabilidad es la función de densidad

$f(x)={\dfrac {dF(x)}{dx}}$

Al mismo tiempo, por la propiedad de continuidad absoluta

$F(x)=\int _{-\infty }^{x}f(u)\;du$

y de esta definición se genera el término equivalente función acumulativa de distribución de $F(x)$. Nótese que éstas definiciones son perfectamente generales: si la variable aleatoria se encuentra definida en algún intervalo de la recta de los números reales, fuera de ese intervalo $f(x)$ es igual a cero, y $F(x)$ es igual a cero y a uno a la izquierda y a la derecha del intervalo, respectivamente.

Cuando se consideran con respecto a su función de densidad $f(x)$, los valores en el intervalo se definen como la población o el universo de la variable aleatoria $X$.

Las propiedades de una variable aleatoria se visualizan comúnmente en términos de su función de densidad. Los términos exponencial o rectangular se refieren a estas representaciones. La figura muestra las funciones de densidad de tres variables conocidas. Las funciones de distribución de las variables exponenciales y normales se obtienen integrando directamente. La función de distribución normal sólo se puede evaluar por integración numérica.

Las funciones de la figura involucran parametros que determinan sus posiciones y sus formas. Los parámetros ${\frac {1}{2}}(a+b)$ de la densidad rectangular y $\mu$ de la función normal son parámetros de ubicación porque sus valores especifican las posiciones de las densidades en el eje real. El intervalo $b-a$ de la variable rectangular, el parámetro exponencial $\beta$, y el parámetro $\sigma$ de la densidad normal, son parámetros de escala porque cambios en sus valores equivalen a cambios en las unidades de las variables. Los incrementos en los valores de estos parámetros implican un mayor extendimiento de la función de densidad y por lo tanto más variación en la variable aleatoria. En general, si la densidad $f(x;\alpha ,\beta )$ se puede expresar como

$f\left({\dfrac {x-\alpha }{\beta }}\right)$

entonces $\alpha$ es un parámetro de ubicación y $\beta$ es un parámetro de escala. Se hace notar que todas las funciones de densidad de esta forma contienen el factor $1/\beta$ asociado con el elemento diferencial $dx$.

Las variables aleatorias y sus densidades tienen una interpretación física: considérese la función de densidad $f(x)$ de la variable aleatoria $X$ como la función que define la densidad de una barra sólida que se encuentra posicionada en el eje de las $x$. El k-ésimo momento alrededor del origen es

$\mu _{k}^{'}=\int _{-\infty }^{\infty }x^{k}f(x)\;dx$

El primer momento es la media o el valor esperado

$E(X)=\int _{-\infty }^{\infty }xf(x)\;dx$

de la variable aleatoria y corresponde al concepto físico del centro de gravedad de la vara. El símbolo $E$ denota la operación de calcular el valor esperado, y se conoce como el operador de esperanza. Si $c$ y $k$ son cantidades no aleatorias, entonces podemos aseverar las siguientes propiedades de las esperanzas:

${\begin{aligned}E(c)&=c\\E(cX)&=cE(X)\\E(k+cX)&=k+cE(X)\end{aligned}}$

Como ejemplo, considérese la función de densidad rectangular de la figura

Figura 1 - Funciones de Densidad Rectangular y Exponencial

. Por definición:

${\begin{aligned}E(X)&={\dfrac {1}{b-a}}\int _{a}^{b}x\;dx\\&={\dfrac {1}{b-a}}{\dfrac {1}{2}}x^{2}{\Bigr |}_{a}^{b}\\&={\frac {1}{2}}(b+a)\end{aligned}}$

Se puede intuir que el valor medio es el punto medio de los límites de la función de densidad, y si se desea, se puede reparametrizar la función de densidad en términos de la media como parámetro de ubicación y el intervalo $w=b-a$ como el parámetro de escala. Entonces, la nueva densidad será:

$h(x)={\begin{cases}0,&-\infty <x<E(X)-{\frac {1}{2}}w\\{\frac {1}{w}},&E(X)-{\frac {1}{2}}w\leq x\leq E(X)+{\frac {1}{2}}w\\0,&E(X)+{\frac {1}{2}}w<x<\infty \end{cases}}$

La varianza de una variable aleatoria es el valor esperado de las desviaciones cuadradas alrededor de la media, conocido también como el segundo momento central. Se denota la varianza de $X$ como

${\text{var}}(X)$

o también como $\sigma _{x}^{2}$

Por definición

${\begin{aligned}{\text{var}}(X)=&\int _{-\infty }^{\infty }\left[x-E(X)\right]^{2}f(x)\;dx\\=&\int _{-\infty }^{\infty }x^{2}f(x)\;dx-\left[E(X)\right]^{2}\\=&E(X^{2})-\left[E(X)\right]^{2}\end{aligned}}$

y el simbolo var aplicado a cualquier variable aleatoria denotará el cálculo de la varianza. Si $c$ y $k$ son constantes, la varianza posee las siguientes propiedades:

${\begin{aligned}{\text{var}}(c)&=0\\{\text{var}}(cX)&=c^{2}{\text{var}}(X)\\{\text{var}}(X+c)&={\text{var}}(X)\end{aligned}}$

La primera propiedad simplemente afirma que una variable no aleatoria tiene cero varianza. La naturaleza cuadrática de la varianza se encuentra reflejada en la segunda propiedad ya que un cambio de escala $c$ altera la varianza en $c^{2}$. El tercer resultado proviene de la definición y afirma que la varianza no se ve afectada al cambiar el origen en el eje de las $X$. Se puede demostrar que las varianzas de las funciones de densidad rectangular, exponencial, y normal

Figura 2 - Funciones de Densidad Rectangular, Exponencial, y Gaussiana

son ${\frac {1}{12}}(b-a)^{2}$, $\beta ^{2}$, y $\sigma ^{2}$ respectivamente.

En muchas ocasiones es deseable tener una medida de dispersión que posea las mismas unidades de la variable. La desviación estándar de $X$ es la raiz cuadrada positiva de la varianza:

$\sigma _{x}=+{\sqrt {{\text{var}}(X)}}$

Hacemos notar que los parámetros naturales de la densidad normal en la Figura 2 son la media y la desviación estándar.

Variables Independientes

En las primeras partes de este capítulo nos referimos vagamente a las variables aleatorias "independientes". Ahora proporcionaremos una definición precisa de la independencia de variables. Como se verá en el Capitulo 3, se puede extender el concepto de distribución y de funciones de densidad a múltiples variables, y si se puede expresar la función de distribución conjunta de $X_{1},\ldots ,X_{p}$ como

$F(x_{1},\ldots ,x_{p})=\int _{-\infty }^{x_{p}}\cdots \int _{-\infty }^{x_{1}}f(u_{1},\ldots ,u_{p})du_{1}\cdots du_{p}$

donde $f(u_{1},\ldots ,u_{p})$ es la densidad conjunta, entonces las variables se dice que son independientes si y solo si

$F(x_{1},\ldots ,x_{p})=F_{1}(x_{1})\cdots F_{p}(x_{p})$

donde $F_{i}(x_{i})$ es la función de distribución de la variable sencilla $X_{i}$. Alternativamente, la independencia es válida si y sólo si la factorización de

$f(x_{1},\ldots ,x_{p})=f_{1}(x_{1})\cdots f_{p}(x_{p})$

de la densidad conjunta en el producto de las funciones individuales de densidad es válida.

El momento de productos

$E(X_{1}^{k_{1}}\cdots X_{p}^{k_{p}})=\int _{-\infty }^{\infty }\cdots \int _{-\infty }^{\infty }x_{1}^{k_{1}}\cdots x_{p}^{k_{p}}f(x_{1},\ldots ,x_{p})\;dx_{1}\cdots dx_{p}$

de las variables $X_{1},\ldots ,X_{p}$ se factoriza como el producto

$E(X_{1}^{k_{1}})\cdots E(X_{p}^{k_{p}})$

de los $k$ momentos individuales si las variables son independientes. De este resultado se puede concluír que la varianza de la suma de variables aleatorias independientes es simplemente la suma de las varianzas individuales: ${\text{var}}(X_{1}+\cdots +X_{p})={\text{var}}(X_{1})+\cdots +{\text{var}}(X_{p}).$

Variables Aleatorias Normales (Gaussianas)

En esta sección se describirán algunas propiedades de una sola variable aleatoria normal como preparación para los resultados y técnicas que se basarán en la distribución normal multivariable. Recuérdese que la función de densidad normal es

$f(x)={\dfrac {1}{\sqrt {2\pi \sigma }}}\exp \left[-{\dfrac {1}{2}}\left({\dfrac {x-\mu }{\sigma }}\right)^{2}\right]\;\;-\infty <x<\infty$

y que su función de distribución se encuentra dada por la integral

$P(x)={\dfrac {1}{\sqrt {2\pi \sigma }}}\int _{-\infty }^{x}\exp \left[-{\dfrac {1}{2}}\left({\dfrac {x-\mu }{\sigma }}\right)^{2}\right]dx={\dfrac {1}{\sqrt {2\pi }}}\int _{-\infty }^{(x-\mu )/\sigma }\exp(-{\frac {1}{2}}u^{2})du=\Phi \left({\dfrac {x-\mu }{\sigma }}\right)$

donde $\Phi (z)$ denota la función de distribución normal estándar con media igual a cero y varianza igual a uno. Se denotará la distribución de una variable normal por medio del llamado símbolo de Wilks: $N(\mu ,\sigma ^{2})$; la distribución estandar será indicada por $N(0,1)$.

Los valores de $\Phi (z)$ se pueden encontrar fácilmente en cualquier tabla de funciones estadísticas. Por otro lado, el porcentaje 100α de la distribución normal unitaria se define como el valor $z_{\alpha }$ tal que

$\alpha =P(Z>z_{\alpha })=1-\Phi (z_{\alpha })$

Estas relaciones entre porcentajes y sus probabilidades se ilustran en la figura 3

Fig 3-Relaciones entre funciones de densidad y probabilidades

para la densidad normal y la función de distribución.

Ahora, permítasenos ofrecer la justificación donde mucha de la metodología estadística se basa para justificar la suposición de una población normal. Existe el llamado teorema del límite central el cual afirma que las variables aleatorias que son la suma de muchos efectos independientes tienden a encontrarse normalmente distribuídas a medida que el número de efectos va creciendo. Formalmente:

Si las variables aleatorias $X_{1},\ldots ,X_{N}$ se encuentran distribuídas independientemente de acuerdo a alguna función común de distribución con media $E(X_{i})=\mu$ y varianza finita ${\text{var}}(X_{i})=\sigma ^{2}$, entonces a medida que el número de variantes $N$ aumenta sin cotas, la variante

$Z={\dfrac {\sum _{i=1}^{N}X_{i}-N_{\mu }}{\sqrt {N\sigma ^{2}}}}$

converge en su distribución a una variante normal con media cero y varianza igual a uno.

Se puede encontrar una demostración de ésta versión del teorema del límite central en cualquier texto de probabilidad o estadística matemática. Una consecuencia inmediata del teorema es que la función ${\sqrt {N}}({\bar {X}}-\mu )$ de la media de las muestras de una secuencia de variables aleatorias independientes

${\bar {X}}={\frac {1}{N}}(X_{i}+\ldots +X_{N})$

cuya distribución común tiene media $\mu$ y varianza $\sigma ^{2}$ tiende a distribuírse como una variante de tipo $N(0,\sigma ^{2})$ a medida que $N$ se vuelve más grande. Una y otra vez estaremos invocando una propiedad muy importante de la distribución normal: Si $X_{1},\ldots ,X_{n}$ son variantes independientes con distribuciones $N(\mu _{1},\sigma _{1}^{2}),\ldots ,N(\mu _{n},\sigma _{n}^{2})$ entonces la combinación lineal

$Y=a_{1}X_{1}+\ldots +a_{n}X_{n}$

también se encuentra normalmente distribuída con media $\sum _{i=1}^{n}a_{i}\mu _{i}$ y varianza $\sum _{i=1}^{n}a_{i}^{2}\sigma _{i}^{2}$. Veremos más adelante que este resultado se puede extender a una combinación lineal de variantes normales dependientes.

Nota: Muchas de las siguientes funciones hacen uso de la función gamma $\left(\Gamma (x)\right)$. En el recuadro al final de esta sección se hace una breve presentación de esta función y sus propiedades.

La distribución ji cuadrada

Muchas distribuciones se pueden derivar a partir de diferentes transformaciones de un conjunto de variantes normalmente distribuídas. Una de las más importantes es la de la variante ji cuadrada. Si las variantes $X_{1},\ldots ,X_{n}$ son independientes y se encuentran normalmente distribuídas con media cero y varianza igual a uno, entonces

$\chi ^{2}=X_{1}^{2}+\ldots +X_{n}^{2}$

tiene la función de densidad

$f(\chi ^{2})={\frac {2^{-n/2}}{\Gamma (n/2)}}(\chi ^{2})^{(n-2)/2}\exp(-{\frac {1}{2}}\chi ^{2})\ \ 0\leq \chi ^{2}<\infty$

y se dice que es una variante ji-cuadrada con $n$ grados de libertad. Debido a que comenzamos con variantes normales estśndar, la densidad no contiene parámetros de ubicación o de escala sino sólo el parámetro $n$. Aún más, si las $X_{i}$ se encuentran independientemente distribuídas como variantes $N(\mu _{i},\sigma _{i}^{2})$, entonces la cantidad

$\sum _{i=1}^{n}{\dfrac {(X_{i}-\mu _{i})^{2}}{\sigma _{i}^{2}}}$

tiene una distribucion ji-cuadrada con $n$ grados de libertad, ya que cada término cuadrado en la suma ha sido normalizado con su media y su varianza. En el apéndice mostraremos una pequeña rutina en lenguaje ${\textsf {R}}$ para obtener una tabla de los valores de la distribución ji-cuadrada con varios grados de libertad. Expresaremos el porcentaje 100α de la distribución ji-cuadrada con $n$ grados de libertad como:

$\chi _{\alpha ;n}^{2}$

donde $\alpha =P(\chi ^{2}>\chi _{\alpha ;n}^{2})$.

La distribución t

La variable aleatoria $t$ con $n$ grados de libertad se define como el cociente

$t={\dfrac {z}{\sqrt {\chi ^{2}/n}}}$

de una variante normal estándar $z$ y la raíz cuadrada de una variante ji cuadrada independiente dividida esta por sus $n$ grados de libertad. La variable $t$ es adimensional y su función de densidad

$f(t)={\dfrac {\Gamma \left[(n+1)/2\right]}{{\sqrt {\pi n}}\ \Gamma (n/2)}}{\dfrac {1}{(1+t^{2}/n)^{(n+1)/2}}}\quad -\infty <t<\infty$

depende solamente del parámetro de grados de libertad $n$. Veremos más adelante como generar valores de una distribución $t$ con $n$ grados de libertad usando el lenguaje ${\textsf {R}}$.

Se expresara el porcentaje 100α como $t_{\alpha ;n}$ o

$\alpha =P(t>t_{\alpha ;n})$

La distribución F

El cociente

$F={\dfrac {\chi _{1}^{2}/m}{\chi _{2}^{2}/n}}$

de variantes ji-cuadradas independientes divididas por sus respectivos grados de libertad, tiene la distribución $F$ o cociente de varianzas con función de densidad

$f(F)={\dfrac {\Gamma \left[(m+n)/2)\right]}{\Gamma (m/2)\Gamma (n/2)}}\left({\dfrac {m}{n}}\right)^{m/2}F^{(m-2)/2}\left(1+{\frac {m}{n}}F\right)^{-(m+n)/2}\ 0\leq F<\infty$

Se denotará el porcentaje superior 100α de la distribución $F$ con $m,n$ grados de libertad por medio de $F_{\alpha ;m,n}$:

$\alpha =P(F>F_{\alpha ;m,n})$

De la definición de la variante $F$ se deduce que los puntos porcentuales inferiores se pueden obtener de los recíprocos de los valores superiores con los grados de libertad invertidos:

$F_{1-\alpha ;m,n}={\dfrac {1}{F_{\alpha ;m,n}}}$

Una Nota Breve Sobre la Funcion $\Gamma (n)$: La función $\Gamma (n)$ se define como una extensión de la función factorial a los números complejos y reales. La relación con la función factorial es:
$\Gamma (n)=(n-1)!$
La función gamma se puede definir como una integral definida en $R[z]>0$ (La forma integral de Euler):
$\Gamma (z)\equiv \int _{0}^{\infty }t^{z-1}e^{-t}\;dt$
o también
$\Gamma (z)=\int _{0}^{\infty }e^{-t^{2}}t^{2z-1}\;dt.$

Muestras Aleatorias y Estimación

Hasta este punto sólo se ha hablado de las variables aleatorias en función de su origen. Este origen son poblaciones abstractas las cuales se encuentran especificadas por su distribución o función de densidad. De vez en cuando estas funciones se les encuentra subyaciendo en algún fenómeno aleatorio, y entonces se puede crear una descripción del fenómeno a partir de su modelo matemático. Por lo general, lo que sucede realmente es que no se conocen ni la forma matemática de la distribuición ni sus parámetros. Entonces es necesario ir más alla de la teoría de la probabilidad y llegar al dominio de la inferencia estadística para obtener estimaciones de $F(x)$ o de sus parámetros a partir de muestras finitas de valores de la variable aleatoria. En esta sección se considerará un enfoque heurístico el cual lleva a obtener estimaciones con algunas propiedades requeridas.

Empecemos suponiendo que los valores de la variable aleatoria continua $X$ pueden ser observados y registrados. Aún cuando este requisito pudiera parecer obvio, no siempre es posible observar directamente la variable de interés. Pueden interferir varios factores tales como limitaciones del equipo de medición, o el fenómeno bajo observación no expresa la variable en forma directamente medible. De cualquier modo, supondremos, para simplificar las cosas, que nuestros datos son observaciones de una población continua. Supondremos también que la expresión matemática de la función de densidad $f(x;\theta _{1},\ldots ,\theta _{k})$ de $X$ se sabe a partir de algunas suposiciones anteriores aunque no se sepan los valores de los parámetros $\theta _{i}$.

Definamos ahora el espacio de parámetros de una densidad. Supongamos que la densidad depende de un sólo parámetro real $\theta$, como en el caso de una exponencial descendente. Entonces el espacio de parámetros de $\theta$ es la porción del eje de los reales que contiene todos los valores permitidos de $\theta$. En el caso de la exponencial descendente, el espacio de parámetros sería la parte positiva del eje de los reales, ya que los parámetros negativos anularían la propiedad de densidad, y un parámetro igual a cero convierte a la función en una constante. En el caso de una distribución normal con varianza conocida el espacio sería toda la recta de los reales. Extendiendo este razonamiento, el espacio de parámetros para una densidad de $k$ parámetros sería alguna región del espacio euclideano en $k$ dimensiones. Por ejemplo, el de una distribución normal sería la mitad superior del plano cartesiano con ejes $(\mu ,\sigma ^{2})$ (Fig.4).

Fig. 4 - Ejes del Espacio de Parámetros

Finalmente, se deben definir las unidades experimentales o unidades muestrales sobre las cuales se observaran los valores de $X$. Las unidades disponibles deben contener una coleccion homogénea con respecto a todas las características que puedan afectar los valores de la variable. Por ejemplo, si la variable aleatoria de interés es el nivel de triglicéridos en mujeres mexicanas adolescentes, entonces las unidades muestrales no deben incluír hombres, mujeres adultas, o niños de cualquier sexo. Las unidades muestrales deben ser independientes unas de las otras y no deben poseer cualidades comunes que puedan crear valores dependientes de $X$. Debido a que existe una abundante literatura sobre el tema de mediciones y muestreos no se abundará mas sobre el tema.

Supóngase ahora que se han seleccionado $N$ unidades aleatoriamente. Sus valores observados de $X$ se denotaran como $x_{1},\ldots ,x_{N}$. El primer problema es estimar los parametros $\theta _{i}$ a traves de una función apropiada de las observaciones. Tales estimaciones serán denotadas por un parámetro con capucha: ${\hat {\theta _{i}}}$. Estas cantidades serán llamadas estimaciones puntuales, ya que son valores asignados de $\theta$, a diferencia de los intervalos de estimación que se verán en la siguiente seccio}ón.

El problema de la estimación de parámetros no es trivial y lleva a una serie de consideraciones y cuestionamientos que hay que tomar en cuenta para la solucion del problema de estimación. Estas consideraciones se examinan a continuación:

• Estimadores no sesgados. El valor esperado del estimador debe ser igual al parametro. $E({\hat {\theta }})=\theta$.

• Consistencia. A medida que aumenta el tamaño de la muestra, ${\hat {\theta }}$ debe converger a $\theta$ con mayor probabilidad.

• Varianza minima. Muchas veces el estimador se escoge como aquel que posea la menor varianza de todos los estimadores no sesgados. Si no existe un estimador de varianza mínima para todos los tamaños de la muestra, es posible seleccionar entre diferentes estimadores a través del cociente de sus varianzas asintóticas. El estimador que tenga la mínima varianza relativa se le conoce como estimador eficiente.

• Suficiencia. Se dice que un estimador ${\hat {\theta }}$ es suficiente si contiene toda la información en las observaciones para estimar $\theta$. Esto es, el conocimiento directo de los valores $x_{1},\ldots ,x_{N}$ no proporcionará más información que la que contiene ${\hat {\theta }}$. En términos matemáticos, la distribución conjunta de las $x_{i}$ para un valor fijo de ${\hat {\theta }}$ no depende de $\theta$.

Algunas Observaciones Sobre el Sesgo:

El sesgo es la pesadilla de un estadístico. Es fácil de crear y muy díficil de domar. En ocasiones es hasta imposible de controlar. 
La definición estadística de sesgo es la sobreestima o subestima sistemática del verdadero valor. En términos un poco más claros, quiere decir que los resultados siempre estarán 
"movidos" hacia un lado en una cierta cantidad. Por ejemplo, la balanza donde me peso en la oficina del médico, siempre reporta un peso que es dos kilogramos más de lo correcto 
(en serio). Este desplazamiento consistente de dos kilogramos adicionales representa una sobreestima del peso real.
Lo más importante que hay que recordar con respecto al sesgo, es que hay que prevenirlo, o al menos minimizarlo. El sesgo es como la maleza en un jardín. Una vez que se propaga, 
es muy difícil eliminarla. Por eso es mejor eliminarla desde un principio.

Ahora bien, en cada análisis los criterios mencionados arriba deben ser evaluados con respecto al costo de obtener las muestras, la velocidad de procesamiento de los datos, y las consecuencias que puedan resultar de sesgos pequeños o varianzas grandes sobre el resultado de los experimentos. Como se verá más adelante, algunas veces es posible reducir el sesgo a traves de un escalamiento del estimador o una redefinición. En muestras pequeñas la pérdida de eficiencia en algunos estimadores puede ser compensada por su facilidad de cómputo.

Vamos a examinar ahora un método para generar estimaciones conocido como el "método de máxima posibilidad". Este procedimiento es heurístico pero se puede demostrar que los estimadores obtenidos poseen las propiedades deseadas. Empezamos definiendo la función de posibilidad de la muestra aleatoria de observaciones como la densidad conjunta de las variantes de las unidades muestrales evaluadas en $x_{i},\ldots ,x_{N}$:

$L(\theta _{1},\ldots ,\theta _{k})=\prod _{i=1}^{N}f(x_{i};\theta _{1},\ldots ,\theta _{k})$

La función de posibilidad es una medida relativa de la posibilidad de la muestra en particular $x_{i},\ldots ,x_{N}$. El método de máxima posibilidad demanda que los estimadores de $\theta _{i}$ sean seleccionados de tal manera que maximizen la función para una muestra en particular. Si la función de posibilidad tiene un máximo relativo, esto se puede lograr derivando la función y resolviendo las ecuaciones correspondientes. El valor cero de la derivada es sólo una condición necesaria para un máximo relativo y las condiciones suficientes dadas por las derivadas parciales de segundo orden deben verificarse también. Debido a que en muchas de las posibilidades que se encuentran en la inferencia estadística contienen términos exponenciales, se usara el logaritmo natural de la función de posibilidad:

$l(\theta _{i},\ldots ,\theta _{k})=\ln L(\theta _{1},\ldots ,\theta _{k})$

Entonces, si la función de posibilidad tiene un máximo relativo, las estimaciones de $\theta _{i}$ se pueden encontrar resolviendo el sistema de $k$ ecuaciones simultáneas

${\dfrac {\partial l(\theta _{1},\ldots ,\theta _{k})}{\partial \theta _{j}}}=0;\ \ j=1,\ldots ,k$

para encontrar el valor de los estimadores ${\hat {\theta _{1}}},\ldots ,{\hat {\theta _{k}}}$ Si las ecuaciones poseen raíces múltiples, será necesario escoger aquella que nos conduzca la máxima posibilidad.

Veamos dos ejemplos que describen los pasos a tomar para encontrar los estimadores de máxima posibilidad.

Ejemplo

Vamos a determinar los estimadores de máxima posibilidad de $a$ y $b$ para la densidad rectangular. La función de posibilidad de la muestra $x_{1},\ldots ,x_{N}$, es

$L(a,b)={\dfrac {1}{(b-a)^{N}}}$

Es evidente que $a$ no puede ser menor que la observación más pequeña, y que $b$ no puede ser mayor que la observación más grande. Si se denotan las observaciones ordenadas como $x_{(1)}\leq \cdots \leq x_{(N)}$,

$a\leq x_{(1)}\leq \cdots \leq x_{(N)}\leq b$

La función de posibilidad adquiere su valor máximo cuando $(b-a)$ adquiere su valor m+inimo lo que es consistente con el segundo grupo de desigualdades, y las estimaciones que minimizan el intervalo son

${\hat {a}}=x_{(1)}\;\;{\hat {b}}=x_{(N)}$

Ejemplo: La función de posibilidad de la muestra $x_{1},\ldots ,x_{N}$ de $N$ variables aleatorias normales independientes es

$L(\mu ,\sigma ^{2})={\dfrac {1}{(2\pi )^{N/2}(\sigma ^{2})^{N/2}}}\exp \left[-{\dfrac {1}{2\sigma ^{2}}}\sum _{i=1}^{N}(x_{i}-\mu )^{2}\right]$

tomando el logaritmo

$l(\mu ,\sigma ^{2})=-{\dfrac {N}{2}}\ln(2\pi )-{\dfrac {N}{2}}\ln \sigma ^{2}-{\dfrac {1}{2\sigma ^{2}}}\sum _{i=1}^{N}(x_{i}-\mu )^{2}$

Las derivadas parciales con respecto a $\mu$ y $\sigma ^{2}$ son

${\dfrac {\partial l(\mu ,\sigma ^{2})}{\partial \mu }}={\dfrac {1}{\sigma ^{2}}}\sum _{i=1}^{N}(x_{i}-\mu )$

${\dfrac {\partial l(\mu ,\sigma ^{2})}{\partial \sigma ^{2}}}=-{\dfrac {N}{2\sigma ^{2}}}+{\dfrac {1}{2\sigma ^{4}}}\sum _{i=1}^{N}(x_{i}-\mu )^{2}$

Igualando a cero y cancelando factores, queda el siguiente sistema de ecuaciones simultáneas

$\sum _{i=1}^{N}x_{i}-N\mu =0$

$\sum _{i=1}^{N}(x_{i}-\mu )^{2}-N\sigma ^{2}=0$

Despejar el estimador de $\mu$ primero, y luego usar ese valor para obtener la solución del estimador de $\sigma ^{2}$:

${\hat {\mu }}={\dfrac {1}{N}}\sum _{i=1}^{N}x_{i}={\bar {x}}$

${\hat {\sigma }}^{2}={\dfrac {1}{N}}\sum _{i=1}^{N}(x_{i}-{\bar {x}})^{2}$

$={\dfrac {1}{N}}\left[{\dfrac {1}{N}}\sum _{i=1}^{N}x_{i}^{2}-{\dfrac {1}{N}}(\sum x_{i})^{2}\right]$

Estas estimaciones son intuitivamente plausibles, ya que la media de la población se estima a través de la media de la muestra ${\bar {x}}$, y la estimación de la varianza es el promedio de la desviación cuadrada de la media.

Determinemos ahora si ${\hat {\mu }}$ y ${\hat {\sigma }}^{2}$ son estimadores insesgados. Reemplacemos las observaciones por las variables aleatorias $X$ y tomemos las esperanzas:

$E({\hat {\mu }})$

$={\dfrac {1}{N}}\sum _{i=1}^{N}E(X_{i})$

$={\dfrac {1}{N}}\sum _{i=1}^{N}\mu =\mu$

Esto significa que la media muestral es un estimador insesgado de $\mu$. El calcular la esperanza de ${\hat {\sigma }}^{2}$ involucra un poco más de cómputo:

$E({\hat {\sigma }}^{2})$

$={\dfrac {1}{N}}E\left[{\dfrac {1}{N}}\sum _{i=1}^{N}X_{i}^{2}-{\dfrac {1}{N}}\left(\sum X_{i}\right)^{2}\right]$

$={\dfrac {1}{N}}\left[NE(X_{i}^{2})-{\dfrac {1}{N}}E\left(\sum _{i=1}^{N}\sum _{j=1}^{N}X_{i}X_{j}\right)\right]$

$={\dfrac {1}{N}}\left[NE(X_{i}^{2})-E(X_{i}^{2})-{\dfrac {1}{N}}\sum \sum E(X_{i}X_{j})\right]$

$={\dfrac {1}{N}}\left\lbrace (N-1)E(X_{i}^{2})-{\dfrac {1}{N}}N(N-1)[E(X_{i})]^{2}\right\rbrace$

$={\dfrac {N-1}{N}}\sigma ^{2}$

Entonces ${\hat {\sigma }}^{2}$ no es un estimador insesgado de la varianza. Sin embargo, si reemplazamos el divisor $N$ en la fórmula original por $N-1$, se eliminará el sesgo y se obtendrá la conocida formula de la varianza de la muestra

$s^{2}={\dfrac {1}{N-1}}\left[\sum _{i=1}^{N}x_{i}^{2}-{\dfrac {1}{N}}\left(\sum x_{i}\right)^{2}\right].$

Por lo tanto, los estimadores de máxima posibilidad no son necesariamente insesgados, aunque a veces es posible corregir el sesgo multiplicándolos por un factor apropiado. Tales estimadores por lo general no tienen varianza mínima. Sin embargo, se puede demostrar que todos los estimadores de posibilidad máxima son consistentes y que sus distribuciones asintóticas son normales con parametros conocidos.

Pruebas de Hipótesis para los Parámetros de Poblaciones Normales

Pruebas Estadísticas de Hipótesis.

La inferencia estadística puede dividirse en dos áreas generales: 1) La estimación de funciones de distribución, los parámetros de tales funciones cuando se conoce su formulación matemática, o los parámetros de modelos creados alrededor de variables aleatorias; 2) El problema de probar la validez de hipótesis sobre las funciones de distribución y sus parámetros o los componentes de modelos matemáticos. En las sección anterior se estudió brevemente un aspecto de la estimación usando el principio de posibilidad máxima. Ahora se resumirán algunos principios esenciales de las pruebas de hipótesis, y posteriormente veremos como estas pruebas se pueden invertir para proporcionar intervalos de estimación de parámetros.

Quizás un ejemplo sencillo ayudará a comprender el problema de las pruebas de hipótesis en términos de las regiones del espacio de parámetros. Se sabe por experiencias pasadas que el promedio de calificaciones de los estudiantes de primer semestre en cierta universidad canadiense tienden a hallarse normalmente distribuído con media igual a 2.43 y varianza 0.04 (en algunos países las calificaciones van de 0 a 4 con incrementos decimales). Sin embargo, en este ultimo año escolar, los requisitos de entrada se hicieron más estrictos, y se presentó la hipótesis de que la calificiacion media de la población estudiantil entrante será más alta y la varianza más pequeña. Estas afirmaciones sobre los parámetros poblacionales se pueden resumir de la siguiente manera:

Hipótesis original:	$\mu =2.43$	$\sigma ^{2}=0.04$
Hipotesis alternativa:	$\mu >2.43$	$\sigma ^{2}<0.04$

'

Le llamaremos a la descripción original de la población como la 'hipótesis nula; ésta se denota normalmente como $H_{0}$. La hipótesis alternativa será denotada como $H_{1}$. Se supondrá también que la población se encuentra normalmente distribuída, pero por lo general se asume por omisión. $H_{0}$ se refiere al punto (2.43, 0.04) en el espacio de parámetros y se le conoce como una hipótesis sencilla. La hipótesis alternativa se refiere a la región sombreada de la figura 5; y ya que ese conjunto contiene mas de un punto, $H_{1}$ se conoce como una hipótesis compuesta. Un conjunto importante de hipótesis compuestas se forma por aquellos enunciados donde los valores de uno o más parámetros son completamente desconocidos. Por ejemplo, si la variable aleatoria $X$ se encuentra normalmente distribuída con varianza desconocida, las hipótesis

$H_{0}:\;\mu =\mu _{0}$

y

$H_{1}:\;\mu =\mu _{1}>\mu _{0}$

que sólo abarcan la media, se interpretarían como líneas verticales en el espacio de la figura 4 cortando al eje de las $\mu$ en $\mu _{0}$ y en $\mu _{1}$. Habría entonces que desarrollar pruebas estadísticas para tal hipótesis de tal forma que no se afectaran por el valor desconocido de $\sigma ^{2}$.

Fig. 5

A través del desarrollo de la Estadística, el objetivo de la teoría de inferencias estadísticas ha sido el desarrollar pruebas para la validez de estas hipótesis basadas en observaciones muestrales. Para el caso mas general, sea $X$ una variable aleatoria continua con los números reales como su dominio, y su función de densidad depende de los parámetros $\theta _{1},\ldots ,\theta _{k}$. El espacio muestral $W$ de todos los resultados posibles de $N$ observaciones en $X$ ser el espacio N-dimensional euclideano. Una muestra particular de observaciones se expresará en forma de vector como $\left[x_{1},\ldots ,x_{N}\right]$ y denotará un punto en el espacio muestral. Sean $\omega _{1}$ y $\omega _{2}$ dos regiones disjuntas en el espacio de parámetros $\Omega$, y se proponen las siguientes hipótesis sobre los parámetros de la distribución de $X$.

$H_{0}:[\theta _{1},\ldots ,\theta _{k}]\quad {\text{se encuentra contenida en }}\omega _{1}$
$H_{1}:[\theta _{1},\ldots ,\theta _{k}]\quad {\text{se encuentra contenida en }}\omega _{2}$

Basándonos en las observaciones muestrales deseamos decidir en alguna forma "óptima" cuál de las dos hipótesis es la más viable. Debido a que ciertos conjuntos de observaciones nos llevarían a favorecer $H_{0}$ sobre $H_{1}$ y otros conjuntos apoyarían la alternativa, la regla de decisión debe tener la siguiente forma:

${\text{Aceptar }}H_{0}{\text{ si }}\left[x_{1},\ldots ,x_{N}\right]{\text{cae en }}W-w$
${\text{Aceptar }}H_{1}{\text{ si }}\left[x_{1},\ldots ,x_{N}\right]{\text{cae en }}w$

donde $w$ es una parte en particular del espacio muestral llamada región crítica o región de rechazo de $H_{0}$. Si el verdadero estado de la naturaleza se encuentra descrito por $H_{0}$ o $H_{1}$, el tomador de decisiones puede incurrir en dos tipos de error al aplicar la regla de decisión. Un error del primer tipo, o error Tipo I, consiste en declarar $H_{1}$ como verdadera cuando en realidad $H_{0}$ es la correcta. Un error del segundo tipo, o error Tipo II se comete cuando se toma $H_{0}$ como verdadera cuando es $H_{1}$ la que describe el verdadero estado de la naturaleza. Resumimos los resultados de las acciones anteriores en la siguiente tabla:

	Estado de la naturaleza
Decisión	$H_{0}$ verdadera	$H_{1}$ verdadera
Se acepta $H_{0}$	Correcto	Error Tipo II
Se acepta $H_{1}$	Error Tipo I	Correcto

Las probabilidades de los errores de Tipo I y de Tipo II proporcionan medidas de la eficacia de la regla de decisión. La probabilidad $\alpha$ de un error de Tipo I es la probabilidad de que las observaciones muestrales se encuentren en la región crítica $w$ cuando $H_{0}$ es verdadera, o lo mismo que

$\alpha =P(\left[x_{1},\ldots ,x_{N}\right]\in w|H_{0}{\text{ es verdadera}})$

donde el símbolo $\in$ denota la inclusión en el conjunto $w$ y la barra vertical indica que el enunciado de probabilidad esta condicionado a la certeza de la hipótesis nula. El símbolo $\alpha$ se conoce como el tamaño de la prueba. La probabilidad de un error del Tipo II es

$\beta =P(\left[x_{1},\ldots ,x_{N}\right]\in W-w|H_{1}{\text{ es verdadera}})$

El complemento $1-\beta$ de la probabilidad del segundo error se conoce como la potencia de la prueba o la regla de decisión. Si la potencia se calcula para un continuo de valores de los parámetros, las probabilidades resultantes constituyen la función de potencia.

Si el tamaño de la muestra es fijo, los cambios en la forma de la región crítica que reduzcan $\alpha$ también aumentarán $\beta$. Al mismo tiempo, un aumento de $\alpha$ resultará en una reducción de $\beta$. El enfoque clásico de las pruebas de hipótesis requiere una prueba con $\alpha$ fijo cuya región de rechazo se elige de tal manera que minimize $\beta$, o que maximize la potencia. Si ambas hipótesis son sencillas, esto es, tienen la forma

	$\theta _{1}=\theta _{10}$			$\theta _{1}=\theta _{11}$
$H_{0}:$	$\cdots$		$H_{1}:$	$\cdots$
	$\theta _{k}=\theta _{k0}$			$\theta _{k}=\theta _{k1}$

entonces el lema de Neyman-Pearson afirma que la prueba más poderosa con tamaño $\alpha$ tendrá una región crítica definida por la siguiente regla de decisión:

${\text{Aceptar }}H_{0}{\text{ si }}\ \lambda ={\dfrac {f(x_{1};\theta _{10},\ldots ,\theta _{k0})\cdots f(x_{N};\theta _{10},\ldots ,\theta _{k0})}{f(x_{1};\theta _{11},\ldots ,\theta _{k1})\cdots f(x_{N};\theta _{11},\ldots ,\theta _{k1})}}>c$

y

${\text{Aceptar }}H_{1}{\text{ si }}\ \lambda <c$

donde $c$ es una constante tal que

$P(\lambda <c)=\alpha$

El lema define la región crítica como aquel conjunto de puntos en el espacio de muestras para el cual el cociente de posibilidad $\lambda$ es menor que $c$.

Nótese que el lema de Neyman-Pearson requiere que ambas hipótesis sean simples, esto es, que las regiones $\omega _{1}$ y $\omega _{2}$ del espacio de parámetros sean puntos. Las pruebas de hipótesis compuestas que involucren varios parámetros pueden ser construídas usando el poderoso criterio generalizado del cociente de posibilidad:

$\lambda ={\dfrac {L({\hat {\omega }})}{L({\hat {\Omega }})}}$

donde

$L({\hat {\Omega }})=f(x_{1};{\hat {\theta }}_{10},\ldots ,{\hat {\theta }}_{k0})\cdots f(x_{N};{\hat {\theta }}_{10},\ldots ,{\hat {\theta }}_{k0})$

es la función de posibilidad maximizada bajo la suposición de que

$H_{0}:\left[\theta _{1},\ldots ,\theta _{k}\right]\in \omega$

es cierta, y que $L({\hat {\Omega }})$ es la posibilidad maximizada con $\theta _{1},\ldots ,\theta _{k}$ la cual se permite tomar valores a través del todo el espacio de parámetros $\Omega$. Aceptamos $H_{0}$ si $\lambda >c$ y si no, se acepta la alternativa

$H_{1}:\left[\theta _{1},\ldots ,\theta _{k}\right]\in \Omega -\omega$

La constante $c$ se escoge de tal forma que $P(\lambda <c|H_{0}{\text{ es verdadera}})\leq \alpha$. Tal como en el caso de las hipótesis compuestas, cuando el verdadero tamaño de la prueba es realmente menor o igual a $\alpha$, entonces decimos que la prueba es de nivel $\alpha$. Se puede demostrar que si $H_{0}$ es verdadera, la variante

$\chi ^{2}=-2\ln \lambda$

tiende a encontrarse distribuída en forma ji-cuadrada a medida que aumenta el tamaño de la muestra con un número de grados de libertad igual a la diferencia de dimensiones entre el espacio de parámetros $\Omega$ y el subespacio de la hipótesis nula $\omega$, o lo que es mismo, el número de parámetros determinados por $H_{0}$.

Pruebas de la Media de una Variante Normal con Varianza Conocida

Sea $x_{1},\ldots ,x_{N}$ una muestra de observaciones independientes de la variable aleatoria con distribución $N(\mu ,\sigma ^{2})$. La varianza $\sigma ^{2}$ es conocida aunque $\mu$ no. Basándonos sobre las observaciones muestrales deseamos probar la hipótesis

$H_{0}:\mu =\mu _{0}$

de que la media poblacional tiene un cierto valor determinado $\mu _{0}$ contra la hipótesis alternativa:

$H_{1}:\mu =\mu _{1}>\mu _{0}$

que la media es de un valor mayor $\mu _{1}$. Usando el lema de Neyman-Pearson la prueba mas poderosa de tamaño $\alpha$ se basa en la estadística de prueba:

$z={\dfrac {({\bar {x}}-\mu _{0}){\sqrt {N}}}{\sigma }}$

cuya region crítica es

$z>z_{\alpha }$

donde $z_{\alpha }$ es el punto que divide el $100\alpha \%$ superior del resto de la distribucion normal unitaria $N(0,1)$. En términos de la media de las observaciones originales debemos rechazar $H_{0}$ en favor de $H_{1}$ si

${\bar {x}}>\mu _{0}+{\dfrac {\sigma }{\sqrt {N}}}z_{\alpha }$

y si no, entonces aceptar la hipótesis nula. Si la hipotesis alternativa hubiera sido

$H_{1}':\mu =\mu _{1}<\mu _{0}$

entonces la misma estadística de prueba hubiera sido empleada, aunque la región crítica hubiera sido

$z<-z_{\alpha }$

o, equivalentemente

${\bar {x}}<\mu _{0}-{\dfrac {\sigma }{\sqrt {N}}}z_{\alpha }$

Las pruebas anteriores y las hipótesis alternativas $H_{1},H'_{1}$ se conocen como hipótesis unilaterales, debido a que se indica claramente la dirección del cambio de $\mu _{0}$ a $\mu _{1}$. En la medida en que es posible hacer tales predicciones tomando en cuenta el contexto o por previas investigaciones, la potencia de las pruebas será considerablemente mayor que la de las pruebas bilaterales las cuales tienen la hipótesis alternativa:

$H_{1}^{''}:\mu =\mu _{1}\neq \mu _{0}$

lo cual permite mayores o menores valores alternativos de la media. La estadística de prueba sigue siendo , pero la región de rechazo para una prueba de tamaño $\alpha$ es de

$|z|>z_{\alpha /2}$

donde $z_{\alpha /2}$ es el punto que divide el $50\alpha \%$ superior del resto de la distribución normal unitaria. Equivalentemente, se rechaza $H_{0}$ si

${\bar {x}}>\mu _{0}+{\dfrac {\sigma }{\sqrt {N}}}z_{\alpha /2}\;\;\;\;{\text{ó}}\;\;\;\;{\bar {x}}<\mu _{0}-{\dfrac {\sigma }{\sqrt {N}}}z_{\alpha /2}$

son satisfechas.

Las curvas de potencia para las pruebas contra las tres alternativas pueden ser construídas a partir de tablas de la distribución normal, y ya que pueden ser encontradas en muchos textos de estadística, se dejará como ejercicio para la práctica del lenguaje ${\textsf {R}}$. En la práctica uno selecciona un valor apropiado de $\alpha$ y una o más probabilidades $\beta$ que sean tolerables. El tamaño de la muestra $N$ se escoge entonces para lograr que la $\beta$ mínima sea consistente con las limitaciones del presupuesto o el laboratorio.

Finalmente, se hace notar que la función de densidad normal es un miembro de la familia exponencial mencionada anteriormente en esta sección, y por lo tanto las pruebas de hipótesis compuestas tales como

$H_{0}:\mu =\mu _{0}\;\;\;\;H_{1}:\mu >\mu _{0}$

o

$H_{0}:\mu <\mu _{0}\;\;\;\;H_{1}:\mu >\mu _{0}$

son las más poderosas uniformemente.

Pruebas de Medias Cuando se Desconoce la Varianza

En muchas aplicaciones, tanto científicas, tecnológicas, financieras o de otro tipo, no es común que se conozca la varianza de la población por adelantado. Se logró un avance importante en la inferencia estadística cuando W.S. Gosset (en un artículo publicado bajo el seudónimo de "Student") obtuvo la distribución de la estadística de prueba para hipótesis compuestas sobre la media de una distribución normal con varianza desconocida. El criterio generalizado del cociente de posibilidades para probar la hipótesis contra la hipótesis alternativa basándose en $N$ observaciones independientes $x_{1},\ldots ,x_{N}$ con media ${\bar {x}}$ y varianza $s^{2}$ lleva a la estadística de prueba

$t={\dfrac {({\bar {x}}-\mu _{0}){\sqrt {N}}}{s}}$

Si $H_{0}$ es verdadera, $t$ tiene la distribución Student-Fisher con $N-1$ grados de libertad y rechazamos la hipótesis nula en una prueba con tamaño $\alpha$ si

$t>t_{\alpha ;N-1}$

donde $t_{\alpha ;N-1}$ es el punto que separa el $100\alpha \%$ superior de la distribución $t$ definida en la sección anterior. De la misma forma, si la hipótesis alternativa hubiera sido $H_{1}^{'}$, la región de rechazo estaria definida por

$t<-t_{\alpha ;N-1}$

y para la alternativa bilateral, \ref{H1biprime} la hipótesis nula será rechazada si

$|t|>t_{\alpha /2;N-1}$

La potencia de las pruebas que involucran la estadística $t$ se puede calcular de la rutina Pearson-Hartley con $\nu _{1}=1$ y $\nu _{2}=N-1$ grados de libertad y parametro de descentralidad

$\phi ={\dfrac {|\mu -\mu _{0}|}{\sigma }}\left({\dfrac {N}{2}}\right)^{1/2}$

Entonces se puede determinar, dado un $\alpha$ fijo, un tamanño de muestra que garantizará una probabilidad de potencia por encima de un mínimo especificado.

Ahora considérese que dos muestras independientes se han tomado de poblaciones normalmente distribuídas con una varianza común de $\sigma ^{2}$ pero con medias posiblemente distintas $\mu _{1}$ y $\mu _{2}$. Denótense las observaciones de las muestras como $x_{1},\ldots ,x_{N_{1}},y_{1},\ldots ,y_{N_{2}}.$ Entonces el criterio generalizado del cociente de posibilidades para la prueba de hipótesis

$H_{0}:\mu _{1}=\mu _{2}$

de medias poblacionales iguales contra la alternativa

$H_{1}:\mu _{1}>\mu _{2}$ 

lleva a la estadística de prueba

$t={\dfrac {{\bar {x}}-{\bar {y}}}{\sqrt {{\dfrac {\sum _{i=1}^{N_{1}}(x_{i}-{\bar {x}})^{2}+\sum _{i=1}^{N_{2}}(y_{i}-{\bar {y}})^{2}}{N_{1}+N_{2}-2}}\left({\dfrac {1}{N_{1}}}+{\dfrac {1}{N_{2}}}\right)}}}$

Para una prueba con tamaño $\alpha$ se rechaza $H_{0}$ en favor de $H_{1}$ si

$t>t_{\alpha ;N_{1}+N_{2}-2}$

Las regiones de rechazo para las otras hipótesis alternativas $H'_{1}:\mu _{1}<\mu _{2}$ y $H''_{1}:\mu _{1}\neq \mu _{2}$ son similares a \ref{Trej1} y a \ref{Trej2} de las pruebas con una sola muestra.

Intervalos de Confianza para Medias.

La investigadora que ha llevado a cabo un experimento costoso o complejo, rara vez se siente satisfecha al escuchar que sus observaciones han rechazado tal o cual hipótesis. Si los hallazgos muestran que el nuevo medicamento o el nuevo tratamiento tiene algún efecto significativo mas allá del efecto placebo o de las normas anteriores, la investigadora y la comunidad científica no solo preferirían saber la mejor estimación de la magnitud del efecto sino también un intervalo de valores razonables del parámetro de interés. Tales aseveraciones de los posibles valores se conocen como intervalos de confianza o, a diferencia de las estimaciones puntuales de la Sección 1.4, estimaciones de intervalo.

Supóngase que una muestra aleatoria de $N$ observaciones ha sido tomada de alguna población con densidad continua $f(x;\theta )$. El intervalo de confianza del $100(1-\alpha )\%$ para $\theta$ es aquel conjunto de valores

$t_{1}<\theta <t_{2}$

cuyos límites fueron calculados a partir de los puntos porcentuales de la función de distribución de la estimación ${\hat {\theta }}$ de $\theta$ tal que

$P(t_{1}\leq \theta \leq t_{2})=1-\alpha$

donde $\alpha$ generalmente toma el valor 0.05 o alguna probabilidad más pequeña. $1-\alpha$ se conoce como el coeficiente de confianza del intervalo. Es esencial que el enunciado de probabilidad sea interpretado como "la probabilidad de que el intervalo con extremos $t_{1},t_{2}$ contenga a $\theta$ es de $1-\alpha$" ya que en este contexto el parámetro no es una variable aleatoria. También hacemos notar que existe un número infinito de intervalos de confianza que satisfacen \ref{ConfiNt}; en muchos casos se eligirán $t_{1}$ y $t_{2}$ de tal manera que la longitud de $t_{2}-t_{1}$ sea la más corta.

Los intervalos de confianza para la media de una población normal o la diferencia de las medias de dos poblaciones normales se puede deducir de las anteriores pruebas de hipótesis. Si se han tomado $N$ observaciones independientes con media ${\bar {x}}$ de la variante $N(\mu ,\sigma ^{2})$ con varianza conocida, $H_{0}:\mu =\mu _{0}$ se rechaza en favor de

$H_{1}:\mu =\mu _{1}>\mu _{0}$

usando la prueba de tamaño $\alpha _{1}$ cuando el valor de $\mu _{0}$ está dado por

$\mu _{0}={\bar {x}}-{\dfrac {\sigma }{\sqrt {N}}}z_{\alpha _{1}}$

De la misma forma, el valor más grande de $\mu _{0}$ para el cual, usando una prueba con tamaño $\alpha -\alpha _{1}$, $H_{0}$ es rechazada en favor de $H'_{1}$ cuando $H'_{1}:\mu =\mu _{1}<\mu _{0}$, entonces ese valor de $\mu _{0}$ será

$\mu _{0}={\bar {x}}+{\dfrac {\sigma }{\sqrt {N}}}z_{\alpha -\alpha _{1}}$

Donde $z_{\alpha _{1}}$y $z_{\alpha -\alpha _{1}}$ son los porcentajes superiores $100\alpha _{1}$ y $100(\alpha -\alpha _{1})$ de la distribución normal unitaria. Entonces el intervalo de confianza del $100(1-\alpha )$ por ciento para $\mu$ es

${\bar {x}}-{\dfrac {\sigma }{\sqrt {N}}}z_{\alpha _{1}}\leq \mu \leq {\bar {x}}+{\dfrac {\sigma }{\sqrt {N}}}z_{\alpha -\alpha _{1}}$

y su longitud es de $(\sigma /{\sqrt {N}})(z_{\alpha _{1}}+z_{\alpha -\alpha _{1}})$. Pero se puede demostrar que la longitud mínima se obtiene cuando

$z_{\alpha _{1}}=z_{\alpha -\alpha _{1}}$ 

o si

$\alpha _{1}={\frac {1}{2}}\alpha$. 

Entonces el intervalo más corto con coeficiente $1-\alpha$ tiene la forma simétrica

${\bar {x}}-{\dfrac {\sigma }{\sqrt {N}}}z_{\alpha /2}\leq \mu \leq {\bar {x}}+{\dfrac {\sigma }{\sqrt {N}}}z_{\alpha /2}$

De la misma forma, si $\sigma ^{2}$ no es conocida, el intervalo de confianza del $100(1-\alpha )$ por ciento para $\mu$ está dado por

${\bar {x}}-{\dfrac {s}{\sqrt {N}}}t_{\alpha /2;N-1}\leq \mu \leq {\bar {x}}+{\dfrac {s}{\sqrt {N}}}t_{\alpha /2;N-1}$

donde $s$ es la desviación estándar de la muestra, y $t_{\alpha /2;N-1}$ es el $50\alpha \%$ superior de la distribución $t$ con $N-1$ grados de libertad.

Muchas veces es necesario obtener un intervalo de confianza para la diferencia de medias de dos poblaciones normales con una varianza común aunque desconocida. Las observaciones de la primera población quizá fueron obtenidas como control para aquellas en la segunda muestra que fueron obtenidas bajo un nuevo tratamiento o condición experimental. Bajo las suposiciones de normalidad y varianza común, el intervalo de confianza del $100(1-\alpha )$ por ciento para el cambio atribuible a la condición experimental es

${\bar {x}}-{\bar {y}}-s_{{\bar {x}}-{\bar {y}}}t_{\alpha /2;N_{1}+N_{2}-2}\leq \mu _{1}-\mu _{2}\leq {\bar {x}}-{\bar {y}}+s_{{\bar {x}}-{\bar {y}}}t_{\alpha /2;N_{1}+N_{2}-2}$

donde ${\bar {x}}$ y ${\bar {y}}$ son las medias de la primera y segunda muestras compuestas de $N_{1}$ y $N_{2}$ observaciones y

$s_{{\bar {x}}-{\bar {y}}}={\sqrt {{\dfrac {\sum (x_{i}-{\bar {x}})^{2}+\sum (y_{i}-{\bar {y}})^{2}}{N_{1}+N_{2}+2}}\left({\dfrac {1}{N_{1}}}+{\dfrac {1}{N_{2}}}\right)}}$

es la estimación muestral de la desviación estándar de la diferencia de medias.

Pruebas e Intervalos de Confianza para la Varianza

Las pruebas multivariables y declaraciones de confianza de los siguientes capitulos generalmente se van a aplicar a las medias y a otros parámetros de ubicación. Sin embargo, para efectos de integridad, veremos brevemente algunas hipótesis e intervalos de estimación de la varianza de una población normalmente distribuída. Si las observaciones $x_{1},\ldots ,x_{N}$ constituyen una muestra aleatoria de $N(\mu ,\sigma ^{2})$, entonces la cantidad

${\dfrac {(N-1)s^{2}}{\sigma ^{2}}}={\dfrac {\sum _{i=1}^{N}(x_{i}-{\bar {x}})^{2}}{\sigma ^{2}}}$ 

se encuentra distribuida como una variante ji-cuadrada con $N-1$ grados de libertad. El criterio generalizado del cociente de posibilidad para la prueba de la hipótesis

$H_{0}:\sigma ^{2}=\sigma _{0}^{2}$ 

contra la hipótesis alternativa

$H_{1}:\sigma ^{2}>\sigma _{0}^{2}$ 

especifica que la región de rechazo para una prueba de tamaño $\alpha$ es

${\dfrac {(N-1)s^{2}}{\sigma _{0}^{2}}}>\chi _{\alpha ;N-1}^{2}$ 

donde $\chi _{\alpha ;N-1}^{2}$ es el $100\alpha$ porciento superior de la distribución ji-cuadrada con $N-1$ grados de libertad. De la misma forma, las regiones de rechazo para probar la hipótesis nula contra las alternativas $H'_{1}:\sigma ^{2}<\sigma _{0}^{2}$ y $H''_{1}:\sigma ^{2}\neq \sigma _{0}^{2}$ sería

${\dfrac {(N-1)s^{2}}{\sigma _{0}^{2}}}<\chi _{1-\alpha ;N-1}^{2}$ 

y

${\dfrac {(N-1)s^{2}}{\sigma _{0}^{2}}}<\chi _{1-\alpha ;N-1}^{2}\;\;o\;\;{\dfrac {(N-1)s^{2}}{\sigma _{0}^{2}}}>\chi _{\alpha _{2};N-1}^{2}$ 

respectivamente, donde en el segundo caso $\alpha _{1}+\alpha _{2}=\alpha$. Estrictamente hablando, $\alpha _{1}$ y $\alpha _{2}$ deben elegirse de tal forma que la prueba sea insesgada, esto es, su función de potencia no debe ser menor que su tamaño $\alpha$, pero para muchas aplicaciones con muestras grandes es válido partir $\alpha$ en dos partes iguales. Los intervalos de confianza para $\sigma ^{2}$ se pueden obtener directamente de la región de rechazo \ref{RejChi} como se demostrará más adelante.

Si las muestras independientes de $N_{1}$ y $N_{2}$ observaciones se han tomado aleatoriamente de las poblaciones $N(\mu _{1},\sigma _{1}^{2})$ y $N(\mu _{2},\sigma _{2}^{2})$, entonces la hipótesis

$H_{0}:\sigma ^{2}=\sigma _{0}^{2}$ 

se puede probar contra la hipótesis alternativa

$H_{1}:\sigma ^{2}>\sigma _{0}^{2}$ 

por medio de la estadística

$F={\dfrac {s_{1}^{2}}{s_{2}^{2}}}$ 

Si la hipótesis nula es verdadera, la estadística tiene una distribución $F$ con $N_{1}-1,N_{2}-1$ grados de libertad. La región de rechazo para una prueba de tamano $\alpha$ es

$F>F_{\alpha ;N_{1}-1,N_{2}-1}$ 

Por otro lado, $H_{0}$ se puede probar contra la otra hipótesis alternativa $H'_{1}:\sigma _{1}^{2}<\sigma _{2}^{2}$ por medio de la estadística

$F={\dfrac {s_{2}^{2}}{s_{1}^{2}}}$ 

cuya region crítica es

$F>F_{\alpha _{1};N_{1}-1,N_{2}-1}\ \ o'\ \ F<{\dfrac {1}{F_{\alpha _{2};N_{1}-1,N_{2}-1}}}$ 

donde la repartición $\alpha _{1}=\alpha _{2}={\frac {1}{2}}\alpha$ será suficiente.

Pruebas de Igualdad de Medias: Análisis de Varianza

Presentaremos ahora el ejemplo más sencillo del modelo lineal general que es la base del diseño experimental estadístico. La teoría de estimación y pruebas de hipótesis en modelos lineales de una sola variable ha sido discutida ampliamente en la literatura y mientras que no se presentará aquí en su totalidad, algunas partes del tema son de particular importancia para el trasfondo teórico de esta sección y servirán como una introducción al tratamiento de modelos multivariables en el capítulo correspondiente

Supóngase que se sabe que cierto compuesto bioquímico es absorbido por el cerebro aunque existe cierta evidencia que la cantidad de la sustancia por gramo de tejido cerebral parece diferir entre las cinco cepas de ratones usados en un laboratorio. Se supondrá que las cantidades relativas extraídas de los ratones sacrificados se encuentran normalmente distribuídas con la misma varianza $\sigma ^{2}$ en cada cepa. Sea $\mu _{j}$ la media poblacional para la j-ésima cepa. Es posible entonces construir una prueba de

$H_{0}:\mu _{1}=\cdots =\mu _{5}$ 

contra la hipótesis alternativa de que algunas medias son distintas usando el criterio generalizado del cociente de posibilidad. Aun mas, si $H_{0}$ es rechazada, se pueden obtener métodos para llevar a cabo pruebas simultáneas sobre las diferencias entre las medias con una probabilidad fija de error Tipo I para todas las comparaciones.

En el caso general de $k$ cepas, tratamientos, categorías de diagnóstico, o condiciones experimentales empezamos postulando que la i-ésima observación del j-ésimo tratamiento (tratamiento es el nombre que en general se le da a la característica que distingue a los $k$ grupos) se puede expresar por medio del modelo matemático

$x_{ij}=\mu +\tau _{j}+e_{ij}\ \ i=1,\ldots ,N_{j},j=1,\ldots ,k.$ 

donde

${\begin{aligned}\mu &={\text{parámetro de ubicación común a todas las observaciones}}\\\tau _{j}&={\text{efecto del j-ésimo tratamiento}}\\e_{ij}&={\text{variable aleatoria normalmente distribuída con media igual a cero y varianza}}\;\;\sigma ^{2}\end{aligned}}$

Los términos variantes $e_{ij}$ se encuentran distribuídos independientemente unos de otros. El modelo se dice que es lineal o aditivo, ya que la aplicación del j-ésimo tratamiento aumenta $\mu$ en la cantidad $\tau _{j}$, y la perturbación aleatoria también agrega o sustrae algo de los parámetros. Hacemos notar que, en una versión sencilla del ejemplo inicial, $\mu _{j}\equiv \mu +\tau _{j}$. La hipótesis de los efectos equivalentes del tratamiento se puede expresar como

$H_{0}:\tau _{1}=\cdots =\tau _{k}$ 

y se tomará como la alternativa a $H_{0}$ el modelo general \ref{GenModel} para las observaciones.

Este modelo se dice que tiene efectos fijos, ya que los $\mu _{j}$ son parámetros y cualquier inferencia a partir de las observaciones sólo se puede hacer con respecto a los $k$ tratamientos en el estudio. Si los tratamientos constituyen una muestra aleatoria de una población mayor (técnicos de laboratorio, experimentos realizados en distintos dias, etc.) entonces los efectos del tratamiento serían una variable aleatoria en sí mismos y se tendría que tomar un enfoque distinto para llevar a cabo el análisis. Esta distinción entre los modelos fijos y aleatorios de análisis de varianza fue hecha por primera vez por Eisenhart (1947) y ha sido desarrollada extensamente por otros autores en el area de diseño de experimentos. Más adelante se usaran técnicas multivariables para el análisis exacto del modelo mixto para un número fijo de tratamientos aplicados repetidamente a cada miembro de una muestra aleatoria de unidades experimentales.

Las observaciones se pueden ordenar como se muestra en la siguiente tabla

$x_{11}$	$\cdots$	$x_{1k}$
$\cdots$	$\cdots$	$\cdots$
$x_{N_{1}1}$	$\cdots$	$x_{N_{k}k}$

Denótese el total de las observaciones para el j-ësimo tratamiento como

$T_{j}=\sum _{i=1}^{N_{j}}x_{ij}$ 

y se denotará la media de ese tratamiento como ${\bar {x}}_{j}=T_{j}/N_{j}$. La suma de todas las observaciones será denotada por

$G=T_{1}+\cdots +T_{k}$ 

el número total de unidades experimentales por

$N=N_{1}+\cdots +N_{k}$ 

y la media general por ${\bar {x}}=G/N$. Expresamos entonces la suma total de los cuadrados

${\begin{aligned}S&=\sum _{j=1}^{k}\sum _{i=1}^{N_{j}}(x_{ij}-{\bar {x}})^{2}\\&=\sum _{j=1}^{k}\sum _{i=1}^{N_{j}}x_{ij}^{2}-{\dfrac {G^{2}}{N}}\end{aligned}}$

como la suma dos componentes independientes

$S=SST+SSE$

donde

${\begin{aligned}SST&=\sum _{j=1}^{k}N_{j}\left({\dfrac {T_{j}}{N_{j}}}-{\bar {x}}\right)^{2}\\&=\sum _{j=1}^{k}{\dfrac {T_{j}^{2}}{N_{j}}}-{\dfrac {G^{2}}{N}}\end{aligned}}$

y

${\begin{aligned}SSE&=\sum _{j=1}^{k}\sum _{i=1}^{N_{j}}(x_{ij}-{\bar {x}}_{j})^{2}\\&=\sum _{j=1}^{k}\sum _{i=1}^{N_{j}}x_{ij}^{2}-\sum _{j=1}^{k}{\dfrac {T_{j}^{2}}{N_{j}}}\end{aligned}}$

Fuente	Suma de cuadrados	Grados de Libertad	Cuadrado Medio
Tratamiento	SST	$k-1$	${\dfrac {SST}{k-1}}$
Intra-tratamiento	SSE	$N-k$	${\dfrac {SSE}{N-k}}$
Total	S	$N-1$

Estos componentes se pueden resumir en la tabla de análisis de varianza mostrada en la tabla de arriba. La estadística para la prueba del cociente de posibilidad de $H_{0}$ es

$F={\dfrac {N-k}{k-1}}{\dfrac {SSE}{SST}}$ 

y, cuando $H_{0}$ es verdadera, $SST/(k-1)\sigma ^{2}$, $SSE/(N-k)\sigma ^{2}$ son variantes independientes con distribucion ji-cuadrada con $k-1$ y $N-k$ grados de libertad, respectivamente. Entonces la estadística de prueba tiene una distribución $F$ con $k-1$ y $N-k$ grados de libertad. Rechazamos $H_{0}$ con una prueba de nivel $\alpha$ si

$F>F_{\alpha ;k-1,N-k}$ 

La función de potencia para esta prueba de análisis de varianza se puede calcular a partir de la distribución F no central Pearson-Hartley. Los grados de libertad son $\nu _{1}=k-1$ y $\nu _{2}=N-k$, y el parámetro de no centralidad que mide la desviación de la hipótesis nula \ref{H0detau} de las medias poblacionales de los tratamientos es

$\phi ={\dfrac {\sqrt {\sum _{j=1}^{k}N_{j}\left(\tau _{j}-{\dfrac {\sum N_{j}\tau _{j}}{N}}\right)^{2}}}{\sigma {\sqrt {k}}}}$ 

o, como en el caso de muestras con tratamientos equivalentes que lo que generalmente se encuentra en problemas de diseño experimental

$\phi ={\dfrac {\sqrt {N_{j}\sum _{j=1}^{k}\tau _{j}^{2}}}{\sigma {\sqrt {k}}}}$ 

En esta última expresión, el término de corrección ha desaparecido de la habitual restricción $\tau _{1}+\ldots +\tau _{k}=0$ que se impone sobre los efectos del tratamiento. El parámetro $\phi$ y las gráficas de Pearson-Hartley se usan normalmente para seleccionar tamaños de muestras.

Comparaciones de Tratamientos Múltiples.

Un análisis de varianza que ha resultado en el rechazo de la hipótesis de efectos equivalentes del tratamiento no indica cuáles de esos efectos pueden ser iguales o cuales son diferentes. Este es el problema de comparaciones múltiples, o inferencias simultáneas sobre los miembros de alguna familia de hipótesis. La pruebas se construyen de tal forma que la probabilidad de error Tipo I para toda la familia será a lo mucho $\alpha$. Dos métodos serán necesarios en los siguientes capítulos, y los presentaremos a continuación en el contexto del análisis de varianza de una vía.

La primera técnica se debe a Scheffé. Defínase un contraste de parámetros llamado $\tau _{j}$ del modelo de una vía como cualquier función lineal

$\sum _{j=1}^{k}c_{j}\tau _{j}$ 

cuyos coeficientes tienen la propiedad de

$\sum _{j=1}^{k}c_{j}=0$ 

Por lo tanto $\tau _{1}-\tau _{2}$ y $3\tau _{1}+(\tau _{2}+\tau _{3}+\tau _{4})$ son contrastes mientras que $\tau _{2}-(\tau _{3}+\tau _{4})$ no lo son. En el caso particular del modelo de análisis de varianza de una vía, Scheffé ha demostrado que los intervalos simultáneos de confianza con coeficiente conjunto $1-\alpha$ para todos los contrastes $\tau _{j}$ tienen la forma

$\sum c_{j}{\bar {x}}_{j}-s{\sqrt {(k-1)F_{\alpha ;k-1,N-k}\sum {\dfrac {c_{j}^{2}}{N_{j}}}}}\leq \sum c_{j}\tau _{j}\leq \sum c_{j}{\bar {x}}_{j}+s{\sqrt {(k-1)F_{\alpha ;k-1,N-k}\sum {\dfrac {c_{j}^{2}}{N_{j}}}}}$ 

donde todas las sumas son sobre los $k$ tratamientos. $\sum c_{j}{\bar {x}}_{j}$ es la estimación muestral del contraste $\sum c_{j}\tau _{j}$ y

$s={\sqrt {\dfrac {SSE}{N-k}}}$ 

Nótese que $s^{2}+\sum c_{j}^{2}N_{j}$ es la estimación de la varianza del contraste estimado. Aceptamos la hipótesis nula

$H_{0}:\sum _{j=1}^{k}c_{j}\tau _{j}=0$ 

al nivel $\alpha$ si el intervalo simultáneo de confianza para ese contraste incluye el valor cero. Si, por otro lado

$\sum c_{j}{\bar {x}}_{j}>s{\sqrt {(k-1)F_{\alpha ;k-1,N-k}\sum {\dfrac {c_{j}^{2}}{N_{j}}}}}$ 

o

$\sum c_{j}{\bar {x}}_{j}<-s{\sqrt {(k-1)F_{\alpha ;k-1,N-k}\sum {\dfrac {c_{j}^{2}}{N_{j}}}}}$ 

rechazamos $H_{0}$ en favor de la hipótesis alternativa unilateral $H_{1}:\sum c_{j}\tau _{j}>0$ o $H'_{1}:\sum c_{j}\tau _{j}<0$. El nivel conjunto de todas las pruebas es de $\alpha$.

La segunda técnica de comparación múltiple se conoce como el método de Bonferroni ya que se encuentra basada en una desigualdad con ese nombre. Para este procedimiento empezamos enfocándonos en una familia de $m$ declaraciones de intervalos de confianza $H_{1},\ldots ,H_{m}$ de los parámetros del modelo lineal. La probabilidad de que $H_{i}$ sea una declaración verdadera (esto es, que $H_{i}$ contenga el valor de la i-ésima función paramétrica) es igual a $P(H_{i})$ y la probabilidad de que todas las declaraciones sean correctas simultáneamente es de $P(H_{1}\cap \cdots \cap H_{m})$, donde la notación de intersección $H_{i}\cap H_{j}$ denota el evento donde ambos $H_{i}$ y $H_{j}$ son verdaderos. Esta probabilidad es el coeficiente simultáneo de confianza para la familia completa de intervalos, y nos gustaría que al menos fuera igual a algún valor especificado $1-\alpha$, donde $\alpha$ será conocida como la proporción de error para la familia de declaraciones. El cálculo de la probabilidad conjunta muchas veces es difícil para problemas prácticos y su valor depende de algunos "parámetros de molestia" que miden las intercorrelaciones de las declaraciones $H_{i}$. En lugar de calcular esta probabiliad, debemos contentarnos con la siguiente cota inferior:

$P(H_{1}\cap \cdots H_{m})\geq 1-\sum _{i=1}^{m}P({\bar {H}}_{i}).$ 

donde $P({\bar {H}}_{i})=1-P(H_{i})$ es la probabilidad de que la i-ésima declaración individual no sea verdadera. La cota es un ejemplo sencillo de las desigualdades de Bonferroni. Para un conjunto de $m$ intervalos de confianza simultáneos, por lo general se le asigna a cada declaración una proporción de error de $\alpha$, de tal manera que el coeficiente para toda la familia es de al menos $1-\alpha$.

Los intervalos de confianza de Bonferroni sobre $m$ contrastes dados

$\psi _{i}=\sum _{j=1}^{k}c_{ij}\tau _{j}\;\;\;i=1,\ldots ,m$ 

de los parámetros del modelo lineal de una vía son

$\sum _{j=1}^{k}c_{ij}{\bar {x}}_{j}-t_{\alpha /(2m);N-k}s{\sqrt {\sum _{j=1}^{k}{\dfrac {c_{ij}^{2}}{N_{j}}}}}\leq \Psi _{i}\leq \sum _{j=1}^{k}c_{ij}{\bar {x}}_{j}+t_{\alpha /(2m);N-k}s{\sqrt {\sum _{j=1}^{k}{\dfrac {c_{ij}^{2}}{N_{j}}}}}$ 

La regla de decisión para probar $H_{0}:\Psi _{i}=0$ es equivalente a aquella para el metodo de Scheffe: si (\ref{Bonf01}) no incluye el valor cero, se rechaza $H_{0}$. Para todas las comparaciones apareadas de los tratamientos, $m={\tfrac {1}{2}}k(k-1)$, mientras que para comparaciones sucesivas (suponiendo que hubo algun ordenamiento a priori), $m=k-1$. Si $m$ es pequeña los intervalos de Bonferroni pueden ser en promedio mas estrechos que los de la técnica de Scheffé, aun cuando el verdadero coeficiente de confianza de la familia es mayor que su valor nominal $1-\alpha$. Para valores de $m$ muy grandes, los intervalos de Bonferroni son imprácticamente grandes.

Ejemplo 3: En una evaluación preliminar de tres medicamentos antidepresivos, las condiciones del experimento decretaron que cada sujeto participante podía recibir solo un medicamento. Dieciocho pacientes con diagnósticos similares fueron calificados en su nivel de angustia usando una escala de cero a siete y se tomaron tres grupos de seis aleatoriamente. A cada grupo se le administró uno de los tres medicamentos, y despues de tres días los pacientes fueron calificados de nuevo usando la misma escala. Estos fueron los cambios en los niveles de angustia:

	A	B	C
	4	0	1
	5	2	0
	3	1	0
	3	2	1
	4	2	2
	2	2	2
Media	3.5	1.5	1.0

Las sumas pertinentes son

$T_{1}=21\;\;\;T_{2}=9\;\;\;T_{3}=6\;\;\;G=36\;\;\;\sum _{j=1}^{3}\sum _{i=1}^{6}x_{ij}^{2}=106$

La suma de cuadrados debida a los medicamentos es

$SST={\frac {1}{6}}(21^{2}+9^{2}+6^{2})-{\dfrac {36^{2}}{18}}=21$

y la suma total de cuadrados es $S=106-72=34$. Los valores se resumen en la tabla de análisis de varianza:

Fuente	Suma de cuadrados	Grados de libertad	Media cuadrada	F
Medicamentos	21	2	10.5	12.1
Entre Medicamentos	12	15	0.867
Total	34	17

Ya que $F_{0.01;2,15}=6.36$ concluímos que la hip[otesis de efectos equivalentes del medicamento no es válida al nivel del uno por ciento. Consultando las tablas de la distribución $F$ observamos que la $F$ excede el valor crítico para $\alpha =0.001$

Ahora usaremos el procedimiento de comparaciones múltiples de Scheffé para determinar cuáles medicamentos son los distintos. Se elige un coeficiente simultáneo de confianza de 0.99 y asi:

${\sqrt {(k-1)F_{0.01;2.15}}}=3.567$

Primero compárense los medicamentos B y C calculando el intervalo de confianza para $\tau _{2}-\tau _{3}$. Aqui $c_{1}=0,c_{2}=1,c_{3}=-1$ y la desviación estándar estimada de la población ${\bar {x_{2}}}-{\bar {x_{3}}}$ es de

${\sqrt {0.867({\frac {1}{6}}+{\frac {1}{6}})}}=0.5376$

El intervalo de confianza es $-1.42\leq \tau _{1}-\tau _{2}\leq 2.42$ y ya que esto cubre el valor cero con facilidad, podemos concluír que la hipótesis de efectos equivalentes de los medicamentos B y C no puede ser rechazada al nivel 0.01 de confianza.

De la misma forma, el intervalo simultáneo de confianza para el efecto de la diferencia entre A y B es

$0.08\leq \tau _{1}-\tau _{2}\leq 3.92$

y es posible aceptar la hipótesis alternativa

$H_{1}:\tau _{1}>\tau _{2}$

al nivel 0.01. También será interesante determinar si el medicamento A es distinto del efecto promedio de los medicamentos B y C. Aqui $c_{1}=1,\ c_{2}=-{\frac {1}{2}},\ c_{3}=-{\frac {1}{2}}$ y la desviación estándar estimada de ese contraste es

${\sqrt {0.867({\frac {1}{6}}+{\frac {1}{6}}+{\frac {1}{24}})}}=0.4653$

La estimación del contraste es de 2.25, y el intervalo de confianza es de

$0.59\leq \tau _{1}-{\frac {1}{2}}(\tau _{2}+\tau _{3})\leq 3.91$

El medicamento A parece ser distinto de los demas medicamentos restantes que acabaron siendo esencialmente los mismos.

Los intervalos de confianza de Bonferroni para los contrastes apareados con coeficiente de confianza del 0.99 son

0.13	$\leq \tau _{1}-\tau _{2}\leq$	3.87
-1.37	$\leq \tau _{2}-\tau _{3}\leq$	2.37
0.63	$\leq \tau _{1}-\tau _{3}\leq$	4.37

Introducción al lenguaje R

¿Qué es R?

El sistema ${\textsf {R}}$ para cómputos estadísticos es un ambiente para análisis de datos y gráficos. El origen de ${\textsf {R}}$ es el lenguaje ${\textsf {S}}$ el cual fue desarrollado por John Chambers y sus colegas en los laboratorios Bell por la década de los 1960s. El lenguaje ${\textsf {S}}$ fue diseñado y desarrollado como un lenguaje de programación para análisis de datos pero es, en su implementación presente, un lenguaje completo de programación con todas las características.

La distribución básica de ${\textsf {R}}$ es mantenida por un pequeño grupo de estadísticos, quienes forman el grupo central de desarrollo. Asimismo, un inmenso grupo de voluntarios le agregan al paquete una gran cantidad de funciones que aumentan enormemente su utilidad. La fuente principal de información sobre el sistema ${\textsf {R}}$ se encuentra localizada en el siguiente URL:

http://www.R-project.org

Todos los recursos de ${\textsf {R}}$ se encuentran en esta página: el sistema mismo, los paquetes accesorios, los manuales y la documentación.

La intención de este capítulo es presentar informalmente los conceptos básicos y las técnicas de manipulación de datos del lenguaje ${\textsf {R}}$ para los principiantes. No entraremos en muchos detalles técnicos. Más bien, nos concentraremos en ejemplos prácticos para que el lector pueda empezar a resolver problemas de inmediato.

Instalación del lenguaje

El sistema ${\textsf {R}}$ consiste de dos partes principales: el sistema base y la colección de paquetes accesorios que han sido contribuídos por los usuarios. El lenguaje se encuentra implementado en el sistema base. Las implementaciones de los procedimientos estadísticos y gráficos están separados del sistema base y se encuentran organizados en forma de "paquetes". Un paquete es una colección de funciones, ejemplos y documentación. La función de un paquete muchas veces se concentra en un método estadístico en particular. Ambos el sistema base y la paquetería se distribuyen a traves del "Comprehensive ${\textsf {R}}$ Archive Network" (CRAN), el cual se encuentra en

http://CRAN.R-project.org

El Sistema Base y Primeros Pasos

El sistema base se encuentra disponible en forma de código fuente y en forma binaria para varios sistemas Unix, Linux, Windows y Mac OS X. Para el analista de datos, será suficiente bajar la distribución binaria e instalarla localmente. Los usuarios de Windows pueden seguir el enlace:

http://CRAN.R-project.org/bin/windows/base/release.htm

bajar el archivo ejecutable, correrlo localmente en su máquina, y seguir las instrucciones proporcionadas por el instalador automático. Dependiendo del sistema operativo, ${\textsf {R}}$ puede ejecutarse ya sea escribiendo "R" en una terminal (Unix, Linux), o haciendo doble click en el logotipo del sistema (ver figura). %TODO: Poner logo del sistema R Para verificar que el sistema fue instalado correctamente, teclear el siguiente mando en una terminal de mandos (no incluír el signo de pesos ya que representa el apronte del sistema operativo):

 $ R --version

Aparecerá entonces en la pantalla un corto mensaje y el apronte del sistema R '>' que será algo parecido a lo siguiente:

R version 2.14.1 (2011-12-22)
Copyright (C) 2011 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: x86_64-redhat-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

>

Lo más probable es que la versión varíe un poco como también la plataforma, dependiente ésta de donde se ejecute el programa. Este mensaje nos indica que estamos corriendo la versión 2.14 del sistema R. Si el mensaje no aparece en la pantalla, hubo algún error y no se realizó la instalación correctamente. Consulte a su administrador de sistemas.

Se puede cambiar la apariencia del apronte del sistema ${\textsf {R}}$ con el siguiente mando:

 > options(prompt = "R> ")

Usaremos de ahora en adelante el apronte R> para mostrar los ejemplos de código a través de este libro. Un símbolo + al principio de la línea indica la continuación de un mando después de un regreso de carro.

Básicamente, el sistema ${\textsf {R}}$ evalúa los mandos proporcionados en el apronte y regresa los resultados de los cálculos. El fin de un mando se indica con la tecla de regreso (Return). Casi todos los textos de introducción a ${\textsf {R}}$ empiezan con un ejemplo usando al sistema como una calculadora:

 R> x <- sqrt(25) + 2

Este mando simple le pide al interpretador que evalúe la raíz cuadrada de 25 y le agregue 2. El resultado se lo asigna a un objeto de ${\textsf {R}}$ que tiene el nombre de variable x. El operador de asignación <- afianza el valor a su lado derecho con el nombre de la variable a su lado izquierdo. El valor del objeto x puede ser inspeccionado escribiendo

 R> x
   [1] 7

esto llama implícitamente al método print

 R> print(x)
   [1] 7

Paquetería

La distribución base incluye un conjunto de paquetes de alta prioridad como gráficos, estadísticas básicas, y algunos paquetes estándar tales como modelos lineales, regresiones, y análisis de varianza. Los paquetes que no se encuentren incluídos en la distribución base se pueden instalar directamente desde el apronte de ${\textsf {R}}$. Suponiendo que el usuario se encuentra conectado al Internet, un paquete se puede instalar proporcionando el nombre del paquete a la función install.packages. Por ejemplo, si se requiriera la estimación robusta de matrices de covarianza a través de estimadores tipo "sandwich", el paquete sandwich se puede bajar e instalar a través de:

 R> install.packages("sandwich")

Las funciones del paquete se encontrarán disponibles después de afianzar el paquete por medio de

 R> library("sandwich")

Se puede obtener una lista completa de los paquetes disponibles en

http://CRAN.R-project.org/web/packages/

Nótese que los paquetes bajados a sistemas en Windows estarán precompilados antes de bajarse e instalarse. Esto, a diferencia de los paquetes en Unix/Linux donde los paquetes se compilan localmente antes de instalarse.

Ayuda y Documentación

Existen tres formas de documentación para el sistema ${\textsf {R}}$: la ayuda en línea que viene con la distribución base, manuales electrónicos, y documentación en forma de libros y publicaciones.

La ayuda en línea consiste en una colección de páginas tipo manual de Unix (man) las cuales describen las funciones y los datos que se incluyen con ${\textsf {R}}$. Una página de manual se muestra cuando llamamos a la función help:

 R> help("mean")

o, de forma más corta:

 R> ?mean

Cada página electrónica del manual consiste de una descripción general, la lista de argumentos que requiere la función junto con una descripción de cada argumento, información sobre el valor que regresa la función, y, a veces, referencias, enlaces, y ejemplos ejecutables. La función help.search es útil para buscar dentro de las páginas manuales. Por ejemplo, se puede obtener una panorámica de los temas documentados en un paquete por medio de

 R> help(package="sandwich")

Este mando nos ayuda a obtener un resumen de las funciones del paquete sandwich.

Muchas veces un paquete viene con documentación adicional que describe las funciones del paquete y muestra ejemplos del uso del mismo. Tal documento se conoce como viñeta. Por ejemplo, la viñeta del paquete sandwich se puede abrir usando

 R> vignette("sandwich", package="sandwich")

Se puede encontrar una extensa colección de documentación en

http://CRAN.r-project.org/manuals.html

Otras funciones útiles son find y apropos. La función find nos dice en que paquete se encuentra alguna función:

> find("lowess")
[1] "package:stats"

mientras que apropos regresa con una matriz conteniendo los nombres de todos los objetos que se encuentran en la búsqueda:

  > apropos("lm")
 [1] ".__C__anova.glm"      ".__C__anova.glm.null" ".__C__glm"
 [4] ".__C__glm.null"       ".__C__lm"             ".__C__mlm"
 [7] ".__C__optionalMethod" "anova.glm"            "anova.glmlist"
[10] "anova.lm"             "anova.lmlist"         "anova.mlm"
[13] "anovalist.lm"         "colMeans"             "contr.helmert"
[16] "getAllMethods"        "glm"                  "glm.control"
[19] "glm.convert"          "glm.fit"              "glm.fit.null"
[22] "glm.nb"               "glmmPQL"              "hatvalues.lm"
[25] "KalmanForecast"       "KalmanLike"           "KalmanRun"
[28] "KalmanSmooth"         "kappa.lm"             "lm"

Ejemplos Preparados y Demostraciones

Si se desea ver el ejemplo preparado, solo hay que usar la función example usando como argumento el nombre de la función que se desea. En este ejemplo, usaremos la función de modelos lineales, o lm en esta instancia:

 R> example(lm)

El sistema presentará los resultados y gráficos incluídos en el ejemplo donde se ilustra la función lm.

Para acceder a las demostraciones de algunos paquetes el usuario puede ejecutar los siguientes mandos:

 R> demo(persp)
 R> demo(Hershey)
 R> demo(graphics)
 R> demo(plotmath)

Objetos en R

El paquete básico de ${\textsf {R}}$ incluye algunos conjuntos de datos que se pueden usar para efectos de prueba, calibración o meramente didácticos. Si se desea echar un vistazo a todos los grupos de datos que ofrece ${\textsf {R}}$ basta teclear

 R> data()

y el sistema desplegará una lista de todos los diferentes conjuntos de datos incluídos.

Para ayudar a explicar las técnicas de manipulación de datos, usaremos un grupo llamado HairEyeColor que consiste en una distribución del color de cabello, color de ojos y género en 592 estudiantes provenientes de los diferentes grupos de un curso de estadística en una universidad. Para mayor informacion sobre la naturaleza de los datos, se puede teclear

 R> ?HairEyeColor

y se desplegará la información requerida. En forma resumida, este paquete de datos consiste en los resultados de una muestra de 592 estudiantes de quienes se tomaron los datos antes mencionados. Estos datos serán útiles para ilustrar algunas técnicas para el análisis de tablas de contingencia, tales como la prueba ji-cuadrada, algunos modelos lineales y rutinas de graficado usando ${\textsf {R}}$.

Empezamos con dos mandos básicos: imprimir los datos, y analizar su estructura. Para ahorrar espacio, solamente mostraremos aquí los resultados parciales de ejecutar los mandos. Pediremos al lector que siga la lección en su propia computadora para facilitar el aprendizaje. Se teclea primero

 R> print(HairEyeColor)

y se despliega el siguiente resultado

  , , Sex = Male

       Eye
Hair    Brown Blue Hazel Green
  Black    32   11    10     3
  Brown    53   50    25    15
  Red      10   10     7     7
  Blond     3   30     5     8

, , Sex = Female

       Eye
Hair    Brown Blue Hazel Green
  Black    36    9     5     2
  Brown    66   34    29    14
  Red      16    7     7     7
  Blond     4   64     5     8

El sistema nos está mostrando como están distribuídos los colores de ojos y cabello según el género del estudiante. A continuación se teclea

 R> str(HairEyeColor)

y el sistema despliega

   table [1:4, 1:4, 1:2] 32 53 10 3 11 50 10 30 10 25 ...
 - attr(*, "dimnames")=List of 3
  ..$ Hair: chr [1:4] "Black" "Brown" "Red" "Blond"
  ..$ Eye : chr [1:4] "Brown" "Blue" "Hazel" "Green"
  ..$ Sex : chr [1:2] "Male" "Female"

Este mando presenta la estructura interna de un objeto R. En este caso nos está informando que la estructura del conjunto de datos "HairEyeColor" es una tabla de tres dimensiones: $4\times 4\times 2$. Si se desea visualizar este concepto sería como imaginar dos tablas de cuatro filas con cuatro columnas cada una. Cada una de estas dimensiones tiene un atributo al cual se le asigna un nombre. Como podemos ver en el resultado los nombre de los atributos son "Hair", "Eye", y "Sex". Los atributos son los nombres de tres distintos arreglos de caracteres. Por ejemplo, el atributo "Hair" es el nombre del arreglo de cuatro elementos que contiene los nombres de los distintos colores de cabello. Al imprimir la estructura del objeto "HairEyeColor" el sistema nos indica como están ordenados los atributos de la tabla.

Marcos de Datos (Data Frames)

Veamos ahora la otra estructura interna que ${\textsf {R}}$ usa para almacenar datos. Escribamos el siguiente mando:

 R> str(USArrests)

El sistema responde con:

  'data.frame':   50 obs. of  4 variables:
 $ Murder  : num  13.2 10 8.1 8.8 9 7.9 3.3 5.9 15.4 17.4 ...
 $ Assault : int  236 263 294 190 276 204 110 238 335 211 ...
 $ UrbanPop: int  58 48 80 50 91 78 77 72 80 60 ...
 $ Rape    : num  21.2 44.5 31 19.5 40.6 38.7 11.1 15.8 31.9 25.8 ...

${\textsf {R}}$ nos está mostrando la estructura del conjunto de datos "USArrests" el cual es el número de arrestos por cada 100,000 habitantes por los crímenes de robo, homicidio y violación correspondientes a cada uno de los 50 estados de la Unión Americana en 1973. Se proporciona también el porcentaje de la población que vive en áreas urbanas. Estos datos se encuentran estructurados en la forma de un data frame o "marco de datos". Los objetos tipo dataframe son los objetos más importantes dentro de ${\textsf {R}}$. Estos objetos representan a los datos en la forma de una tabla común. Cada fila se encuentra asociada con una sola observación y cada columna representa una variable. Las dimensiones de esta tabla se pueden extraer usando la función dim:

  R> dim(USArrests)
  [1] 50  4

indicando así que el dataframe asociado al conjunto de datos "USArrests" es una tabla con 50 filas y 4 variables. Se pueden obtener también las dimensiones de la tabla usando los mandos nrow y ncolumn:

  R> nrow(USArrests)
  [1] 50
  R> ncol(USArrests)
  [1] 4

Los nombres de las variables que describen las observaciones, esto es, los nombres de las columnas del dataframe se extraen usando

  R> names(USArrests)
[1] "Murder"   "Assault"  "UrbanPop" "Rape"

Los corchetes siempre indican un subconjunto de un objeto mayor, en este caso una sola variable extraída de la tabla completa. Debido a que los dataframes tienen dos dimensiones, siendo éstas las observaciones y las variables, se requiere la coma para especificar que se desea un subconjunto de la segunda dimensión, esto es, de las variables. Por ejemplo las observaciones correspondientes a los homicidios se encuentran representadas en un vector de longitud

  R> length(USArrests[,"Murder"])
  [1] 50

Un vector es la estructura más elemental para manejar datos en ${\textsf {R}}$ y consiste en un conjunto de elementos sencillos, los cuales son todos de la misma clase. Por ejemplo, un vector simple de los números del uno al cinco se puede construír con cualquiera de los tres siguientes mandos:

R> 1:3
[1] 1 2 3

R> c(1,2,3)
[1] 1 2 3

R> seq(from=1, to=3, by=1)
[1] 1 2 3

En el caso del conjunto de datos "USArrests", los datos se encuentran almacenados en cuatro vectores. Por ejemplo, donde se encuentran los homicidios, es un vector numérico:

  R> class(USArrests[,"Murder"])
  [1] "numeric"

El tamaño del vector donde se encuentran almacenados

  R> length(USArrests[,"Murder"])
  [1] 50

Y el primer elemento del vector es

  R> USArrests[,"Murder"][1]
  [1] 13.2

Las mediciones a nivel nominal o categorías se encuentran representadas por variables de tipo factor tales como la categoría a la cual pertenece un cierto tratamiento. Por ejemplo, el conjunto de datos "InsectSprays" representa el conteo de plagas en ciertas unidades agrícolas experimentales las cuales fueron tratadas con diferentes insecticidas nombrados aquí como "A", "B", "C", "D", "E" y "F". Veamos la estructura de este conjunto de datos:

  R> str(InsectSprays)
'data.frame':   72 obs. of  2 variables:
 $ count: num  10 7 20 14 14 12 10 23 17 20 ...
 $ spray: Factor w/ 6 levels "A","B","C","D",..: 1 1 1 1 1 1 1 1 1 1 ...

Los objectos de tipo factor y los de tipo character difieren en la forma en la que son almacenados internamente. Cada elemento de un vector de clase character se almacena como una variable character mientras que un objeto de tipo factor es almacenado como una variable entera indicando el nivel del factor. En otras palabras, el nombre de la categoría usada se convierte en un número entero que va del 1 al número de posibles valores de la categoría. En el caso de los insecticidas, el valor numérico de la variable spray corre del 1 al 6.

  R> nlevels(InsectSprays[,"spray"])
  [1] 6
  R> levels(InsectSprays[,"spray"])
  [1] "A" "B" "C" "D" "E" "F"

Como ejemplo de una estadística sencilla, se puede obtener un conteo de la frecuencia de los niveles de cada variable tipo factor a partir de

  R> table(InsectSprays[,"spray"])

 A  B  C  D  E  F
12 12 12 12 12 12

Importación y Exportación de Datos

En la sección anterior vimos como leer datos en línea y usar los paquetes de datos incluídos con el lenguaje. Ahora veremos formas más relevantes de como importar datos al sistema ${\textsf {R}}$. Los formatos más frecuentes que se encuentran en la práctica son las hojas de trabajo de Microsoft Excel, archivos de texto, y conexiones con bases de datos relacionales usando SQL. El importar datos directamente de una base de datos tal como Oracle o Microsoft SQL Server implica el conocimiento de SQL y por lo tanto dejaremos el tema para un capítulo posterior. Por lo pronto nos abocaremos a las hojas de trabajo y a archivos de texto. Los archivos de texto pueden venir en dos formas: campo fijo (donde cada elemento ocupa un lugar fijo de espacios) o campo variable delimitado por un carácter. Estos últimos por lo general vienen delimitados por una coma, un espacio, un tabulador, o un par de comillas. Supongamos que tenemos un archivo de texto el cual contiene información que queremos procesar. El archivo lleva el nombre de users.csv y es un archivo que contiene los datos de un grupo de usuarios. El archivo sólo contiene un identificador numérico, la edad del usuario, género, profesión, y el código postal de su domicilio. Los elementos se encuentran separados por comas y los datos de tipo texto se encuentran delimitados por comillas. Mostramos algunos renglones del archivo para obtener una idea de su estructura:

880,13,"M","student","83702"
881,39,"M","marketing","43017"
882,35,"M","engineer","40503"
883,49,"M","librarian","50266"
884,44,"M","engineer","55337"

Usaremos entonces la función read.table para leer el archivo:

R> users <- read.table("users.csv", header=TRUE, sep=",")

El argumento header=TRUE indica que la primera línea del archivo users.csv debe interpretarse como los nombres de las variables a leerse. Las columnas se encuentran separadas por una coma (sep=",") y los datos se almacenarán en el dataframe users. Se puede usar también la función read.csv para leer archivos separados por comas pero la función read.table es un poco más general y es mejor que nos acostumbremos a usarla. La función read.table trata de adivinar cuál es la clase de cada variable al leer el archivo. Veamos como fue estructurado dentro de ${\textsf {R}}$:

  R> str(users)
'data.frame':   943 obs. of  5 variables:
 $ UserID    : int  1 2 3 4 5 6 7 8 9 10 ...
 $ Age       : int  24 53 23 24 33 42 57 36 29 53 ...
 $ Gender    : Factor w/ 2 levels "F","M": 2 1 2 2 1 2 2 2 2 2 ...
 $ Occupation: Factor w/ 21 levels "administrator",..: 20 14 21 20 14 7 1 1 19 10 ...
 $ Zipcode   : Factor w/ 795 levels "00000","01002",..: 623 690 271 332 134 759 652 49 2 643 ...

Podemos observar que cada observación (cada línea del archivo) está compuesta por cinco variables: dos de clase int (enteros) y tres variables de tipo factor cada una con distinto número de niveles.

La exportación de datos desde ${\textsf {R}}$ se hace de forma similar. Se construye un archivo delimitado con comas para que pueda ser leído por Microsoft Excel a partir de un objeto tipo dataframe a través de:

  R> write.table(users, file="USERS.csv", sep=",", col.names=NA)

Si se desea guardar los datos para que solo sean procesados por ${\textsf {R}}$, cualquier objeto en ${\textsf {R}}$ se puede salvar en un archivo binario externo a través de:

  R> save(users, file="USERS.rda")

donde la extensión .rda es estándar. Si deseamos cargar los objetos de un archivo tipo rda a ${\textsf {R}}$.

 R> load("USERS.rda")

Manejo Básico de Datos

Los ejemplos mostrados en la sección anterior muestran la importancia de usar los dataframes para el manejo y almacenamiento de datos tabulares en ${\textsf {R}}$. Internamente, un dataframe es una lista de vectores con una longitud común $n$, donde $n$ es el número de filas en la tabla. Cada uno de estos vectores representa las medidas de una variable y hemos vista que se puede acceder tal variable por medio de su nombre. Por ejemplo, los tipos de profesiones en el juego de datos users.csv:

  R> occ <- users[,"Occupation"]
  R> occ[1]
  [1] technician

El vector occ es de clase Factor con 21 niveles y tiene una longitud de 943 elementos. Se puede extraer un subconjunto de los elementos del vector occ usando el operador []. Por ejemplo, para extraer los elementos 125 al 128, usamos el siguiente operador:

  R> occ[125:128]
  [1] lawyer    lawyer    none      marketing
  21 Levels: administrator artist doctor educator engineer ... writer

${\textsf {R}}$ nos recuerda que la variable que estamos leyendo es de tipo Factor y que posee 21 niveles que van desde "administrator" a "writer". Nota: El llamarles "niveles" a los posibles valores que puede tomar la variable Occupation no implica que exista una jerarquía entre ellos. Es simplemente un asunto de nomenclatura.

Una característica interesante de los índices en ${\textsf {R}}$ es el uso de índices negativos. El uso de éstos indica todos los elementos que no son parte del vector indicado por los corchetes. Por ejemplo, si queremos indicar los primeros cinco elementos del vector occ' lo podemos hacer de la siguiente forma usando índices negativos:

  R> occ[-(6:943)]
  [1] technician other      writer     technician other

Y en forma convencional con índices positivos:

  R> occ[1:5]
  [1] technician other      writer     technician other

Se puede imprimir una tabla de los primeros tres registros (filas) del juego de datos usando el hecho de que los dataframes manejan el concepto de filas y columnas. Entonces, separamos los subconjuntos correspondientes a las filas y a las columnas por medio de una coma. Por ejemplo:

  R> users[1:3, c("UserID", "Age", "Gender", "Occupation", "Zipcode")]
    UserID Age Gender Occupation Zipcode
  1      1  24      M technician   85711
  2      2  53      F      other   94043
  3      3  23      M     writer   32067

Con esta instrucción, le estamos diciendo a ${\textsf {R}}$ que imprima los tres primeros registros (filas) de los datos con las columnas indicadas. Al mismo tiempo, si quisiéramos extraer una sola variable del dataframe usaríamos:

  R> occ <- users$Occupation
  R> occ[1:5]
  [1] technician other      writer     technician other
  21 Levels: administrator artist doctor educator engineer ... writer

La instrucción mostrada arriba es equivalente a la instrucción

  R> occ <- users[,"Occupation"]

Ejemplo: Nos interesa saber cuáles son los datos correspondientes a cinco usuarios de mínima edad. Primero, ordenaremos las edades de menor a mayor, usando los índices correspondientes del vector Age para usarlos como índices en el vector users.

  R> ages <- order(users$Age)
  R> users[ages[1:3], c("UserID", "Age", "Gender", "Occupation", "Zipcode")]
    UserID Age Gender Occupation Zipcode
30      30   7      M    student   55436
471    471  10      M    student   77459
142    142  13      M      other   48118
609    609  13      F    student   55106
289    289  11      M       none   94619

La función order ordena el vector Age de menor a mayor y regresa los índices del vector en este nuevo orden. Usamos entonces estos índices para extraer los cinco elementos de users que poseen las cinco mínimas edades.

Otra forma de extraer datos es usar un vector lógico que adquiere el valor TRUE cuando selecciona el elemento deseado y FALSE cuando no. Por ejemplo, para obtener los datos de aquellos usuarios menores de 15 años, expresaríamos:

  R>  users[users$Age<15,  c("UserID", "Age", "Gender", "Occupation", "Zipcode")]
      UserID Age Gender Occupation Zipcode
30      30   7      M    student   55436
142    142  13      M      other   48118
206    206  14      F    student   53115
289    289  11      M       none   94619
471    471  10      M    student   77459
609    609  13      F    student   55106
628    628  13      M       none   94306
674    674  13      F    student   55337
813    813  14      F    student   02136
880    880  13      M    student   83702
887    887  14      F    student   27249
 R>

donde la expresion users$Age <15 indica un vector lógico de longitud 943 donde

R> table(users$Age<15)

FALSE  TRUE
  932    11

los elementos son falsos(FALSE) or verdaderos (TRUE).

Si tuviéramos datos faltantes, éstos estarían indicados por medio de un símbolo especial: NA indicando así que esa medición no se encuentra disponible. Si queremos saber cuáles elementos faltan en alguna columna usamos la siguiente función

 R> na_occ <- is.na(users$Occupation)
 R> table(na_occ)
na_occ
FALSE
  943

La función regresó con un vector lógico de 943 elementos todos ellos igual a FALSE. Lo que quiere decir que no hay elementos faltantes en la columna Occupation.

El tomar subconjuntos de los dataframes usando expresiones lógicas nos puede ahorrar mucho tecleado. La función subset toma un dataframe como su primer argumento, y una expresión lógica como el segundo. Por ejemplo, podemos seleccionar un subconjunto de todos los usuarios mayores de 40 años por medio de

 R> mayores50 <- subset(users, Age> 50)
 R> dim(mayores50)
[1] 105   5

El vector mayores50 es un vector de 105 filas con 5 columnas. Nótese que no es necesario especificar en el primer argumento la variable Age bajo consideración.

Cálculos con Datos

Algunas Estadísticas Sencillas

Dos funciones son muy útiles para extraer información general sobre los objetos en ${\textsf {R}}$: str y summary. La primera es para obtener información sobre el tipo y cantidad de los objetos que se encuentran almacenados en el dataframe; y la segunda es para obtener estadísticas de sumarización del conjunto de datos.

  R> summary(users)
     UserID           Age        Gender          Occupation     Zipcode
 Min.   :  1.0   Min.   : 7.00   F:273   student      :196   55414  :  9
 1st Qu.:236.5   1st Qu.:25.00   M:670   other        :105   55105  :  6
 Median :472.0   Median :31.00           educator     : 95   10003  :  5
 Mean   :472.0   Mean   :34.05           administrator: 79   20009  :  5
 3rd Qu.:707.5   3rd Qu.:43.00           engineer     : 67   55337  :  5
 Max.   :943.0   Max.   :73.00           programmer   : 66   27514  :  4
                                         (Other)      :335   (Other):909

Podemos ver que para la variable numérica Age, la función arrojó los valores extremos, el primer y tercer cuartil, la mediana y la media. Para el caso de las variables nominales Gender, Occupation, y Zipcode, la función calculó el número correspondiente a los niveles de la variable. Esto es, la variable Gender tiene 273 valores de F, y 670 de M. En el caso de las variables Occupation y Zipcode, la función proporciona los 6 valores más frecuentes.

Programando con R

En las secciones anteriores hemos estado haciendo todos nuestros cálculos con las funciones internas de ${\textsf {R}}$. Sin embargo, el verdadero poder de ${\textsf {R}}$ proviene de la expresividad de este lenguaje en forma de funciones que se usan para las tareas de análisis. Estudiaremos dos ejemplos pequeños a continuación.

Digamos que queremos agregar una medida robusta de variabilidad a las medidas de ubicación calculadas en la sección anterior. Aparte de la media y la mediana, queremos calcular el intervalo intercuartil. Esto es, la diferencia entre el tercer y el primer cuartil. Haciendo una búsqueda rápida del manual en línea (usando help("interquartile")) obtenemos el resultado de IQR, pero en lugar de usar directamente esta función, usaremos la función quartile para calcular cuantiles.

Una función en ${\textsf {R}}$ no es más que un objeto, y todos los objetos se crean de la misma manera. Esto es, tenemos que asignarle un objeto de tipo function a una variable. Un objeto function consiste de una lista de argumentos, valores por omisión, y un cuerpo que define los cálculos. El cuerpo empieza y termina con llaves {} y, por supuesto, se entiende que el cuerpo consiste de código en ${\textsf {R}}$

Regresando a nuestro ejemplo, le llamaremos a nuestra función iqr. La función iqr deberá operar sobre vectores numéricos, por lo que tendrá como argumento de entrada un vector x. Este vector numérico le será pasado a la función quantile para calcular los cuartiles en la muestra. La diferencia requerida entre el tercer y el primer cuartil se calculará usando diff. La definición de nuestra función se presenta a continuación:

R> iqr <- function(x) {
       q <- quantile(x, prob = c(0.25, 0.75), names = FALSE)
	return(diff(q))
}

Una prueba usando datos simulados provenientes de una distribución normal estándar muestra que nuestra función sí opera como debido. Si la comparamos con la función IQR mostramos que el resultado es el correcto

R> xdata <- rnorm(100)
R> iqr(xdata)
[1] 0.9977416
R> IQR(xdata)
[1] 0.9977416

Sin embargo, cuando el vector numérico contiene valores faltantes, la función falla como se muestra en el siguiente ejemplo:

R> xdata[1] <- NA
R> iqr(xdata)
Error in quantile.default(x, prob = c(0.25, 0.75), names = FALSE) :
  missing values and NaN's not allowed if 'na.rm' is FALSE

Para hacer la función un poco más flexible, sería conveniente agregar todos los argumentos de quantile a la lista de argumentos de iqr. Si lo hacemos por simple copia podría llevar a inconsistencias y errores como, por ejemplo, cuando la lista de argumentos de quantile cambia. En lugar de eso, usamos el argumento de puntos, el cual sirve de comodín para cualquier argumento de la función

R> iqr <- function(x, ...) {
+   q <- quantile(x, prob = c(0.25, 0.75), names=FALSE,
+                 ...)
+   return(diff(q))
+ }
R> iqr(xdata, na.rm=TRUE)
[1] 1.019106
R> IQR(xdata, na.rm=TRUE)
[1] 1.019106

El argumento na.rm=TRUE detecta si existen valores faltantes (NA) en el vector de entrada y los extrae para efectos del cálculo. El vector de entrada no es alterado al salir de la función.

Gráficas Sencillas

El grado de oblicuidad (skewness) de una distribución se puede obtener construyendo histogramas usando la función hist (Existen alternativas un poco más sofisticadas pero se verán en capítulos posteriores). Por ejemplo, el código para crear la figura 6 primero divide la región de graficando horizontalmente en dos partes igualmente espaciadas (layout) y luego grafica el histograma de las edades en la parte superior y el histograma del logaritmo de la edad en el segundo renglón. Este último parece ser un poco más

simétrico.

Histogramas de edades

Las relaciones bivariantes de dos variables continuas se grafican generalmente como una gráfica de puntos. En ${\textsf {R}}$, las relaciones de regresión se especifican a través de fórmulas modelo los cuales, en el caso de dos variables, pueden tomar una forma parecida a

  R> fm <- ingresos ~ ventas
  R> class(fm)
   [1] "formula"

donde la variable dependiente es la izquierda y la independiente la derecha (ventas). El tilde separa las partes izquierda y derecha. Esta fórmula se puede pasar como argumento a una función modelo como parte de un cálculo más complicado.

Supongamos que ahora deseamos las distribuciones de las edades pero separadas por género. Esto es, queremos las distribuciones de las edades masculinas y las femeninas por separado. El código ${\textsf {R}}$ para obtener una gráfica de cajas de las edades es:

R> tmp <- subset(users, Gender %in% c("M", "F"))
R> plot(Age ~ Gender, data = tmp, ylab="Edades", xlab="Genero", varwidth=TRUE)

Podemos observar el resultado en la figura 7

Gráfica de cajas (Boxplots) por edades

.

Organizando un Análisis

Aún cuando es posible realizar un análisis estadístico tecleando todos los mandos directamente en el apronte de ${\textsf {R}}$, es mucho más cómodo tener un archivo de texto por separado que contenga todos los pasos necesarios para llevar a cabo el análisis requerido. El archivo, que podemos llamar análisis.R, se puede correr dentro del ambiente de ${\textsf {R}}$ usando el mando source:

  R> source("analisis.R", echo = TRUE)

Si ponemos todos los mandos juntos necesarios para un análisis, podemos reproducir el análisis a voluntad con mayor número de datos o con distintos datos según la necesidad.

Algebra Matricial

Introducción

En un capítulo posterior veremos que una variable aleatoria multidimensional es simplemente un conjunto ordenado $[X_{1},\ldots ,X_{p}]$ de variantes sencillas. Cuando decimos "ordenado" queremos decir que la variante que describe la i-ésima faceta de cada unidad muestreada de la población siempre aparecerá en la i-ésima posición de la secuencia. El número de variantes $p$ en el arreglo será constante a través del problema o el análisis que se esté llevando a cabo. Por ejemplo, supóngase que los componentes de la variante son los pesos al nacer de los ratones hembra de una cierta cepa y los demás componentes son los pesos a los 10, 20 y 30 días de nacidas. Los pesos entonces son descritos por una variable aleatoria $[X_{1},X_{2},X_{3},X_{4}]$ con alguna distribución en el espacio de 4 dimensiones. Si esos pesos fueron registrados con una muestra de $N$ ratones, las observaciones se pueden resumir por medio del arreglo

${\begin{bmatrix}x_{11}&\cdots &x_{14}\\\cdots &\cdots &\cdots \\x_{N1}&\cdots &x_{N4}\end{bmatrix}}$

donde las filas corresponden a los distintos ratones. Tales arreglos rectangulares y lineales de números se conocen como matrices y vectores respectivamente, y las reglas para su manipulación se conocen como álgebra lineal. Esta álgebra es el lenguaje virtual del análisis multivariable, especialmente en la parte más común e importante la cuál se encuentra basada en la distribución normal multidimensional. Incluso, se puede afirmar que las técnicas de análisis multivariable no hubieran podido ser desarrolladas sin la facilidad, la conveniencia y la elegancia de las matrices.

En este capítulo se resumirán algunas propiedades, operaciones y teoremas del álgebra matricial necesarias para el material que se avecina. Si se desea profundizar en el tema, se proporcionarán referencias para consulta y aprendizaje.

Algunas Definiciones

Supongamos que tenemos a nuestra disposición una serie de elementos que se comportan de acuerdo a un conjunto de axiomas, por ejemplo, los números reales o complejos. Definimos una matriz

$A={\begin{bmatrix}a_{11}&\cdots &a_{1c}\\\cdots &\cdots &\cdots \\a_{r1}&\cdots &a_{rc}\end{bmatrix}}$

como un arreglo rectangular ordenado de los elementos. El término general de la matriz se expresará como $a_{ij}$ donde el primer subíndice se refiere a la i-ésima fila y el segundo subíndice se refiere a la j-ésima columna. Las dimensiones de una matriz son importantes, de modo que una matriz con $r$ filas y $c$ columnas será referida como una matriz de $r\times c$ dimensiones.

Para diferenciarlos de las matrices, les llamaremos a los números y variables unidimensionales que manejamos habitualmente, escalares. Un escalar es considerado como una matriz de $1\times 1$. En las primeras secciones de este capítulo supondremos que los elementos de las matrices son escalares reales. Posteriormente trataremos con matrices cuyos elementos son matrices de dimensiones más pequeñas.

Un vector es una matriz con una sola fila o columna. Normalmente, expresaremos el vector columna de $n$ componentes como

$x={\begin{bmatrix}x_{1}\\\vdots \\x_{n}\end{bmatrix}}$

De la misma forma, el vector fila lo expresaremos como

${\textbf {x}}^{\prime }=[x_{1},\ldots ,x_{n}]$

y consiste de una sola fila de $n$ elementos. Cada vector especifica las coordenadas de un punto en el espacio euclideano $n$-dimensional y la conexión con el concepto físico o analítico de un vector es aparente si se considera ese punto como el término de un segmento de recta que empieza en el origen de los ejes coordenados.

El apóstrofe que acompaña al símbolo del vector fila indica que ${\textbf {x}}^{\prime }$ es la transpuesta del vector columna x. En general, si A es cualquier matriz de $r\times c$ dimensiones, la transpuesta de A es la matriz ${\textbf {A}}^{\prime }$ de $c\times r$ dimensiones formada al intercambiar las filas y las columnas:

${\textbf {A}}^{\prime }={\begin{bmatrix}a_{11}&\cdots &a_{1c}\\\cdots &\cdots &\cdots \\a_{r1}&\cdots &a_{rc}\end{bmatrix}}={\begin{bmatrix}a_{11}&\cdots &a_{r1}\\\cdots &\cdots &\cdots \\a_{1c}&\cdots &a_{cr}\end{bmatrix}}$

Si una matriz es cuadrada e igual a su transpuesta, se dice que es simétrica. Entonces $a_{ij}=a_{ji}$ para todos los pares de $i$ y $j$. Por ejemplo

${\textbf {A}}={\begin{bmatrix}3&0&-1\\0&1&2\\-1&2&-4\end{bmatrix}}$

es simétrica, mientras que

${\textbf {B}}={\begin{bmatrix}2&3\\1&2\end{bmatrix}}$

no lo es. Por brevedad, se omitirán con frecuencia los elementos inferiores duplicados. Los elementos $a_{ii}$ de una matriz cuadrada ocupan lo que se conoce como las posiciones de la diagonal principal. La suma de estos elementos diagonales se conoce como la traza de A, y se denotará por

${\text{tr}}\;{\textbf {A}}=\sum _{i=1}^{n}a_{ii}$

Ciertas matrices son particularmente importantes. La matriz identidad

${\textbf {I}}={\begin{bmatrix}1&\cdots &0\\\cdots &\cdots &\cdots \\0&\cdots &1\end{bmatrix}}$

es una matriz cuadrada con unos en la diagonal principal y ceros en las posiciones restantes. La matriz diagonal de $p\times p$

${\textbf {D}}(a_{i})={\begin{bmatrix}a_{1}&\cdots &0\\\cdots &\cdots &\cdots \\0&\cdots &a_{p}\end{bmatrix}}$

tiene los elementos $a_{1},\ldots ,a_{p}$ en las posiciones de la diagonal principal y ceros en las posiciones restantes. Algunos de los elementos $a_{i}$ pueden ser cero. De aquí en adelante el término diag(A) = D($a_{ii}$) denotará la matriz diagonal formada por la matriz cuadrada A. Una matriz triangular de $p\times p$ tiene el patrón

${\begin{bmatrix}t_{11}&t_{12}&\cdots &t_{1p}\\0&t_{22}&\cdots &t_{2p}\\\cdots &\cdots &\cdots &\cdots \\0&0&\cdots &t_{pp}\end{bmatrix}}$

La matriz nula de $r\times c$

${\begin{bmatrix}0&\cdots &0\\\cdots &\cdots &\cdots \\0&\cdots &0\end{bmatrix}}$

tiene un cero en cada una de sus posiciones. A veces necesitaremos el vector

${\textbf {j}}^{\prime }=[1,\ldots ,1]$

y la matriz

${\textbf {E}}={\begin{bmatrix}1&\cdots &1\\\cdots &\cdots &\cdots \\1&\cdots &1\end{bmatrix}}$

la cual tiene la unidad en cada posición.

Operaciones Elementales con Matrices y Vectores

Las operaciones aritméticas escalares de suma, resta y multiplicación se pueden llevar a cabo si se siguen ciertas reglas. La división matricial es un poco más complicada y se dejará para una sección posterior. Igualdad. Dos matrices de $r\times c$ A y B se dice que son iguales si y sólo si

$a_{ij}=b_{ij}$

Suma. La suma de dos matrices de iguales dimensiones es la matriz de la suma de los elementos correspondientes. Si

${\textbf {A}}={\begin{bmatrix}a_{11}&\cdots &a_{1c}\\\cdots &\cdots &\cdots \\a_{r1}&\cdots &a_{rc}\end{bmatrix}}\;\;{\textbf {B}}={\begin{bmatrix}b_{11}&\cdots &b_{1c}\\\cdots &\cdots &\cdots \\b_{r1}&\cdots &b_{rc}\end{bmatrix}}$

entonces

${\textbf {A+B}}={\begin{bmatrix}a_{11}+b_{11}&\cdots &a_{1c}+b_{1c}\\\cdots &\cdots &\cdots \\a_{r1}+b_{r1}&\cdots &a_{rc}+b_{rc}\end{bmatrix}}$

Una matriz se puede sustraer de otra con iguales dimensiones formando la matriz de diferencias de los elementos individuales. Entonces

${\begin{aligned}\mathbf {A+B} &=\mathbf {B+A} \\\mathbf {A+(B+C)} &=\mathbf {(A+B)+C} \\\mathbf {A-(B-C)} &=\mathbf {A-B+C} \end{aligned}}$

Si las dimensiones de las matrices no son iguales, sus sumas o diferencias se encuentran indefinidas.

Multiplicación por un Escalar. La matriz A se multiplica por el escalar $c$ multiplicando cada elemento de A por $c$:

$c{\textbf {A}}={\begin{bmatrix}ca_{11}&\cdots &ca_{1k}\\\cdots &\cdots &\cdots \\ca_{r1}&\cdots &ca_{rk}\end{bmatrix}}$

Multiplicación de Matrices. Para que el producto de dos matrices AB esté definido es necesario que número de columnas de A sea igual al número de filas de B. Se dice entonces que las dimensiones de éstas matrices son conformables. Si A tiene dimensiones de $p\times r$, y B tiene dimensiones de $r\times q$, entonces el $ij$-ésimo elemento del producto ${\textbf {C=AB}}$ se calcula como

$c_{ij}=\sum _{k=1}^{r}a_{ik}b_{kj}$

Esta es la suma de los productos de los elementos correspondientes a la i-ésima fila de A y la j-ésima columna de B. Las dimensiones de AB son $p\times q$. Por ejemplo, si

$\mathbf {A} ={\begin{bmatrix}1&2&3\\-1&0&1\end{bmatrix}}\;\;\;\mathbf {B} ={\begin{bmatrix}6&5&4\\-1&1&-1\\0&2&0\end{bmatrix}}$

entonces

$\mathbf {AB} ={\begin{bmatrix}1(6)+2(-1)+3(0)&1(5)+2(1)+3(2)&1(4)+2(-1)+3(0)\\-1(6)+0(-1)+1(0)&-1(5)+0(1)+1(2)&-1(4)+0(-1)+1(0)\end{bmatrix}}$

El producto BA se encuentra indefinido, ya que las dos filas de A no se conforman a las tres columnas de B.

Las leyes distributivas y asociativas se aplican a la multiplicación de matrices:

${\begin{aligned}\mathbf {A(B+C)} &=\mathbf {AB+AC} \\\mathbf {A(BC)} &=\mathbf {AB(C)} \end{aligned}}$

Sin embargo, la ley conmutativa no se puede aplicar a la multiplicación matricial ya que en general, no es cierto que AB = BA. Por esto, el orden de la multiplicación es crucial y nos referiremos a la premultiplicación por A en el producto AB o la posmultiplicación por B en el producto AB. Por ejemplo, sea

$\mathbf {A} ={\begin{bmatrix}1&2\\2&3\end{bmatrix}}\;\;\mathbf {B} ={\begin{bmatrix}3&-1\\-1&1\end{bmatrix}}$

Entonces

$\mathbf {AB} ={\begin{bmatrix}1&1\\3&1\end{bmatrix}}\;\;\mathbf {BA} ={\begin{bmatrix}1&3\\1&1\end{bmatrix}}$

Ni los productos son iguales ni se conserva la simetría de las matrices originales en la multiplicación.

La multiplicación de una matriz por una matriz identidad conformable deja a la matriz inmutada. La premultiplicación de una matriz por la matriz diagonal con elementos $d_{1},\ldots ,d_{r}$ posee el efecto de multiplicar cada elemento en la i-ésima fila por $d_{i}$:

$\mathbf {D} (d_{i})\mathbf {A} ={\begin{bmatrix}d_{1}a_{11}&\cdots &d_{1}a_{1c}\\\cdots &\cdots &\cdots \\d_{r}a_{r1}&\cdots &d_{r}a_{rc}\end{bmatrix}}$

La posmultiplicación por una matrix diagonal similar de $c\times c$ multiplica cada elemento en la j-ésima columna por $d_{j}$.

Una matriz puede considerarse como una transformación lineal de los vectores en un espacio a vectores en otro espacio. Si x tiene $m$ componentes y y tiene $n$ componentes, es posible expresar una transformación del sistema de coordenadas en $m$ dimensiones de los elementos de x al espacio $n$-dimensional de y en forma matricial como

${\begin{bmatrix}y_{1}\\\vdots \\y_{n}\end{bmatrix}}={\begin{bmatrix}a_{11}&\cdots &a_{1m}\\\cdots &\cdots &\cdots \\a_{n1}&\cdots &a_{nm}\end{bmatrix}}{\begin{bmatrix}x_{1}\\\vdots \\x_{m}\end{bmatrix}}=\mathbf {Ax}$

La transformación a un tercer conjunto de variables especificado por el vector ${\textbf {z}}$ de $p$ componentes puede representarse como

${\begin{aligned}\mathbf {z} &=\mathbf {By} \\&=\mathbf {BAx} \end{aligned}}$

y así el producto de dos o más matrices se puede considerar como la matriz del resultado de una serie de transformaciones lineales.

Ejemplo: Muchas veces es necesario en estudios e investigaciones el convertir varias mediciones recabadas de los sujetos a una escala con un origen e unidades comunes. Si $x_{ij}$ es la medición del $i$-ésimo individuo en $ji$-ésima característica, y ${\bar {x}}_{j}$ y $s_{j}$ son la media y la desviación estándar de la muestra respectivamente, una transformación común que puede usarse es la medida z: $z_{ij}={\dfrac {z_{ij}-{\bar {x}}_{j}}{s_{j}}}$ Las transformaciones observadas se pueden calcular con facilidad a través de algunas operaciones matriciales sencillas. Si expresamos las mediciones originales como una matriz de $N\times p$

${\textbf {X}}={\begin{bmatrix}x_{11}&\cdots &x_{1p}\\\cdots &\cdots &\cdots \\x_{N1}&\cdots &x_{Np}\end{bmatrix}}$

y formamos la matriz diagonal

${\textbf {D}}={\begin{bmatrix}{\dfrac {1}{s_{1}}}&\cdots &0\\\cdots &\cdots &\cdots \\0&\cdots &{\dfrac {1}{s_{p}}}\end{bmatrix}}$

a partir de las desviaciones estándard de la muestra. Si introducimos la matriz de $N\times N$

${\textbf {E}}={\begin{bmatrix}1&\cdots &1\\\cdots &\cdots &\cdots \\1&\cdots &1\end{bmatrix}}$

con un uno en cada posición, entonces la matriz Z de $N\times p$ dimensiones correspondiente a las mediciones normalizadas se puede calcular como

${\textbf {Z}}=\left({\textbf {I}}-{\frac {1}{N}}{\textbf {E}}\right){\textbf {XD}}$

El resultado de posmultiplicar E por X es obtener la suma de cada columna de la matriz de datos.

Productos Internos Vectoriales. El producto interno de dos vectores con el mismo número de elementos se define como la suma de los productos de los elementos correspondientes:

${\textbf {x}}^{\prime }{\textbf {y}}=[x_{1},\ldots ,x_{p}]{\begin{bmatrix}y_{1}\\\vdots \\y_{p}\end{bmatrix}}$

Ya que el producto interno es escalar, entonces ${\textbf {x}}^{\prime }{\textbf {y}}={\textbf {y}}^{\prime }{\textbf {x}}$. El producto interno de x consigo mismo es la suma de los cuadrados de los elementos de x.

Los productos internos tienen importantes interpretaciones geométricas. El producto interno de x consigo mismo se conoce como la longitud cuadrada de x, ya que es el cuadrado de la distancia desde el origen al punto especificado por los elementos de x en el sistema de coordenadas en $p$ dimensiones. En general, la distancia entre los puntos con coordenadas dadas por x y y es

$d=\left[\sum _{i=1}^{p}(x_{i}-y_{i})^{2}\right]^{\frac {1}{2}}$

El coseno $\theta$ del ángulo entre los vectores x y y es

$\cos \theta ={\frac {{\textbf {x}}^{\prime }{\textbf {y}}}{({\textbf {x}}^{\prime }{\textbf {x}})^{\frac {1}{2}}({\textbf {y}}^{\prime }{\textbf {y}})^{\frac {1}{2}}}}$

Esta división de vectores entre sus respectivas longitudes se conoce como normalización y es fácil ver que a medida que los vectores normalizados van coincidiendo, su producto interno tiende a la unidad. Al mismo tiempo, aquellos vectores que se encuentran perpendiculares tienen un producto interno igual a cero.

Ejemplo. En el espacio tridimensional, los vectores

${\begin{aligned}{\textbf {x}}^{\prime }&=[1,0,0]\\{\textbf {y}}^{\prime }&=[0,1,0]\\{\textbf {z}}^{\prime }&=[0,0,1]\end{aligned}}$

se pueden interpretar como si fueran los ejes coordenados. El vector

${\textbf {e}}^{\prime }=[1,1,1]$

hace un ángulo de ${\sqrt {3}}/3$ con ${\textbf {x}}^{\prime },{\textbf {y}}^{\prime }$ y ${\textbf {z}}^{\prime }$ y caracteriza la línea equiangular en el espacio tridimensional. Las cantidades ${\sqrt {3}}/3$ se conocen como los cosenos direccionales} de la línea. Similarmente, los vectores

${\begin{aligned}{\textbf {u}}^{\prime }&=[1,1,0]\\{\textbf {v}}^{\prime }&=[1,-1,0]\end{aligned}}$

forman un ángulo de 90 grados en el plano $xy$ del espacio. Los vectores ${\textbf {u}}^{\prime }$ y ${\textbf {v}}^{\prime }$ forman cada uno un ángulo de 45 grados con ${\textbf {x}}^{\prime }$, y sus ángulos con ${\textbf {y}}^{\prime }$ son 45 y -45 grados, respectivamente.

Transpuesta de un Producto Matricial. Si la matriz A es conformable para posmultiplicación por otra matriz B, se puede verificar fácilmente que la transpuesta de su producto es igual al producto de sus transpuestas tomadas en el orden inverso:

$({\textbf {AB}})^{\prime }={\textbf {B}}^{\prime }{\textbf {A}}^{\prime }$

En general, si ${\textbf {A}}_{1},\ldots ,{\textbf {A}}_{k}$ son conformables,

$({\textbf {A}}_{1},\cdots {\textbf {A}}_{k})^{\prime }={\textbf {A}}_{k}^{\prime }\cdots {\textbf {A}}_{1}^{\prime }$

Estas propiedades van a ser muy necesarias en capítulos posteriores.

El Determinante de una Matriz Cuadrada

Con toda matriz cuadrada se encuentra asociada una cantidad escalar llamada el determinante. Un determinante generalizado para una matriz $A$ tiene el valor

$\det {\textbf {A}}\equiv |{\textbf {A}}|=\sum _{i=1}^{k}a_{ij}C_{ij}$

donde se supone que no hay suma sobre el índice $j$ y donde $C_{ij}$ es el cofactor de $a_{ij}$ definido por

$C_{ij}\equiv (-1)^{i+j}M_{ij}$

y $M_{ij}$ es el Error al representar (Error de conversión. El servidor ("https://wikimedia.org/api/rest_") informó: "Cannot get mml. TeX parse error: Undefined control sequence \emph"): {\displaystyle {\emph {menor}}} de la matriz ${\textbf {A}}$ formado al eliminar la fila $i$ y la columna $j$ de ${\textbf {A}}$. Este proceso se conocido como la expansión Laplaciana por menores o también conocido como Desarrollo Laplaciano. Un determinante también puede ser calculado expresando todas las permutaciones de $(1,\ldots ,n)$ tomando cada permutación como los subíndices de las letras $a,b,\ldots$ y sumando con los signos dados por $\epsilon _{p}=(-1)^{i(p)}$ donde $i(p)$ es el número de inversiones de la permutación en la permutación $p$, y $\epsilon _{n_{1}n_{2}}$ es el símbolo de permutación. Por ejemplo, para $n=3$, las permutaciones y el número de inversiones que contienen son 123 (0 inversiones), 132 (1 inversion), 213(1), 231(2), 312(2), y 321 (3), de modo que el determinante de una matriz de $3\times 3$ es

${\begin{vmatrix}a_{1}&a_{2}&a_{3}\\b_{1}&b_{2}&b_{3}\\c_{1}&c_{2}&c_{3}\\\end{vmatrix}}=a_{1}b_{2}c_{3}-a_{1}b_{3}c_{2}-a_{2}b_{1}c_{3}+a_{2}b_{3}c_{1}+a_{3}b_{1}c_{2}-a_{3}b_{2}c_{1}.$

Para determinantes mayores de $3\times 3$ calcularlos a mano es oneroso y existen métodos numéricos comprobados que los calculan. Mostraremos más adelante en el capítulo dedicado a programación en ${\textsf {R}}$ como calcular determinantes.

Ciertas propiedades de los determinantes nos serán útiles en los capítulos subsiguientes:

1. El determinante de una matriz diagonal es simplemente el producto de los elementos de la diagonal. Igualmente, el determinante de una matriz triangular es el producto de los elementos diagonales.
2. Si los elementos de una sola fila o columna de la matriz \textbf{A} de $n\times n$ son multiplicados por el escalar $c$, el determinante de la nueva matriz es igual a $c|{\textbf {A}}|$. Si cada elemento es multiplicado por $c$, entonces $|c{\textbf {A}}|=c^{n}|{\textbf {A}}$.
3. Si dos filas o columnas de una matriz cuadrada son intercambiadas, entonces el signo del determinante se invierte.
4. Si dos columnas o dos filas de la matriz son iguales, entonces el determinante debe ser cero. Por lo tanto, las filas o columnas proporcionales de una matriz indican un determinante igual a cero.
5. El determinante de una matriz no cambia agregando un múltiplo de alguna columna a otra columna. Lo mismo vale para las filas.
6. Si todos los elementos de una fila o una columna son iguales a cero, entonces el determinante es cero.
7. Si ${\textbf {A}}$ y ${\textbf {B}}$ son matrices de $n\times n$, el determinante de ${\textbf {AB}}$ es igual al producto de los determinantes individuales.
8. La suma de los productos de los elementos de una fila en una matriz cuadrada con los cofactores correspondientes de una fila distinta es igual a cero. Lo mismo se aplica a las columnas.

La Matriz Inversa

Ahora ya podemos definir la operación matricial correspondiente a la división escalar. La inversa de la matriz cuadrada ${\textbf {A}}$ es la matriz única ${\textbf {A}}^{-1}$ con elementos tales que

${\textbf {AA}}^{-1}={\textbf {A}}^{-1}{\textbf {A}}={\textbf {I}}$

Es posible que ${\textbf {A}}^{-1}$ no exista, en la misma forma que no es posible dividir entre cero con escalares. Entonces se dice que ${\textbf {A}}$ es Error al representar (Error de conversión. El servidor ("https://wikimedia.org/api/rest_") informó: "Cannot get mml. TeX parse error: Undefined control sequence \emph"): {\displaystyle {\emph {singular}}} . Las matrices cuyas inversas existen se conocen como Error al representar (Error de conversión. El servidor ("https://wikimedia.org/api/rest_") informó: "Cannot get mml. upstream connect error or disconnect/reset before headers. reset reason: connection termination"): {\displaystyle {\emph {nosingulares}}} .

Los elementos de ${\textbf {A}}^{-1}$ pueden ser calculadas a partir de dos resultados de la sección anterior. Se forma la matriz de cofactores $C_{ij}$:

${\textbf {C}}={\begin{bmatrix}C_{11}&\cdots &C_{1n}\\\cdots &\cdots &\cdots \\C_{n1}&\cdots &C_{nn}\end{bmatrix}}$

conocida como la Error al representar (Error de conversión. El servidor ("https://wikimedia.org/api/rest_") informó: "Cannot get mml. TeX parse error: Undefined control sequence \emph"): {\displaystyle {\emph {adjunta}}} de ${\textbf {A}}$ Entonces el producto interno de la $i$-ésima fila de ${\textbf {C}}$ y la $h$-ésima fila de ${\textbf {A}}$ es igual a $|{\textbf {A}}|$ si $i=j$ y es igual a cero si $i\neq h$. Si tomamos la transpuesta de ${\textbf {C}}$ y dividimos cada elemento entre $|{\textbf {A}}|$, tendremos entonces la inversa deseada:

${\textbf {A}}^{-1}={\begin{bmatrix}{\frac {1}{|{\textbf {A}}|}}C_{11}&\cdots &{\frac {1}{|{\textbf {A}}|}}C_{n1}\\\cdots &\cdots &\cdots \\{\frac {1}{|{\textbf {A}}|}}C_{1n}&\cdots &{\frac {1}{|{\textbf {A}}|}}C_{nn}\end{bmatrix}}$

Es evidente de esta última expresión que la matriz inversa existe si y sólo si el determinante de ${\textbf {A}}$ es distinto de cero. El cálculo de ${\textbf {A}}^{-1}$ por medio de cofactores es muy ineficiente para muchas aplicaciones prácticas y se considerarán otros métodos en la siguiente sección.

Se usarán con frecuencia las siguientes propiedades de las inversas:

1. La inversa de una matriz simétrica también es simétrica.
2. La inversa de la transpuesta de ${\textbf {A}}$ es la transpuesta de ${\textbf {A}}^{-1}$.
3. La inversa del producto de varias matrices cuadradas en igual al producto de las inversas en el orden contrario: $({\textbf {ABC}})^{-1}={\textbf {C}}^{-1}{\textbf {B}}^{-1}{\textbf {A}}^{-1}$
4. Si $c$ es un escalar distinto de cero, $(c{\textbf {A}})^{-1}=1/c{\textbf {A}}^{-1}$.

Ejemplo: La inversa de

${\textbf {A}}={\begin{bmatrix}2&0&0&0\\0&4&0&0\\0&0&3&2\\0&0&1&2\end{bmatrix}}$

es

${\textbf {A}}^{-1}={\begin{bmatrix}{\frac {1}{2}}&0&0&0\\0&{\frac {1}{4}}&0&0\\0&0&{\frac {1}{2}}&-{\frac {1}{2}}\\0&0&-{\frac {1}{4}}&{\frac {3}{4}}\end{bmatrix}}$

El Rango de una Matriz

Dos vectores de $p$ componentes se dice que son linealmente dependientes si los elementos de un vector son proporcionales a los del segundo. Por lo tanto, los vectores fila ${\textbf {x}}^{\prime }=[1,0,-1]\;\;\;{\textbf {y}}^{\prime }=[4,0,-4]$ son linealmente dependientes, mientras que ${\textbf {u}}^{\prime }=[2,-1,0,7]\;\;\;{\textbf {v}}^{\prime }=[6,2,0,0]$ son linealmente independientes. Un conjunto de $k$ vectores equidimensionales se conoce como linealmente independiente si no es posible expresar cualquier vector del conjunto como alguna combinación lineal de los vectores restantes. Esto es, los vectores ${\textbf {x}}_{1},\ldots ,{\textbf {x}}_{k}$ forman un conjunto linealmente independiente si no existen escalares $c_{1},\ldots ,c_{k}$ tales que para algun vector ${\textbf {x}}_{i}$ en el conjunto se da que $c_{i}{\textbf {x}}_{i}=c_{1}{\textbf {x}}_{1}+\cdots +c_{i-1}{\textbf {x}}_{i-1}+c_{i+1}{\textbf {x}}_{i+1}+\cdots +c_{k}{\textbf {x}}_{k}\;.$ Los vectores ${\textbf {x}}^{\prime }=[1,-1,2]\;\;\;{\textbf {y}}^{\prime }=[2,0,-1]\;\;\;{\textbf {z}}^{\prime }=[0,-2,5]$ constituyen un conjunto linealmente dependiente ya que ${\textbf {z}}^{\prime }=2{\textbf {x}}^{\prime }-{\textbf {y}}^{\prime }$, mientras que los vectores unitarios ${\textbf {t}}^{\prime }=[1,0,0]\;\;\;{\textbf {u}}^{\prime }=[0,1,0]\;\;\;{\textbf {v}}^{\prime }=[0,0,1]$ son linealmente independientes.

Ejemplo: Supóngase que una cierta prueba de habilidades cognoscitivas consiste de seis subescalas. Las tres primeras miden ciertas habilidades verbales, y las calificaciones se suman para obtener lo que se llama la medida verbal. Las otras tres reflejan ciertas habilidades motrices y espaciales y el total se conoce como la medida de desempeño. Finalmente, la suma de las seis pruebas se usa como una medida general de habilidad cognoscitiva. La prueba le ha sido proporcionada a un cierto número de individuos y se ha propuesto que las nueve medidas se usen para tratar de relacionar la inteligencia con otras medidas obtenidas de los individuos. Sin embargo, si cada muestra de observaciones sobre una medida se considera como un vector de $N$ componentes, se puede ver de inmediato que sólo seis de los vectores son linealmente independientes y que no se ha obtenido información nueva alguna de las tres medidas derivadas. Se verá más adelante que el incluír éstas composiciones lineales pueden impedir algunos tipos de análisis estadísticos.

Formalicemos ahora el grado de independencia lineal en un conjunto de vectores como el rango de una matriz formada por los vectores. Supóngase por el momente que el número $k$ de vectores no excede la dimensionalidad $p$ de los vectores. Entonces el rango de

${\textbf {X}}={\begin{bmatrix}{\textbf {x}}^{\prime }\\\vdots \\{\textbf {x}}_{k}^{\prime }\end{bmatrix}}$

es el número de vectores fila linealmente independientes de la matriz. El rango puede variar desde cero para una matriz nula, hasta $k$ para una matriz de rango completo. Si, por otro lado, el número de filas de ${\textbf {X}}$ es mayor que el número de columnas, el rango de \textbf{X} será el número de columnas linealmente independientes. En cualquiera de los casos, se puede mostrar que el rango de una matriz es un número único, sin importar si viene del cálculo con las filas o las columnas. Tenemos entonces las siguiente propiedad: La matriz ${\textbf {A}}$ tiene un rango igual a r si contiene al menos un menor distinto de cero de $r\times r$ dimensiones y ningún menor distinto de cero con dimensionalidad mayor que r.

Ejemplo: La matriz

${\begin{bmatrix}1&2&3&4\\3&6&9&12\\4&3&2&1\\-1&3&7&11\\8&6&4&2\end{bmatrix}}$

tiene un rango de dos, ya que se puede mostrar que la segunda fila es igual a tres veces la primera, la fila cuatro es igual a la diferencia de las filas dos y tres, y la fila cinco es igual a dos veces la fila tres. Las únicas filas linealmente independientes son la uno y la tres, y se dice que forman una base de los cinco vectores fila.

El rango tiene estas importantes propiedades:

1. El rango de ${\textbf {A}}^{\prime }$ es igual al rango de ${\textbf {A}}$. Esto viene de la equivalencia de las definiciones de rango para filas y columnas.
2. El rango de ${\textbf {A}}^{\prime }{\textbf {A}}$ es igual al de ${\textbf {A}}$. Similarmente, el rango de ${\textbf {A}}{\textbf {A}}^{\prime }$ es equivalente al rango de ${\textbf {A}}$.
3. El rango de ${\textbf {A}}$ no cambia si se pre- o posmultiplica ${\textbf {A}}$ por una matriz no singular.

Operaciones Elementales de Filas y Columnas

El rango de una matriz no cambia al realizar las siguientes operaciones elementales de filas y columnas:

1. Intercambio de dos filas o dos columnas.
2. La multiplicación de todos los elementos de una fila o una columna por una constante escalar.
3. La suma de una fila o columna cuyos elementos han sido multiplicados por una constante escalar a otra fila o columna.

Estas transformaciones se pueden representar por matrices no singulares. Las operaciones de fila se llevan a cabo por premultiplicación de la matriz dada por la matriz de la operación elemental, mientras que las operaciones de columna se llevan a cabo por posmultiplicación. Por la propiedad 3 del párrafo anterior tales transformaciones no singulares dejan al rango invariante. Eligiendo la secuencia apropiada de operaciones de fila o columna se puede reducir la matriz a una forma que consiste solamente de columnas o filas linealmente independientes, y así se determina su rango.

Ejemplo: Las matrices

${\textbf {E}}_{1}={\begin{bmatrix}0&1&0\\1&0&0\\0&0&1\end{bmatrix}}\;\;\;{\textbf {E}}_{2}={\begin{bmatrix}c&0&0\\0&1&0\\0&0&1\end{bmatrix}}\;\;\;{\textbf {E}}_{3}={\begin{bmatrix}1&0&0\\d&1&0\\0&0&1\end{bmatrix}}$

se pueden usar para llevar a cabo transformaciones elementales sobre las primeras dos filas de

${\textbf {A}}={\begin{bmatrix}5&1\\-2&3\\3&2\end{bmatrix}}$

dando así

${\textbf {E}}_{1}{\textbf {A}}={\begin{bmatrix}-2&3\\5&1\\3&2\end{bmatrix}}\;\;\;{\textbf {E}}_{2}{\textbf {A}}={\begin{bmatrix}5c&c\\-2&3\\3&2\end{bmatrix}}\;\;\;{\textbf {E}}_{3}{\textbf {A}}={\begin{bmatrix}5&1\\5d-2&d+3\\3&2\end{bmatrix}}$

Usemos ahora las operaciones elementales de fila para calcular el rango de la matriz

${\begin{bmatrix}1&0&2\\-3&2&1\\-1&2&5\\3&0&6\end{bmatrix}}$

Nótese que la última fila es igual a la primera fila multiplicada por 3. La fila tres es igual a la fila dos mas dos veces la fila uno. Operaciones sucesivas de fila apropiadas reducen $A$ en la siguiente forma:

${\begin{bmatrix}1&0&2\\-3&2&1\\-1&2&5\\0&0&0\end{bmatrix}}\;\;\;{\begin{bmatrix}1&0&2\\-3&2&1\\0&0&0\\0&0&0\end{bmatrix}}$

Ya que las filas restantes son linealmente independientes, se concluye que $A$ es de rango igual a dos.

Ambas operaciones de fila y columna se pueden aplicar simultáneamente a cualquier matriz de $m\times n$ para reducirla a su forma canónica

${\textbf {F}}={\begin{bmatrix}1&\cdots &0&\cdots &0\\\cdots &\cdots &\cdots &\cdots &\cdots \\0&\cdots &1&\cdots &0\\\cdots &\cdots &\cdots &\cdots &\cdots \\0&\cdots &0&\cdots &\cdots \end{bmatrix}}$

la cual contiene un uno en las primeras $r$ posiciones diagonales y cero en todas las demás. Esto quiere decir que para cada matriz $A$ existen matrices elementales de transformación de filas y columnas llamadas ${\textbf {R}}_{1}\ldots ,{\textbf {R}}_{p},{\textbf {C}}_{1},\ldots ,{\textbf {C}}_{q}$ tales que

${\begin{aligned}{\textbf {R}}_{p}\cdots {\textbf {R}}_{1}{\textbf {A}}{\textbf {C}}_{1}\cdots {\textbf {C}}_{q}&={\textbf {PAQ}}\\&={\textbf {F}}\end{aligned}}$

es la forma canónica de $A$. Pero si ${\textbf {A}}$ es cuadrada y de rango completo $n$, entonces las operaciones de fila producirán

${\textbf {R}}_{p}\cdots {\textbf {R}}_{1}{\textbf {A}}={\textbf {I}}$

esto significa que

${\textbf {A}}^{-1}={\textbf {R}}_{p}\cdots {\textbf {R}}_{1}{\textbf {I}}$

Por lo tanto, la misma secuencia de operaciones elementales de fila que reduce $A$ a la matriz identidad transforma la matriz identidad de $n\times n$ a la inversa de $A$. Esta propiedad de las transformaciones elementales nos proporciona con una herramienta poderosa para invertir matrices: simplemente se expresa la matriz de $n\times 2n$

${\begin{bmatrix}{\textbf {A}}&&{\textbf {I}}\end{bmatrix}}$

y se aplican operaciones elementales a cada fila con el fin de transformar las primeras $n$ columnas en la matriz identidad. Cuando esto ha sido logrado, aparecerá ${\textbf {A}}^{-1}$ en las últimas $n$ columnas.

Ejemplo: La transformaciones sucesivas en la inversión de

${\textbf {A}}={\begin{bmatrix}1&1&2\\2&3&3\\3&-1&8\end{bmatrix}}$

se muestran a continuación:

${\begin{bmatrix}1&1&2&1&0&0\\2&3&3&0&1&0\\3&-1&8&0&0&1\end{bmatrix}}$
${\begin{bmatrix}1&1&2&1&0&0\\0&1&-1&-2&1&0\\0&-4&2&-3&0&1\end{bmatrix}}$
${\begin{bmatrix}1&0&3&3&-1&0\\0&1&-1&-2&1&0\\0&0&-2&-11&4&1\end{bmatrix}}$
${\begin{bmatrix}1&0&0&-13.5&5&1.5\\0&1&0&3.5&-1&-0.5\\0&0&1&5.5&-2&-0.5\end{bmatrix}}$

y

${\textbf {A}}^{-1}={\begin{bmatrix}-13.5&5&1.5\\3.5&-1&-0.5\\5.5&-2&-0.5\end{bmatrix}}$

El determinante de una matriz también puede ser calculado por medio de operaciones elementales de filas o columnas ya que por la propiedad 5 de la sección 2.4 estas transformaciones no alteran su valor. Este proceso, llamado condensación pivotal, consiste en reducir la matriz a forma triangular, de donde el determinante se puede calcular como el producto de los elementos diagonales. Si se usan operaciones de fila, la matriz se expresa con el elemento más grande de la primera columna en la posición (1,1). Ese elemento se usa como el pivote para reducir a cero los elementos restantes de la columna. Estos reposicionamientos y reducciones se continúan hasta que se llega a la forma triangular de la matriz. Se elige el elemento más granda por cuestiones de exactitud numérica.

Ejemplo: Vamos a evaluar el determinante de la matriz

${\textbf {A}}={\begin{bmatrix}10&5&-2&0\\6&3&2&1\\4&5&12&3\\5&1&3&8\end{bmatrix}}$

Empezamos restando múltiplos apropriados de la primera fila de las filas restantes:

${\begin{bmatrix}10&5&-2&0\\0&0&3.2&1\\0&3&12.8&3\\0&-1.5&4&8\end{bmatrix}}$

Ya que el elemento (2,2) es igual a cero, intercambiamos la segunda y tercera fila de la matriz, y compensamos el cambio de signo del determinante multiplicando la nueva segunda fila por -1:

${\begin{bmatrix}10&5&-2&0\\0&-3&-12.8&-3\\0&0&3.2&1\\0&-1.5&4&8\end{bmatrix}}$

Una operación de fila pivoteando en la segunda columna transforma la matriz anterior en:

${\begin{bmatrix}10&5&-2&0\\0&-3&-12.8&-3\\0&0&3.2&1\\0&0&10.4&9.5\end{bmatrix}}$

Finalmente tenemos

${\begin{bmatrix}10&5&-2&0\\0&-3&-12.8&-3\\0&0&3.2&1\\0&0&0&6.25\end{bmatrix}}$

y el determinante es:

${\begin{aligned}|{\textbf {A}}|&=(10)(-3)(3.2)(6.25)\\&=-600\end{aligned}}$

Matrices Inversas Generalizadas

Se han propuesto varias generalizaciones de la matriz inversa para matrices rectangulares de cualquier rango. Una inversa generalizada que es útil para resolver sistemas de ecuaciones lineales es la matriz ${\textbf {G}}$ la cual satisface

${\textbf {AGA}}={\textbf {A}}$

donde ${\textbf {A}}$ es de $p\times q$ dimensiones y de rango $r$, mientras que ${\textbf {G}}$ tiene que ser necesariamente de $q\times p$ dimensiones y del mismo rango que ${\textbf {A}}$. Existen otros tipos de inversas generalizadas, las que satisfacen las siguientes condiciones:

${\textbf {GAG}}={\textbf {G}}\;\;\;({\textbf {GA}})^{\prime }={\textbf {GA}}\;\;\;({\textbf {AG}})^{\prime }={\textbf {AG}}$

Nuestra matriz ${\textbf {G}}$ se puede calcular a partir de la reducción canónica de ${\textbf {A}}$. Expresamos

${\textbf {PAQ}}={\begin{bmatrix}{\textbf {D}}&{\textbf {0}}\\{\textbf {0}}&{\textbf {0}}\end{bmatrix}}$

donde ${\textbf {D}}$ es una matriz diagonal de $r\times r$ (que no es necesariamente la identidad) y las matrices nulas tienen las dimensiones apropiadas. Las matrices ${\textbf {P}}$ y ${\textbf {Q}}$ son los productos respectivos definidos anteriormente. Presentamos ahora la matriz ${\textbf {F}}^{-}$ de $q\times p$:

${\textbf {F}}^{-}={\begin{bmatrix}{\textbf {D}}^{-1}&{\textbf {0}}\\{\textbf {0}}&{\textbf {0}}\end{bmatrix}}$

Entonces ${\textbf {G}}={\textbf {QF}}^{-}{\textbf {P}}$ es una inversa generalizada de la matriz ${\textbf {A}}$ ya que uno puede verificarlo reemplazando ${\textbf {A}}$ en la expresión por su representación ${\textbf {P}}^{-1}{\textbf {FQ}}^{-1}$. La matriz ${\textbf {G}}$ no es única. a menos que ${\textbf {A}}$ sea cuadrada y no singular.

${\textbf {Ejemplo:}}$ Vamos a calcular la inversa generalizada de la matriz ${\textbf {A}}$ vista en un ejemplo anterior. Las matrices de operadores de fila que reducen la tercera y la cuarta fila a los vectores nulos son:

${\textbf {R}}_{1}={\begin{bmatrix}1&0&0&0\\0&1&0&0\\0&0&1&0\\-3&0&0&1\end{bmatrix}}\;\;\;\;{\textbf {R}}_{2}={\begin{bmatrix}1&0&0&0\\0&1&0&0\\-2&-1&1&0\\0&0&0&1\end{bmatrix}}$

Entonces

${\textbf {P}}={\textbf {R}}_{2}{\textbf {R}}_{1}={\begin{bmatrix}1&0&0&0\\0&1&0&0\\-2&-1&1&0\\-3&0&0&1\end{bmatrix}}$

Las transformaciones de columna que reducen a ${\textbf {PA}}$ a la forma diagonal canónica son

${\textbf {C}}_{1}={\begin{bmatrix}1&0&0\\{\frac {3}{2}}&1&-{\frac {1}{2}}\\0&0&1\end{bmatrix}}\;\;\;\;{\textbf {C}}_{2}={\begin{bmatrix}1&0&-2\\0&1&0\\0&0&1\end{bmatrix}}$

y

${\textbf {Q}}={\textbf {C}}_{1}{\textbf {C}}_{2}={\begin{bmatrix}1&0&-2\\{\frac {3}{2}}&1&-{\frac {7}{2}}\\0&0&1\end{bmatrix}}$

Entonces

${\textbf {F}}={\textbf {PAQ}}={\begin{bmatrix}1&0&0\\0&2&0\\0&0&0\\0&0&0\end{bmatrix}}\;\;\;\;{\textbf {F}}^{-}={\begin{bmatrix}1&0&0&0\\0&{\frac {1}{2}}&0&0\\0&0&0&0\end{bmatrix}}$

y

${\textbf {G}}={\begin{bmatrix}1&0&0&0\\{\frac {3}{2}}&{\frac {1}{2}}&0&0\\0&0&0&0\end{bmatrix}}$

es una inversa generalizada de ${\textbf {A}}$.

Ecuaciones Lineales Simultáneas

El conjunto de ecuaciones con las incógnitas $x_{1},\ldots ,x_{n}$

${\begin{aligned}a_{11}x_{1}+\cdots +a_{1n}x_{n}&=c_{1}\\\vdots \\a_{n1}x_{1}+\cdots +a_{nn}x_{n}&=c_{n}\end{aligned}}$

se conoce como un sistema de ecuaciones lineales simultáneas en n incógnitas o sistema lineal y se puede expresar en forma compacta usando matrices como

${\textbf {Ax}}={\textbf {c}}$

donde ${\textbf {A}}$ es la matriz de coeficientes de $m\times n$, ${\textbf {x}}^{\prime }=[x_{1},\ldots ,x_{n}]$, y ${\textbf {c}}^{\prime }=[c_{1},\ldots ,c_{m}]$ El sistema general ${\textbf {Ax=c}}$ de $m$ ecuaciones en $n$ incógnitas tiene una solución si y solo si la matriz aumentada de $m\times n$ $\left[{\textbf {A}}\;\;{\textbf {c}}\right]$ tiene el mismo rango $r$ que ${\textbf {A}}$. Si no es así, entonces se dice que el sistema es inconsistente. Es esencial distinguir entre sistemas homogéneos, para los cuales ${\textbf {c}}=0$, y sistemas no homogéneos. Se considerarán tres tipos de ecuaciones que son particularmente relevantes:

Sistemas No Homogéneos: A cuadrada y no singular

Debido a que ${\textbf {A}}^{-1}$ existe, la solución única del sistema es

${\textbf {x}}={\textbf {A}}^{-1}{\textbf {c}}$

Esto viene de premultiplicar ambos lados de la ecuación por ${\textbf {A}}^{-1}$.

${\textbf {Ejemplo:}}$ La matriz inversa del sistema

${\begin{aligned}x_{1}+x_{2}-x_{3}&=1\\-x_{1}+x_{2}+x_{3}&=-1\\x_{1}-x_{2}+x_{3}&=1\end{aligned}}$

es

${\begin{bmatrix}{\frac {1}{2}}&0&{\frac {1}{2}}\\{\frac {1}{2}}&{\frac {1}{2}}&0\\0&{\frac {1}{2}}&{\frac {1}{2}}\\\end{bmatrix}}$

La solución es

${\begin{aligned}x_{1}&={\frac {1}{2}}(1)+0(-1)+{\frac {1}{2}}(1)\\&=1\\x_{2}&={\frac {1}{2}}(1)+{\frac {1}{2}}(-1)+0(1)\\&=0\\x_{3}&=0(1)+{\frac {1}{2}}(-1)+{\frac {1}{2}}(1)\\&=0\end{aligned}}$

Sistemas No Homogéneos: A $m\times n$ y Rango $r$

Si se satisface la condición sobre los rangos de las matrices aumentadas y de coeficientes, se puede encontrar una solución seleccionando cualesquier $r$ ecuaciones linealmente independientes y obteniendo $r$ de las $n$ incógnitas en términos de las constantes $c_{i}$ y de las $n-r$ variables restantes.

${\textbf {Ejemplo:}}$ Las ecuaciones

${\begin{aligned}x_{1}+x_{2}&=2\\x-1+x_{2}&=1\end{aligned}}$

se pueden interpretar geométricamente como dos líneas paralelas en el plano $x_{1}x_{2}$. Debido a que no intersectan, el sistema no puede tener solución. El rango de la matriz de coeficientes

${\begin{bmatrix}1&1\\1&1\end{bmatrix}}$

es uno, mientras que el de la matriz aumentada

${\begin{bmatrix}1&1&2\\1&1&1\end{bmatrix}}$

es dos.

${\textbf {Ejemplo:}}$ El sistema

${\begin{aligned}2x_{1}+3x_{2}-x_{3}&=1\\x_{2}+x_{3}&=2\end{aligned}}$

tiene matrices aumentadas y de coeficientes con rango igual a dos. Se expresan las ecuaciones como

${\begin{aligned}2x_{1}+3x_{2}-x_{3}&=1\\x_{2}&=2-x_{3}\end{aligned}}$

y se obtiene $x_{1}$ por sustitución. La solución general está dada por el vector

$[{\frac {1}{2}}(-5+4x_{3}),2-x_{3},x_{3}]$

donde $x_{3}$ puede tomar cualquier valor.

Sistemas Homogéneos de Ecuaciones

No es posible que un sistema homogéneo de ecuaciones lineales sea inconsistente, ya que el rango de $[{\textbf {A}}\;{\textbf {0}}]$ es el mismo que el de ${\textbf {A.}}$ Todos los sistemas homogéneos se ven satisfechos por la solución trivial ${\textbf {x}}^{\prime }=[0,\ldots ,0]$, y si el rango de la matriz de coeficientes es igual al número de incógnitas $n$, ésta será la única solución. Por lo tanto un sistema homogéneo tentra una solución no trivial si y sólo si el rango $r$ de ${\textbf {A}}$ es estrictamente menor que $n$, y siempre será posible encontrar $n-r$ soluciones linealmente independientes del sistema tal que cualquier combinación lineal de éstas soluciones es también una solución. En aplicaciones estadísticas se encontrarán frecuentemente sistemas que contienen tantas incógnitas como ecuaciones, tales sistemas tienen soluciones no triviales si y sólo si los determinantes de los coeficientes desaparecen.

${\textbf {Ejemplo:}}$ El sistema

${\begin{aligned}3x_{1}-2x_{2}&=0\\5x_{1}+x_{2}&=0\end{aligned}}$

define dos líneas en el plano $x_{1}x_{2}$ que intersectan en el origen. Por lo tanto la única solución es la trivial: $[0,0]$.

${\textbf {Ejemplo:}}$ En el sistema

${\begin{aligned}x_{1}-x_{2}+2x_{3}&=0\\x_{1}+3x_{2}-2x_{3}&=0\\3x_{1}+x_{2}+2x_{3}&=0\end{aligned}}$

la tercera ecuación es igual a la suma de dos veces la primera mas la segunda, y entonces el rango es igual a dos. La única solución independiente se puede determinar en unidades de $x_{3}$ resolviendo el sistema

${\begin{aligned}x_{1}-x_{2}&=-2x_{3}\\x_{1}+3x_{2}&=2x_{3}\end{aligned}}$

La solución general del conjunto original es, en unidades de $x_{3}$

${\begin{aligned}x_{1}&=-x_{3}\\x_{2}&=x_{3}\end{aligned}}$

${\textbf {Ejemplo:}}$ El sistema

${\begin{aligned}x_{1}+x_{2}-x_{3}-x_{4}&=0\\-x_{2}+3x_{3}+x_{4}&=0\\x_{1}-x_{2}+2x_{3}&=0\\3x-1+x_{2}-2x_{4}&=0\end{aligned}}$

también es de rango dos. La primera y tercera ecuaciones se pueden resolver convenientemente en términos de $x_{3}$ y $x_{4}$ para dar el vector solución

$[{\frac {1}{2}}(-x_{3}+x_{4}),{\frac {1}{2}}(3x_{3}+x_{4}),x_{3},x_{4}]$

Si hacemos $x_{3}$ igual a cero y $x_{4}$igual a uno, y alternativamente hacemos $x_{3}$ igual a uno y $x_{4}$ igual a cero, tendremos los siguientes vectores solución linealmente independientes:

$[{\frac {1}{2}},{\frac {1}{2}},0,1]$
$[-{\frac {1}{2}},{\frac {3}{2}},1,0]$

Ya que cualquier combinación lineal de éstas soluciones es también una solución, la solución más general al sistema original es

${\begin{aligned}x_{1}&={\frac {1}{2}}(a-b)\\x_{2}&={\frac {1}{2}}(a+3b)\\x_{3}&=b\\x_{4}&=a\end{aligned}}$

para $a$ y $b$ arbitrarios.

Solución por medio de inversas generalizadas

Se ha demostrado que el sistema ${\textbf {Ax}}={\textbf {c}}$ de $m$ ecuaciones consistentes en $n$ incógnitas tiene la solución

${\textbf {x}}^{*}={\textbf {Gc}}+({\textbf {GA-I}}){\textbf {z}}$

donde ${\textbf {G}}$ es la inversa generalizada de ${\textbf {A}}$ definida por , ${\textbf {I}}$ es la matriz identidad de $n\times n$, y ${\textbf {z}}$ es cualquier vector de $n\times 1$ constantes arbitrarias. Hacemos notar que el término que contine a ${\textbf {z}}$ desaparece cuando ${\textbf {G}}={\textbf {A}}^{-1}$.

${\textbf {Ejemplo:}}$ Obtendremos soluciones a los sistemas de los ejemplos X y Y por el método de las inversas generalizadas. Para el primer sistema, ${\textbf {P}}={\textbf {I}}$ y

${\textbf {Q}}={\begin{bmatrix}1&-3/2&2\\0&1&-1\\0&0&1\end{bmatrix}}\;\;\;\;{\textbf {F}}={\begin{bmatrix}2&0&0\\0&1&0\end{bmatrix}}$
${\textbf {G}}={\begin{bmatrix}1/2&-3/2\\0&1\\0&0\end{bmatrix}}\;\;\;\;{\textbf {GA}}={\begin{bmatrix}1&0&-2\\0&1&1\\0&0&0\end{bmatrix}}$

La solución más general es

${\textbf {x}}^{\prime }=[-1/2(5+z_{3}),2+z_{3},-z_{3}]$

Para el ejemplo Y

${\textbf {P}}={\begin{bmatrix}1&0&0\\-1&1&0\\-2&-1&1\end{bmatrix}}\;\;\;\;{\textbf {Q}}={\begin{bmatrix}1&1&-1\\0&1&1\\0&0&1\end{bmatrix}}\;\;\;\;{\textbf {F}}^{-}={\begin{bmatrix}1&0&0\\0&1/4&0\\0&0&0\end{bmatrix}}$
${\textbf {P}}={\begin{bmatrix}3/4&1/4&0\\-1/4&1/4&0\\0&0&0\end{bmatrix}}\;\;\;\;{\textbf {GA}}={\begin{bmatrix}1&0&1\\0&1&-1\\0&0&0\end{bmatrix}}$

Si hacemos que $z_{3}$ sea igual a $-x_{3}$la solución general es ${\textbf {x}}^{\prime }=[-x_{3},x_{3},x_{3}]$.

Métodos Numéricos para Sistemas Lineales

Existen también una variedad de métodos computacionales iterativos y matriciales para obtener la solución de un sistema lineal. Ya que una cátedra sobre métodos numéricos para sistemas lineales está fuera del alcance de este texto, no abundaremos sobre el tema. Referiremos al lector a alguno de los bien conocidos textos sobre análisis numérico.

Vectores y Matrices Ortogonales

En la sección mostramos la expresión

$\cos \theta ={\frac {{\textbf {x}}^{\prime }{\textbf {y}}}{({\textbf {x}}^{\prime }{\textbf {x}})^{\frac {1}{2}}({\textbf {y}}^{\prime }{\textbf {y}})^{\frac {1}{2}}}}$

para el coseno del ángulo entre los vectores ${\textbf {x}}$ y ${\textbf {y}}$ y observamos que si su producto interno desaparecía, entonces los vectores formaban un ángulo recto. Tales vectores se llaman ortogonales u ortonormales si sus longitudes han sido normalizadas a la unidad. Una matriz ortogonal ${\textbf {T}}$ es una matriz cuadrada cuyas files son un conjunto de vectores ortonormales. Por lo tanto

${\begin{aligned}{\textbf {TT}}^{\prime }&={\textbf {T}}^{\prime }{\textbf {T}}\\&={\textbf {I}}\end{aligned}}$

y la inversa de ${\textbf {T}}$ es simplemente su transpuesta ${\textbf {T}}^{\prime }$. Una matriz ortogonal sencilla de $2\times 2$ es

${\begin{bmatrix}0&1\\1&0\end{bmatrix}}$

Esta es la matriz de transformación de una reflexión de los puntos en el plano $xy$ alrededor de la línea de 45 grados. Una matriz ortogonal más general de esa dimensión es

${\textbf {T}}={\begin{bmatrix}\cos \theta &\sin \theta \\-\sin \theta &\cos \theta \end{bmatrix}}$

y se interpreta como la matriz de transformación para la rotación de los ejes coordenados $xy$ en un ángulo $\theta$. Esto es, si $x$ y $y$ fueran las coordenadas de un punto bajo los ejes originales, el punto después de haber girado los ejes, tendría las coordenadas

${\begin{aligned}u&=x\cos \theta +y\sin \theta \\v&=-x\sin \theta +y\cos \theta \end{aligned}}$

donde los nuevos ejes coordenados han sido girados un ángulo de $\theta$ grados. Se pueden construír matrices ortogonales más grandes por medio del proceso de ortogonalización Gram-Schmidt empezando con una fila normalizada y armando la matriz según el requisito de mutua ortogonalidad de filas. Por ejemplo, la llamada matriz Helmert

${\textbf {T}}={\begin{bmatrix}{\frac {1}{\sqrt {3}}}&{\frac {1}{\sqrt {3}}}&{\frac {1}{\sqrt {3}}}\\{\frac {1}{\sqrt {2}}}&-{\frac {1}{\sqrt {2}}}&0\\{\frac {1}{\sqrt {6}}}&{\frac {1}{\sqrt {6}}}&-{\frac {2}{\sqrt {6}}}\end{bmatrix}}$

es una matriz ortogonal de $3\times 3$.

La matriz ortogonal ${\textbf {T}}$ de $3\times 3$ se interpreta como la matriz de una transformación lineal equivalente a una rotación rígida o a una rotación seguida de una reflexión de los $n$ ejes coordenados alrededor del origen. Las distancias en este espacio permanecen invariantes ante esta rotación, ya que si ejecutamos la transformación

${\textbf {y}}={\textbf {Tx}}$

entonces

${\begin{aligned}{\textbf {x}}^{\prime }{\textbf {x}}&={\textbf {y}}^{\prime }{\textbf {T}}^{\prime }{\textbf {Ty}}\\&={\textbf {y}}^{\prime }{\textbf {Iy}}\end{aligned}}$

y las suma de los cuadrados será invariante. Las matrices ortogonales tienen estas útiles propiedades:

1. Las columnas de una matriz ortogonal son ortogonales.
2. El determinante de una matriz ortogonal siempre será 1 ó -1. Ya que $|{\textbf {TT}}^{\prime }|=|{\textbf {I}}|=1$ y $|{\textbf {T}}|=|{\textbf {T}}^{\prime }|$, $|{\textbf {T}}|=\pm 1\,$.
3. El producto de matrices ortogonales de la misma dimensión es también ortogonal. Esto es, una secuencia de rotaciones rígidas y reflexiones de los ejes coordenados se puede expresar como una sola rotación rígida con las reflexiones apropiadas

Formas Cuadráticas

Una forma cuadrática en las variables $x_{1},\ldots ,x_{n}$ es una expresión del tipo

$f(x_{1},\ldots ,x_{n})=a_{11}x_{1}^{2}+a_{22}x_{2}^{2}+\cdots +a_{nn}x_{n}^{2}+2a_{12}x_{1}x_{2}+\cdots +2a_{2n}x_{1}x_{n}+\cdots +2a_{n-1,n}x_{n-1}x_{n}=\sum _{i=1}^{n}\sum _{j=1}^{n}a_{ij}x_{i}x_{j}$

donde $a_{ij}=a_{ji}$. Algunas de las $a_{ij}$ pueden ser iguales a cero. Hacemos notar de inmediato que la forma cuadrática se puede expresar en notación matricial como ${\textbf {x}}^{\prime }{\textbf {Ax}}$ donde ${\textbf {x}}^{\prime }=[x_{1},\ldots ,x_{n}]$y ${\textbf {A}}$ es la matriz simétrica de $n\times n$ de coeficientes. La forma cuadrática más sencilla es simplemente $f(x)=a_{11}x^{2}$, la ecuación de una parábola en la variable $x$. Las formas cuadráticas en donde las $x_{i}$ son variables aleatorias desempeñan un papel muy importante en la estadística de una y muchas variables. Por ejemplo, la suma del cuadrado de las desviaciones alrededor de la media muestral

$\sum _{i=1}^{N}(x_{i}-{\bar {x}})^{2}=\sum _{i=1}^{N}x_{i}^{2}-{\frac {1}{N}}\left(\sum x_{i}\right)^{2}$

se puede expresar como una forma cuadrática en términos de las observaciones $x_{i}$ usando la matriz

${\textbf {A}}={\begin{bmatrix}{\frac {N-1}{N}}&{\frac {-1}{N}}&\cdots &{\frac {-1}{N}}\\{\frac {-1}{N}}&{\frac {N-1}{N}}&\cdots &{\frac {-1}{N}}\\\cdots &\cdots &\cdots &\cdots \\{\frac {-1}{N}}&{\frac {-1}{N}}&\cdots &{\frac {N-1}{N}}\end{bmatrix}}$

Una matriz simétrica ${\textbf {A}}$ y su forma simétrica asociada se llaman definidas positivas si ${\textbf {x}}^{\prime }{\textbf {Ax}}>0$ para toda ${\textbf {x}}$ no nula. Si ${\textbf {x}}^{\prime }{\textbf {Ax}}\leq 0$ la forma y su matriz se conocen como positivas semidefinidas. Las definiciones son similares para las formas negativas definidas y las formas negativas semidefinidas. Una forma indefinida puede tomar valores positivos, negativos o cero. Vamos a necesitar en los capítulos posteriores las propiedades especiales de las matrices positivas definidas y semidefinidas como también de las formas cuadráticas.

Las formas positivas definidas tienen matrices de rango completo. Es posible, a través de completar el cuadrado sucesivamente, el reducir tales formas en $n$ variables a la forma

$d_{1}y_{1}^{2}+\cdots +d_{n}y_{n}^{2}$

la cual solo contiene cuadrados de las nuevas variables $y_{i}$ y con coeficientes $d_{i}>0$. En forma similar, una forma positiva semidefinida cuadrática se puede reducir a

$d_{1}y_{1}^{2}+\cdots +d_{r}y_{r}^{2}$

donde todos los coeficientes son positivos, y $r\leq n$ es el rango de la matriz de la forma. Sin embargo, estas propiedades no son medios convenientes para determinar la naturaleza de la forma y procederemos a declarar una condición necesaria y suficiente para la definitez o semidefinitez positiva. A partir de la matriz, se forma la secuencia de determinantes menores principales

$p_{0}=1\;\;\;\;p_{1}=a_{11}\;\;\;\;p_{2}={\begin{vmatrix}a_{11}&&a_{12}\\a_{12}&&a_{22}\end{vmatrix}}\;\;\;\;p_{i}={\begin{vmatrix}a_{11}&\cdots &a_{1i}\\\cdots &\cdots &\cdots \\a_{1i}&\cdots &a_{ii}\end{vmatrix}}\;\;\;\;\cdots \;\;\;\;p_{n}=|{\textbf {A}}|$

Si ${\textbf {A}}$ es de rango $r$ se dice que es regular si $p_{r}\neq 0$ y no hay $p_{i}$ consecutivos que sean cero. Siempre será posible reexpresar cualquier matriz simétrica en forma regular intercambiando las filas y, al mismo tiempo, las columnas correspondientes. Entonces, si ${\textbf {A}}$ es una matriz arreglada regularmente,

1. Una condición necesaria y suficiente para la definitez positiva es que $p_{i}>0$ para todo $i=1,\ldots ,n$ .
2. Una condición necesaria y suficiente para la semidefinitez positiva es que $p_{1}>0,\ldots ,p_{r}>0$ y los $n-r\,\,p_{i}$ restantes sean iguales a cero, donde $r$ puede ser igual a $n$.

${\textbf {Ejemplo:}}$ La matriz

${\begin{bmatrix}4&6&0\\6&9&0\\0&0&0\end{bmatrix}}$

tiene la secuencia de determinantes principales menores $p_{0}=1$, $p_{1}=4$, $p_{2}=0$,$p_{3}=0$, y por eso no se encuentra regularmente arreglada. Podemos reordenar las filas y las columnas para obtener la nueva matriz:

${\begin{bmatrix}2&0&0\\0&4&6\\0&6&9\end{bmatrix}}$

la cual tiene la secuencia de determinantes $p_{0}=1$, $p_{1}=2$, $p_{2}=8$, $p_{3}=0$. La matriz y su forma cuadrática son positivas semidefinidas con rango igual a dos.

${\textbf {Ejemplo:}}$ La matriz de la ecuación X que es la suma de cuadrados de la muestra es un caso particular de una matriz en patrón con una elemento común diagonal $a$ y elementos $b$ que se encuentran fuera de la diagonal. El determinante de una matriz así de $i\times i$ dimensiones es

$(a-b)^{i-1}[a+(i-1)b]$

de tal forma que los determinantes menores principales de son

$p_{i}={\frac {1}{N}}[N-1-(i-1)]$

Los primeros $N-1$ determinantes son positivos, mientras que el $N$-ésimo siempre será cero. Entonces la suma de las desviaciones cuadradas es una forma cuadrática positiva semidefinida de rango $N-1$.

Las raíces y vectores característicos de una matriz

Las raíces características de la matriz ${\textbf {A}}$ de $p\times p$ son las soluciones de la ecuación determinante

$|{\textbf {A}}-\lambda {\textbf {I}}|=0$

El determinante es un polinomio de $p$-ésimo grado en la variable $\lambda$, y por eso ${\textbf {A}}$ tiene $p$ raíces características. El desarrollo de Laplace del determinante característico nos permite expresar el polinomio característico como

$|{\textbf {A}}-\lambda {\textbf {I}}|=(-\lambda )^{p}+S_{1}(-\lambda )^{p-1}+S_{2}(-\lambda )^{p-2}+\cdots -S_{p-1}\lambda +|{\textbf {A}}|$

donde $S_{i}$ es la suma de todos los $i\times i$ determinantes menores. $S_{1}$ es simplemente la suma de los elementos diagonales de ${\textbf {A}}$, o ${\text{tr}}{\textbf {A}}$. De la teoría de ecuaciones polinomiales se puede obtener que:

1. El producto de las raíces características de ${\textbf {A}}$ es igual a $|{\textbf {A}}|$.
2. La suma de las raíces características de ${\textbf {A}}$ es igual a la traza de ${\textbf {A}}$.

${\textbf {Ejemplo:}}$ La matriz

${\textbf {A}}={\begin{bmatrix}2&1&1\\1&2&1\\1&1&2\end{bmatrix}}$

tiene una traza de 6 y tres determinantes menores principales de $2\times 2$ cada uno igual a

${\begin{vmatrix}2&1\\1&2\end{vmatrix}}=3$

mientras que $|{\textbf {A}}|=4$. La ecuación característica de ${\textbf {A}}$ es

$\lambda ^{3}-6\lambda ^{2}+9\lambda -4=0$

y sus raíces son 1, 1, y 4.

De aquí en adelante usaremos frecuentemente las siguientes propiedades de las raíces características:

1. Las raíces características de una matriz simétrica con elementos reales son todas reales.
2. Las raíces características de una matriz positiva definida son todas positivas.
3. Si una matriz de simétrica de $n\times n$ es positiva semidefinida con rango $r$, contiene exactamente $r$ raíces características positivas y $n-r$ raíces cero.
4. Las raíces características distintas de cero del producto ${\textbf {AB}}$ son iguales a las raíces distintas de cero del producto ${\textbf {BA}}$. Por lo tanto las trazas de ${\textbf {AB}}$ y ${\textbf {BA}}$ son iguales.
5. Las raíces características de una matriz diagonal son los elementos diagonales mismos.

Con cada raíz característica $\lambda _{i}$ de la matriz cuadrada ${\textbf {A}}$ se encuentra asociado un vector característico ${\textbf {x}}_{i}$ cuyos elementos satisfacen el sistema homogéneo de ecuaciones

$\left[{\textbf {A}}-\lambda _{i}{\textbf {I}}\right]{\textbf {x}}_{i}={\textbf {0}}$

Por definición de la raíz característica el determinante del sistema desaparece, y siempre existe una solución no trivial ${\textbf {x}}_{i}$. Muchos de los vectores característicos que encontraremos en el futuro serán calculados a partir de matrices simétricas. Las raíces y vectores característicos de tales matrices tienen las siguientes propiedades:

1. Si $\lambda _{i}$ y $\lambda _{j}$ son raíces características distintas de la matriz simétrica ${\textbf {A}}$, sus vectores asociados ${\textbf {x}}_{i}$ y ${\textbf {x}}_{j}$ serán ortogonales. Esto es evidente si premultiplicamos las definiciones de los vectores

${\textbf {Ax}}_{i}=\lambda _{i}{\textbf {x}}_{i}\;\;\;\;\;{\textbf {Ax}}_{j}=\lambda _{j}{\textbf {x}}_{j}$

por ${\textbf {x}}_{j}^{\prime }$ y ${\textbf {x}}_{i}^{\prime }$ respectivamente. Pero esto implica que

$\lambda _{i}{\textbf {x}}_{j}^{\prime }{\textbf {x}}_{i}=\lambda _{j}{\textbf {x}}_{i}^{\prime }{\textbf {x}}_{j}$

y ya que $\lambda _{i}\neq \lambda _{j}$ podemos concluír que los vectores son ortogonales.

1. Para cada matriz simétrica real ${\textbf {A}}$ existe una matriz ortogonal ${\textbf {P}}$ tal que

${\textbf {P}}^{\prime }{\textbf {AP}}={\textbf {D}}$

donde ${\textbf {D}}$ es la matriz diagonal de las raíces características de ${\textbf {A}}$. Los vectores característicos normalizados de ${\textbf {A}}$ pueden tomarse como las columnas de ${\textbf {P}}$. Aún si las raíces características no son distintas, puede ser posible seleccionar los elementos de sus vectores para proporcionar un conjunto mutuamente ortogonal de vectores característicos.

Estas propiedades de las matrices simétricas tienen implicaciones importantes para las formas cuadráticas. Si aplicamos la transformación ortogonal

${\textbf {x}}={\textbf {Py}}$

a las $p$ variables en la forma cuadrática ${\textbf {x}}^{\prime }{\textbf {Ax}}$, entonces la forma se convierte en

${\begin{aligned}{\textbf {x}}^{\prime }{\textbf {Ax}}&={\textbf {y}}^{\prime }{\textbf {P}}^{\prime }{\textbf {APy}}\\&={\textbf {y}}^{\prime }{\textbf {Dy}}\\&=\lambda _{1}y_{i}^{2}+\cdots +\lambda _{r}y_{r}^{2}\end{aligned}}$

donde la $\lambda _{i}$ son las raíces características de la matriz de coeficientes, y $r$ es el rango de la forma. Cualquier forma cuadrática real puede reducirse a una suma pesada de cuadrados calculando los vectores y raíces característicos de su matriz.

${\textbf {Ejemplo:}}$ El vector característico de la raíz más grande del ejemplo anterior debe satisfacer las ecuaciones

${\begin{aligned}-2x_{11}+x_{12}+x_{13}&=0\\x_{11}-2x_{12}+x_{13}&=0\\x_{11}+x_{12}-2x_{13}&=0\end{aligned}}$

Si fijamos $x_{13}$ arbitrariamente igual a uno y resolvemos las primeras dos ecuaciones no homogéneas, el vector característico será

${\textbf {x}}^{\prime }=[1,1,1]$

y por supuesto cualquier vector no nulo $[a,a,a]$ es un vector característico para $\lambda _{3}=4$. El vector asociado con la doble raíz $\lambda _{1}=\lambda _{2}=1$ debe satisfacer el sistema con la ecuación linealmente independiente

$x_{21}+x_{22}+x_{23}=0$

La solución más general a esto es el vector

$[a,b,-a-b]$

y seleccionando adecuadamente $a$ y $b$ se pueden generar dos vectores característicos linealmente independientes para la raíz doble. Nótese que cualquier elección de $a$ y $b$ dará vectores ortogonales al primer vector $[1,1,1]$. Podemos elegir $a$y $b$ para continuar esta ortogonalidad. Por ejemplo,

${\begin{aligned}{\textbf {x}}_{1}^{\prime }&=[0,1,-1]\\{\textbf {x}}_{2}^{\prime }&=[-2,1,1]\end{aligned}}$

${\textbf {Ejemplo:}}$ Las matrices cuadradas con la propiedad de

${\textbf {AA}}={\textbf {A}}$

se llaman idempotentes y desempeñan un papel sumamente importante en la teoría del análisis de varianza. Las raíces características de una matriz idempotente son cero o uno, y una forma cuadrática con una matriz idempotente puede reducirse a una suma de $r$ términos cuadrados. Es fácil verificar que la matriz es idempotente y que es posible realizar la transformación

$\sum _{i=1}^{N}(x_{i}-{\bar {x}})^{2}=y_{1}^{2}+\cdots +y_{N-1}^{2}$

a una suma de los cuadrados de $N-1$ variables aleatorias nuevas cuya independencia estadística proviene de la ortogonalidad de la transformación.

Afortunadamente, hoy en día cualquier software matemático incluye rutinas para el cálculo de las raíces y vectores característicos de matrices por lo que no entraremos en el tema de métodos numéricos para éstos cálculos. Sin embargo, las rutinas se cubrirán detenidamente en el capítulo correspondiente al trato de matrices por medio del lenguaje ${\textsf {R}}$.

Particiones de matrices

Con frecuencia es necesario o conveniente agrupar columnas o filas de una matriz debido a que sus variables comparten ciertas características. Tal matriz se puede expresar como un arreglo

${\textbf {A}}={\begin{bmatrix}{\textbf {A}}_{11}&\cdots &{\textbf {A}}_{1n}\\\cdots &\cdots &\cdots \\{\textbf {A}}_{m1}&\cdots &{\textbf {A}}_{mn}\end{bmatrix}}$

de submatrices ${\textbf {A}}_{ij}$ que contienen $r_{i}$ filas y $c_{j}$ columnas, donde es obvio que todas las submatrices en una fila dada de ${\textbf {A}}$ deben tener el mismo número de filas y cada columna debe estar compuesta de matrices con el mismo número de columnas. Las operaciones con matrices partidas en esta forma son similares a las ordinarias, excepto que se debe tomar en cuenta la naturaleza no escalar de los participantes. La suma de las matrices partidas ${\textbf {A}}$ y ${\textbf {B}}$ cuyas submatrices tienen dimensiones similares es

${\textbf {A}}+{\textbf {B}}={\begin{bmatrix}{\textbf {A}}_{11}+{\textbf {B}}_{11}&\cdots &{\textbf {A}}_{1n}+{\textbf {B}}_{1n}\\\cdots &\cdots &\cdots \\{\textbf {A}}_{m1}+{\textbf {B}}_{m1}&\cdots &{\textbf {A}}_{mn}+{\textbf {B}}_{mn}\end{bmatrix}}$

El producto de las matrices partidas ${\textbf {A}}$ y ${\textbf {B}}$ es

${\textbf {AB}}={\begin{bmatrix}\sum _{j=1}^{n}{\textbf {A}}_{1j}{\textbf {B}}_{j1}&\cdots &\sum _{j=1}^{n}{\textbf {A}}_{1j}{\textbf {B}}_{jp}\\\cdots &\cdots &\cdots \\\sum _{j=1}^{n}{\textbf {A}}_{mj}{\textbf {B}}_{j1}&\cdots &\sum _{j=1}^{n}{\textbf {A}}_{mj}{\textbf {B}}_{jp}\end{bmatrix}}$

Las dimensiones dentro de cada producto de las submatrices debe conformarse; si las submatrices de ${\textbf {A}}$ tienen el número de columnas $c_{1},\ldots ,c_{n}$, el número de filas de ${\textbf {B}}$ debe ser $c_{1},\ldots ,c_{n}$.

Cuando los elementos de las matrices partidas son mostrados en ejemplos numéricos, es conveniente separar las submatrices usando líneas punteadas o separándolas apropiadamente.

Derivadas con vectores y matrices

En los siguientes capítulos va a ser necesario usar fórmulas para encontrar las derivadas parciales de posibilidades y otras funciones escalares de vectores y matrices. Sea $f({\textbf {x}})$ una función continua de los elementos del vector ${\textbf {x}}^{\prime }=[x_{1},\ldots ,x_{p}]$ cuya primera y segunda derivadas

${\frac {\partial f({\textbf {x}})}{\partial x_{i}}}\;\;\;\;{\frac {\partial ^{2}f({\textbf {x}})}{\partial x_{i}\partial x_{j}}}$

existen para todos los puntos \textbf{x} en alguna región del espacio euclideanto de $p$ dimensiones. Se define el operador vectorial de derivada parcial como

${\frac {\partial }{\partial {\textbf {x}}}}={\begin{bmatrix}{\frac {\partial }{\partial x_{i}}}\\\vdots \\{\frac {\partial }{\partial x_{p}}}\end{bmatrix}}$

Aplicando este operador a $f({\textbf {x}})$ se obtiene el vector de derivadas parciales

${\frac {\partial f({\textbf {x}})}{\partial {\textbf {x}}}}={\begin{bmatrix}{\frac {\partial f({\textbf {x}})}{\partial x_{1}}}\\\vdots \\{\frac {f({\textbf {x}})}{\partial x_{p}}}\end{bmatrix}}$

Las siguientes funciones y sus derivadas son especialmente importantes:

1. Si $f({\textbf {x}})$ es constante para toda \textbf{x},

${\frac {\partial f({\textbf {x}})}{\partial {\textbf {x}}}}={\begin{bmatrix}0\\\vdots \\0\end{bmatrix}}$

1. $f({\textbf {x}})={\textbf {a}}^{\prime }{\textbf {x}}$:

${\frac {\partial f({\textbf {x}})}{\partial {\textbf {x}}}}={\begin{bmatrix}a_{1}\\\vdots \\a_{p}\end{bmatrix}}$ Nótese que la forma columnar del vector de derivadas no cambia si expresamos $f({\textbf {x}})={\textbf {x}}^{\prime }{\textbf {a}}$.

1. Para calcular el vector de derivadas de la forma cuadrática ${\textbf {x}}^{\prime }{\textbf {AX}}$ se expresa como

Error al representar (función desconocida «\substack»): {\displaystyle \textbf{x}^\prime\textbf{Ax} = \sum_{i=1}^p \sum_{j=1}^p a_{ij}x_ix_j = a_{ii}x_i^2 + 2x_i \sum_{\substack{j=1 \\ j \neq i}}a_{ij}x_j + \sum\sum_{\substack{h \neq i\\j \neq i}}a_{hj}x_hx_j } y se deriva con respecto a $x_{i}$. Si hacemos que ${\textbf {a}}_{i}^{\prime }$ sea la $i$-ésima fila de la matriz simétrica ${\textbf {A}}$, la derivada parcial de la forma cuadrática con respecto a $x_{i}$ será

$2{\textbf {a}}_{i}^{\prime }{\textbf {x}}=2\sum _{j=1}^{p}a_{ij}x_{j}$

y el vector de derivadas parciales es

${\frac {\partial {\textbf {x}}^{\prime }{\textbf {Ax}}}{\partial {\textbf {x}}}}=2{\textbf {Ax}}$

En particular, el vector de derivadas de la suma de cuadrados ${\textbf {x}}^{\prime }{\textbf {x}}$ es simplemente $2{\textbf {x}}$.

1. En la función cuadrática más general

$h({\textbf {x}})=({\textbf {a}}-{\textbf {Cx}})^{\prime }{\textbf {K}}({\textbf {a}}-{\textbf {Cx}})$

sea \textbf{K} una matriz simétrica de $N\times N$, ${\textbf {C}}=[{\textbf {C}}_{1},\ldots ,{\textbf {C}}_{p}]$ es una matriz de constantes $N\times p$, y ${\textbf {a}}$ es un vector de $N\times 1$. Si hacemos que ${\textbf {u}}={\textbf {a}}-{\textbf {Cx}}$, las derivadas de $h({\textbf {x}})$ se pueden calcular por medio de la regla de la cadena:

${\begin{aligned}{\frac {\partial h({\textbf {x}})}{\partial x_{i}}}&=\sum _{j=1}{\frac {\partial ({\textbf {u}}^{\prime }{\textbf {Ku}})}{\partial u_{j}}}{\frac {u_{j}}{\partial x_{i}}}\\&=-2{\textbf {u}}^{\prime }{\textbf {KC}}_{i}\;\;i=1,\ldots ,p\end{aligned}}$

y

${\frac {\partial h({\textbf {x}})}{\partial {\textbf {x}}}}=-2{\textbf {C}}^{\prime }{\textbf {K}}({\textbf {a}}-{\textbf {Cx)}}$

1. La matriz de derivadas parciales de segundo orden de una función de $p$ variables se conoce como la matriz hesiana :

${\textbf {H}}={\frac {\partial ^{2}f({\textbf {x}})}{\partial {\textbf {x}}^{\prime }\partial {\textbf {x}}}}={\begin{bmatrix}{\frac {\partial ^{2}f({\textbf {x}})}{\partial x_{i}^{2}}}&\cdots &{\frac {\partial ^{2}f({\textbf {x}})}{\partial x_{1}\partial x_{p}}}\\\cdots &\cdots &\cdots \\{\frac {\partial ^{2}f({\textbf {x}})}{\partial x_{1}\partial x_{p}}}&\cdots &{\frac {\partial ^{2}f({\textbf {x}})}{\partial x_{p}^{2}}}\end{bmatrix}}$

Por ejemplo, la matriz hesiana de $f({\textbf {x}})={\textbf {x}}^{\prime }{\textbf {Ax}}$ es simplemente $2{\textbf {A}}$. La matriz hesiana es necesariamente simétrica si $f({\textbf {x}})$ satisface todas las condiciones originales de continuidad y existencia de la primera y segunda derivadas.

Determinación de máximos y mínimos

Una condición necesaria y suficiente para un máximo o mínimo de $f({\textbf {x}})$ en ${\textbf {x}}={\textbf {x}}_{0}$ es que

${\frac {\partial f({\textbf {x}})}{\partial {\textbf {x}}}}={\textbf {0}}$

en ese punto. Tal valor de la función se conoce como un máximo o mínimo estacionario, a diferencia de un extremo global el cual puede existir en la frontera de la región admisible de ${\textbf {x}}$ o como una cúspide o alguna otra forma en donde las derivadas se encuentran indefinidas. Una condición suficiente para un máximo en el punto ${\textbf {x}}_{0}$ donde se satisface la ecuación X es que la matriz hesiana evaluada en ${\textbf {x}}_{0}$ sea negativa definida. Al mismo tiempo, una hesiana positiva definida es una condición suficiente para un mínimo estacionario. Si la hesiana es semidefinida, la prueba no es concluyente, y se deben examinar los términos de más alto orden del desarrollo de $f({\textbf {x}})$ en serie de Taylor en la vecindad del punto estacionario. Si la hesiana es una matriz indefinida de ${\textbf {x}}_{0}$, el punto no es ni un máximo ni un mínimo.

Como una aplicación de éstas reglas, determinemos el extremo estacionario de la función $h({\textbf {x}})$ dada por la ecuación(). Haciendo la ecuación() igual al vector nulo, podemos ver que el punto estacionario debe satisfacer el sistema de ecuaciones lineales

${\textbf {C}}^{\prime }{\textbf {KCx}}={\textbf {C}}^{\prime }{\textbf {KA}}$

y si ${\textbf {C}}^{\prime }{\textbf {KC}}$ es de rango completo $p$, entonces

${\textbf {x}}=({\textbf {C}}^{\prime }{\textbf {KC}})^{-1}{\textbf {C}}^{\prime }{\textbf {Ka}}$

Si ${\textbf {K}}$ es positiva definida y ${\textbf {C}}$ es de rango $p$, la hesiana $2{\textbf {C}}^{\prime }{\textbf {KC}}$ será positiva definida y la solución () minimizará $h({\textbf {x}})$.

Maximización sujeta a restricciones

Con frecuencia será necesario maximizar o minimizar una función $f({\textbf {x}})$ sujeta a la restricción $g({\textbf {x}})=c$ sobre los valores de ${\textbf {x}}$. El método más general y eficiente es el de multiplicadores de Lagrange. Los fundamentos matemáticos de ésta técnica se pueden encontrar en cualquier texto de cálculo moderno. Esencialmente, se forma una función nueva

$h({\textbf {x}},\lambda )=f({\textbf {x}})-\lambda [g({\textbf {x)}}-c]$

Para un valor estacionario restringido se deben cumplir estas condiciones:

${\begin{aligned}{\frac {\partial h({\textbf {x}},\lambda )}{\partial {\textbf {x}}}}&={\frac {\partial f({\textbf {x}})}{\partial {\textbf {x}}}}-\lambda {\frac {\partial g({\textbf {x}})}{\partial {\textbf {x}}}}\\&={\textbf {0}}\\{\frac {\partial h({\textbf {x}},\lambda )}{\partial \lambda }}&=-g({\textbf {x}})+c\\&=0\end{aligned}}$

La segunda condición es simplemente la restricción original. En la práctica uno debe resolver las ecuaciones () con respecto a la variable ${\textbf {x}}$ después de haber eliminado $\lambda$ por medio de manipulación algebraica o por medio de la naturaleza de las ecuaciones.

Como un ejemplo, vamos a encontrar el valor máximo de la función cuadrática $f({\textbf {x}})={\textbf {x}}^{\prime }{\textbf {Ax}}$ donde ${\textbf {A}}$ es una matriz positiva definida, sujeta a la restricción ${\textbf {x}}^{\prime }{\textbf {x}}=1$. Entonces

$h({\textbf {x}},\lambda )={\textbf {x}}^{\prime }{\textbf {Ax}}-\lambda ({\textbf {x}}^{\prime }{\textbf {x}}-1)$

y () es el sistema de ecuaciones lineales

$[{\textbf {A}}-\lambda {\textbf {I}}]{\textbf {x}}={\textbf {0}}$

que no es nada más sino la definición de los valores característicos de ${\textbf {A}}$. Premultiplicando por ${\textbf {x}}^{\prime }$ y usando la restricción nos proporciona $\lambda ={\textbf {x}}^{\prime }{\textbf {Ax}}$. Obviamente, si la forma cuadrática va a ser un máximo, entonces $\lambda$ debe ser la raíz característica más grande de ${\textbf {A}}$, y ${\textbf {x}}$ su vector asociado. La restricción establece que el vector debe ser normalizado a una longitud igual a uno. Similarmente, el valor mínimo de $f({\textbf {x}})$ sujeto a la condición estará dado por el vector característico asociado a la raíz característica más pequeña.

La derivada del determinante de la matriz ${\textbf {A}}$ de $n\times n$ con respecto a su elemento $a_{ij}$ se puede encontrar a partir de la expansión de ${\textbf {A}}$ en los cofactores de la $i$-ésima fila o $j$-ésima columna. Usando la primera expansión

${\begin{aligned}{\frac {\partial |{\textbf {A}}|}{\partial a_{ij}}}&={\frac {\partial }{\partial a_{ij}}}(a_{i1}A_{i1}+\cdots +a_{ij}A_{ij}+\cdots +a_{in}A_{in})\\&=A_{ij}\end{aligned}}$

Si ${\textbf {A}}$ es simétrica,

${\frac {\partial |{\textbf {A}}|}{\partial a_{ii}}}=A_{ii}$

y se puede demostrar de los fórmulas más generales para las derivadas de determinantes, cuyos elementos son a su vez funciones de otras variables, que

${\frac {\partial |{\textbf {A}}|}{\partial a_{ij}}}=A_{ij}$

Más adelante será necesario derivar formas cuadráticas y otros escalares con respecto a los elementos de sus matrices constituyentes. Si ${\textbf {X}}$ es una matriz de $m\times n$ con un elemento general $x_{ij}$, su derivada con respecto a ese elemento es

${\frac {\partial {\textbf {X}}}{\partial x_{ij}}}={\textbf {J}}_{ij}$

donde ${\textbf {J}}_{ij}$ es una matriz de $m\times n$ con uno en la posición $ij$ y ceros en las demás. Si ${\textbf {X}}$ es simétrica

${\frac {\partial {\textbf {X}}}{\partial x_{ij}}}=J_{ij}+J_{ji}\;\;\;i\neq j$

La regla para derivar un producto de matrices es parecida a la de los escalares. Supóngase que los elementos de las matrices conformables ${\textbf {X}}$ y ${\textbf {Y}}$ son funciones $x_{ij}(z),y_{ij}(z)$de alguna variable $z$. Entonces

${\frac {\partial {\textbf {XY}}}{\partial z}}={\frac {\partial {\textbf {X}}}{\partial x}}{\textbf {Y}}+{\textbf {X}}{\frac {\partial {\textbf {Y}}}{\partial z}}$

Esta fórmula lleva a una forma de derivar la inversa de la matriz cuadrada no singular ${\textbf {X}}$. Expresamos

${\textbf {I}}={\textbf {XX}}^{-1}$

Entonces

${\begin{aligned}{\frac {\partial {\textbf {I}}}{\partial x_{ij}}}&={\textbf {J}}_{ij}{\textbf {X}}^{-1}+{\textbf {X}}{\frac {\partial {\textbf {X}}^{-1}}{\partial x_{ij}}}\\&={\textbf {0}}\end{aligned}}$

y

${\frac {\partial {\textbf {X}}^{-1}}{\partial x_{ij}}}=-{\textbf {X}}^{-1}{\textbf {J}}_{ij}{\textbf {X}}^{-1}$

Si ${\textbf {X}}$ es simétrica

${\frac {\partial {\textbf {X}}^{-1}}{\partial x_{ij}}}={\begin{cases}-{\textbf {X}}^{-1}{\textbf {J}}_{ii}{\textbf {X}}^{-1}&i=j\\-{\textbf {X}}^{-1}({\textbf {J}}_{ij}+{\textbf {J}}_{ji}){\textbf {X}}^{-1}&i\neq j\end{cases}}$

Avances

El siguiente capítulo tratará la implementación de rutinas para el manejo de matrices y vectores usando el lenguaje ${\textsf {R}}$. Esto concluye el presente capítulo.

Vectores y Matrices en R

Vectores

Como ya hemos visto, ${\textsf {R}}$ opera sobre objetos y estos objetos pueden ser muchos tipos. Uno de los objetos más sencillos es el tipo vector. Existen otros objetos más complejos tales como matrices, arreglos, dataframes y listas. Nos abocaremos a ellos a su debido tiempo. Para trabajar eficientemente con sus datos, debe poseerse la habilidad de manipular objetos en ${\textsf {R}}$. Es necesario saber como crear, unir, y obtener subconjuntos de estos objetos. Dentro de ${\textsf {R}}$, un objecto tipo vector consiste en un conjunto ordenado de los elementos del mismo tipo.

La longitud de un vector es el número de sus elementos. Como ya se ha dicho, todos los elementos deben ser del mismo tipo, sin embargo existe una excepción a esta regla. Cualquier tipo se puede mezclar con el elemento especial ${\texttt {NA}}$ que quiere decir Not Available (No se encuentra disponible). El tipo ${\texttt {NA}}$ no es un valor en sí lo que implica el uso de funciones especiales. Un vector con dimensión igual a cero será el vector vacío.

Creando vectores

Los vectores pueden ser construídos de varias formas:

> v <- 1 : 10  # Se crea un vector llamado 'v' conteniendo los elementos del 1 al 10.
> v <- c(1, 5, 3) # c() encadena los elementos

Se puede también crear un vector llamando la función:

> v <- vector() # Un vector con longitud cero de tipo logico
> (v <- vector(length=10))
[1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE

Del código anterior se concluye que la función ${\texttt {vector()}}$ produce un vector de tipo ${\texttt {logical}}$ por omisión. Si el vector tiene una longitud, todos los elementos se inicializan con el valor de ${\texttt {FALSE}}$. Para crear un vector con un tipo específico, sólo hay que llamar el tipo con una función, con o sin longitud:

> v <- integer()
> (v <- double(10))
[1] 0 0 0 0 0 0 0 0 0 0
> (v <- character(10))
[1] "" "" "" "" "" "" "" "" "" ""

También puede crear vectores asignándole otros vectores. Por ejemplo, la función ${\texttt {c()}}$ regresa con un vector. Cuando se le asigna un objeto al resultado de ${\texttt {c()}}$, el objeto se vuelve un vector:

> v <- c(0,-10, 1000)
[1]    0  -10 1000

Ya que un vector debe incluír elementos de un tipo sencillo, si se encadenan vectores con ${\texttt {c()}}$, ${\textsf {R}}$ forzará a que todos los tipos se reduzcan al más sencillo. Por ejemplo:

> (x <- c('Amor', 'Sin', 'Barreras'))
[1] "Amor"     "Sin"      "Barreras"
> (y <- c(1,2,3))
[1] 1 2 3
> (z <- c(x,y))
[1] "Amor"     "Sin"      "Barreras" "1"    "2"        ``3

Obsérvese este código con atención. En la línea 1 se asignan tres cadenas de caracteres a ${\texttt {x}}$. En la línea 3 se crea un vector numérico de tipo ${\texttt {double}}$ con los valores 1, 2, y 3. Al encadenar los vectores ${\texttt {x}}$ y ${\texttt {y}}$, ${\textsf {R}}$ convierte todo a modo carácter. Por eso terminan encerrados con comillas. Esto puede parecer confuso, pero al tomar en cuenta que ${\textsf {R}}$ trata de ser coherente, deja de serlo.

Funciones con vectores

Al trabajar con datos es necesario trabajar con vectores. Veremos algunas funciones que permiten manipular vectores como la siguiente:

${\texttt {length()}}$: Para obtener la longitud de un vector

> (t.f <- c(TRUE, TRUE, FALSE, TRUE)); length(t.f)
[1] TRUE TRUE FALSE TRUE
[1] 4

${\texttt {sum()}}$: Para obtener la suma de los elementos de un vector

> (unVec <- 1:20)
[1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20
> (promedio <- sum(unVec)/length(unVec))
[1] 10.5

Aritmética vectorial

Las operaciones binarias (esto es, que requieren de dos operandos) como suma, resta, división y multiplicación se usan en operaciones vectoriales. Las operaciones se realizan elemento por elemento usando los símbolos habituales de +, -, \*, y /.

> (x <- 1.2:6.4)
[1] 1.2 2.2 3.2 4.2 5.2 6.2
> x+1
[1] 2.2 3.2 4.2 5.2 6.2 7.2
> x*2
[1]  2.4  4.4  6.4  8.4 10.4 12.4
> x-0.4
[1] 0.8 1.8 2.8 3.8 4.8 5.8
> x -0.2
[1] 1 2 3 4 5 6
> x / 2
[1] 0.6 1.1 1.6 2.1 2.6 3.1
>

Nótese como en cada línea, los elementos de ${\texttt {x}}$ son alterados individualmente por la operación especificada.

Cuando dos vectores son de la misma longitud y del mismo tipo, el resultado de ${\texttt {x/y}}$ es un vector cuyo primer elemento es ${\texttt {x[1]/y[1]}}$, el segundo es ${\texttt {x[2]/y[2]}}$, y así en adelante:

> (x <- seq(2, 10, by = 2))
[1]  2  4  6  8 10
> (y <- 1:5)
[1] 1 2 3 4 5
> x/y
[1] 2 2 2 2 2

Deseamos hacer notar aquí una aparente inconsistencia: Cuando se divide un vector entre un escalar, obtenemos un vector con cada elemento dividido entre el escalar. Pero cuando se divide un vector entre otro de la misma longitud, el resultado es una división de cada elemento entre su correspondiente contraparte. En realidad, esto es consistente una vez que se explican las reglas: Primero, un valor individual es un vector con longitud de 1. Segundo, si ${\texttt {x>y}}$ y se usa aritmética binaria, ${\textsf {R}}$ iterará sobre los elementos de ${\texttt {y}}$ hasta que su longitud sea igual a la de ${\texttt {y}}$. Sin embargo, ${\textsf {R}}$ se quejará de la operación.

Tomando en cuenta esto, la longitud de cualquier vector es un múltiplo entero de la longitud de un vector de tamaño 1. Por ejemplo, un vector de longitud 10 es 10 veces más grande que un vector de longitud 1. Si se realiza una operación binaria entre dos vectores cuyas longitudes son múltiplos respectivos, ${\textsf {R}}$ procederá sin problemas. Por ejemplo, si se divide un vector ${\texttt {x}}$ de longitud 10 entre un vector ${\texttt {y}}$ de longitud 5, ${\textsf {R}}$ usará los elementos del vector ${\texttt {y}}$ dos veces para realizar la operación. Por otro lado, si se trata de realizar una operación binaria entre un vector de longitud 12 y uno de longitud 7, ${\textsf {R}}$ arrojará un mensaje de aviso quejándose de la operación.

> options(digits=3) # Para no imprimir muchos decimales
> (x <- 1:10); (y <- 1:3); x/y
[1]  1  2  3  4  5  6  7  8  9 10
[1] 1 2 3
[1]  1.0  1.0  1.0  4.0  2.5  2.0  7.0  4.0  3.0 10.0
Warning message:
In x/y : longer object length is not a multiple of shorter object length

Se obtuvo el mensaje de aviso debido a que ${\texttt {x}}$ no es un múltiplo entero de ${\texttt {y}}$.

Esta misma regla se aplica con funciones que implementan aritmética vectorial. Por ejemplo, considérese la función de raíz cuadrada: ${\texttt {sqrt()}}$. ${\textsf {R}}$ sigue la regla de operar sobre un elemento a la vez:

> (x <- 1:5); sqrt(x)
[1] 1 2 3 4 5
[1] 1.00 1.41 1.73 2.00 2.24

Vectores de caracteres

Con frecuencia se necesita manipular caracteres. ${\textsf {R}}$ posee un extenso conjunto de funciones para tratar con caracteres. Se verán en esta sección algunas de estas funciones y se presentarán algunas adicionales en la medida que se necesiten. Las cadenas de caracteres se delimitan por comillas dobles o sencillas. Si se necesitan comillas en una cadena de caracteres, se puede alternar entre comillas dobles y sencillas. Por ejemplo:

> (s <- c("Our friend O'Brien", 'the engineer'))
[1] "Our friend O'Brien" "the engineer"

Si queremos crear una sola cadena a partir de los dos elementos en el vector ${\texttt {s}}$ usamos la función ${\texttt {paste()}}$:

> paste(s[1], s[2])
[1] "Our friend O'Brien the engineer"

Por omisión, ${\texttt {paste()}}$ separa sus argumentos con un espacio. Si se desea un carácter diferente para espaciar los elementos, entonces se usa el argumento ${\texttt {sep}}$:

> paste(s[1], s[2], sep = '-')
[1] "Our friend O'Brien-the engineer"

Si se desea más información sobre la función ${\texttt {paste()}}$ úsese ${\texttt {help(paste)}}$.

Subconjuntos de vectores e índices

Una de las operaciones más frecuentes al manejar vectores es el extraer subconjuntos. Para extraer un subconjunto de un vector, hay que especificar los índices de los elemenos que se desea extraer o excluír. Para demostrar vamos a crear un vector y extraer sus elementos en dos formas distintas:

> x <- 15:30
> x
[1] 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
> c(x[2],x[5]) # Extraemos el segundo y quinto valor del vector x
[1] 16 19
> x[c(2,5)] # Esta es otra forma de hacerlo
[1] 16 19

En la segunda forma, usamos algo llamado vector índice. El vector índice especifica los elementos que se deben extraer. Los vectores índice deben ser de uno de los siguientes cuatro modos: Lógico (TRUE, FALSE), enteros positivos, enteros negativos o cadenas de caracteres (strings). Vamos algunos ejemplos

Ejemplo: Vectores índice lógicos y datos faltantes. He aquí un caso típico. Se adquieren valores numéricos con algunos casos faltantes. El vector de datos que tenemos es de esta forma:

> (x <- c(10,20, NA, 4, NA, 2))
[1] 10 20 NA  4 NA  2

Deseamos calcular el promedio de este vector, de modo que realizamos la operación:

> sum(x)/length(x)
[1] NA

Debido a que ejecutamos operaciones sobre elementos ${\texttt {NA}}$, el resultado va a ser ${\texttt {NA}}$. De modo que necesitamos alguna forma de extraer del vector los datos que no son ${\texttt {NA}}$. Vamos a crear el siguiente vector:

> (i <- is.na(x))
[1] FALSE FALSE  TRUE FALSE  TRUE FALSE

La función ${\texttt {is.na(x)}}$ examina cada elemento del vector ${\texttt {x}}$ y si es ${\texttt {NA}}$ le asigna al correspondiente elemento del vector ${\texttt {i}}$ el valor de ${\texttt {TRUE}}$; de lo contrario, le asigna el valor ${\texttt {FALSE}}$. Ahora vamos a usar al vector ${\texttt {i}}$ como un vector índice para extraer los valores que no son ${\texttt {NA}}$ del vector ${\texttt {x}}$. Tomando en cuenta que si ${\texttt {x}}$ tiene un valor numérico, el elemento correspondiente del vector ${\texttt {i}}$ tiene el valor ${\texttt {FALSE}}$. Si ${\texttt {x}}$ no tiene valor, entonces el elemento correspondiente de ${\texttt {i}}$ es ${\texttt {TRUE}}$. Usamos entonces ${\texttt {i}}$ y la negación de ${\texttt {i}}$ para obtener:

> i ; !i
[1] FALSE FALSE  TRUE FALSE  TRUE FALSE
[1]  TRUE  TRUE FALSE  TRUE FALSE  TRUE
> (y <- x[!i]) # Usamos el vector negado de indices para extraer los valores numericos de x
[1] 10 20  4  2

Ahora calculamos el promedio de ${\texttt {x}}$, con los valores faltantes:

> n <- length(x[!i])
> nuevo.x <- x[!i]
> sum(nuevo.x)/n
[1] 9

Se puede obtener el mismo resultado con

> mean(x, na.rm = TRUE)
[1] 9

La función ${\texttt {mean}}$ es una función intrínseca de ${\textsf {R}}$ que calcula el promedio de un vector. El primer argumento de ${\texttt {mean}}$ es el vector de datos, y el segundo es un argumento lógico que si se le asigna el valor de ${\texttt {TRUE}}$ calcula el promedio sin tomar en cuenta los valores ${\texttt {NA}}$. El valor de ${\texttt {na.rm}}$ por omisión es ${\texttt {FALSE}}$.

Ejemplo: Vector índice de enteros positivos. Esta es otra forma para excluír datos tipo ${\texttt {NA}}$:

> (x <- c(160, NA, 175, NA, 180))
[1] 160  NA 175  NA 180
> (no.na <- c(1,3,5))
[1] 1 3 5
> x[no.na] # Usamos el vector no.na como vector indice
[1] 160 175 180

Vamos ahora a extrar los elementos 20 al 30 de un vector ${\texttt {x}}$ de longitud 100:

> x <- runif(100) # 100 numeros aleatorios entre 0 y 1
> round(x[20:30], 3) # Redondear e imprimir los elementos 20 al 30
[1] 0.340 0.262 0.165 0.322 0.510 0.924 0.511 0.258 0.046 0.418 0.854

La instrucción ${\texttt {20:30}}$ crea un vector de enteros que corren del 20 al 30 y se usa para extraer los correspondientes elementos de ${\texttt {x}}$, mismos que se redondean a tres decimales usando la función ${\texttt {round}}$. Una ligera variación sobre este tema sería:

> (i <- c(20:25, 28, 35:49)) # 20 al 25, 28 y del 35 al 49
[1] 20 21 22 23 24 25 28 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49
> round(x[i],2)
[1] 0.34 0.26 0.17 0.32 0.51 0.92 0.05 0.39 0.69 0.69 0.55 0.43 0.45 0.31 0.58
[16] 0.91 0.14 0.42 0.21 0.43 0.13 0.46

Dejamos que el lector interprete este código como ejercicio.

Ejemplo: Vector índice de enteros negativos. El efecto de un vector índice de enteros negativos es lo opuesto de lo que ocurre con los enteros positivos. Para extraer todos los elementos de ${\texttt {x}}$ ${\textbf {excepto}}$ los elementos del 20 al 30, se expresará como: ${\texttt {x[-(20:30)]}}$. Obsérvese como encerramos a ${\texttt {20:30}}$ entre paréntesis, ya que si no lo hubiéramos hecho, ${\textsf {R}}$ supondría que queremos extraer los elementos -20 al 30 y los índices negativos no son permitidos.

Ejemplo: Vector índice de cadenas de caracteres (strings). Medimos el peso ${\texttt {x}}$ y la altura ${\texttt {y}}$ de 5 personas:

> x <- c(160, NA, 175, NA, 180)
> y <- c(50,65,65,80,70)

Ahora, asociamos nombres con los datos:

> names(x) <- c("Ana", "Beto", "Carlos", "Diana", "Eduardo")

La función ${\texttt {names(x)}}$ les da un nombre a los elementos de ${\texttt {x}}$. Para ver los datos los asociamos por columna con la función ${\texttt {cbind()}}$:

> cbind(x,y)
x  y
Ana     160 50
Beto     NA 65
Carlos  175 65
Diana    NA 80
Eduardo 180 70

Ya que los índices tienen nombre, podemos extraer los elementos usando sus nombres

> x[c("Eduardo","Carlos")]
Eduardo  Carlos
180     175

Arreglos y Matrices

Los arreglos generalizan el concepto de vectores. Como vimos anteriormente, un vector tiene dimensión 1 y una longitud de al menos 0. El i-ésimo elemento de un vector se obtiene a través de la notación de subíndices: ${\texttt {v[i]}}$. Las matrices son arreglos rectangulares en dos dimensiones. Los elementos en las matrices se accede a ellos por medio de la notación ${\texttt {m[i,j]}}$. Los arreglos tienen $k$ dimensiones. Cada elemento de un arreglo es accedido por medio de $k$ índices: ${\texttt {a[i1,i2,\dots ,ik]}}$.

Un objeto de tipo arreglo con dimensión 1 es diferente de un vector porque un arreglo tiene un vector dimensión. El vector dimensión es un vector de enteros positivos. La longitud de este vector proporciona la dimensión del arreglo y el vector dimensión es un atributo del arreglo. El nombre de este atributo es ${\texttt {dim}}$. Observemos aquí algunas instrucciones en ${\textsf {R}}$ que ayudarán a aclarar el concepto:

> (v <- 1:10) # v es un vector
[1]  1  2  3  4  5  6  7  8  9 10
> c('vector?' = is.vector(v), árray?' = is.array(v))
vector?  array?
TRUE   FALSE
> dim(v) <- c(10) # Le asignamos el valor de 10 a dim y ahora es un arreglo
> c('vector?' = is.vector(v), árray?' = is.array(v))
vector?  array?
FALSE    TRUE
> v # el arreglos se imprime igual que los vectores
[1]  1  2  3  4  5  6  7  8  9 10

Un objeto tipo matriz es un arreglo de dos dimensiones. Por lo tanto tiene posee el atributo ${\texttt {dim}}$. Su vector dimensión tiene una longitud igual a 2 donde el primer elemento es el número de filas y el segundo es el número de columnas. He aquí un ejemplo de como crear una matriz con tres columnas y dos filas con ${\texttt {matrix()}}$:

> (m <- matrix(0, ncol=3, nrow=2))
[,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    0    0

Ahora, verificamos si ${\texttt {m}}$ es un arreglo con la función ${\texttt {is.array()}}$:

> c('matrix?' = is.matrix(m), árray?' = is.array(m))
matrix?  array?
TRUE    TRUE

Como nos muestra la función ${\texttt {is.matrix()}}$, ${\texttt {m}}$ es al mismo tiempo una matriz y un objeto de tipo arreglo. Dicho de otra manera, toda matriz es un arreglo, pero no todo arreglo es una matriz. Como cualquier otro objeto, se puede obtener la información sobre los atributos de una matriz con

> attributes(m)
$dim
[1] 2 3

Al igual que los vectores, se pueden crear índices y extraer elementos de los arreglos con vectores índice. En el siguiente ejemplo, crearemos una matriz de 5 columnas y 4 filas y extraeremos una submatriz:

> (m <- matrix(1:20, ncol=5, nrow=4))
[,1] [,2] [,3] [,4] [,5]
[1,]    1    5    9   13   17
[2,]    2    6   10   14   18
[3,]    3    7   11   15   19
[4,]    4    8   12   16   20
> i <- c(2,3); j <- 2:4 # Vectores indice
> m[i,j] # Extraemos las filas 2 y 3, y las columnas 2 a 4.
[,1] [,2] [,3]
[1,]    6   10   14
[2,]    7   11   15

Nótese como la matriz se creó con una secuencia de 20 números. La primera columna se llena primero, luego la segunda, y así en adelante. Esto sucederá siempre como regla general. La matrices se llenan columna a columna ya que el índice de la izquierda es el que corre más rápido. Esta regla se aplica a los arreglos con dimensiones mayores que 2. Se pueden poblar las matrices fila a fila usando el argumento ${\texttt {byrow}}$ en la función ${\texttt {matrix()}}$.

Problemas Resueltos

1. Sea

${\textbf {A}}={\begin{bmatrix}1&2&-1\\-1&3&-1\\2&2&4\end{bmatrix}}\;\;\;\;{\textbf {B}}={\begin{bmatrix}3&2&-1\\2&3&1\\-1&1&3\end{bmatrix}}\;\;\;\;{\textbf {C}}={\begin{bmatrix}2&0\\-1&1\\3&2\end{bmatrix}}$

Realícense las siguientes operaciones que se encuentren definidas:

1. 1. ${\textbf {A+B}}$
   2. ${\textbf {A}}-2{\textbf {B}}$
   3. ${\textbf {A}}^{\prime }+{\textbf {B}}$
   4. ${\textbf {A}}+{\textbf {C}}$
   5. (${\textbf {A+B)}}^{\prime }$
   6. ($3{\textbf {A}}^{\prime }-2{\textbf {B}})^{\prime }$

> (A <- matrix(c(1,2,-1,-1,3,-1,2,2,4), ncol=3, nrow=3, byrow=TRUE))
     [,1] [,2] [,3]
[1,]    1    2   -1
[2,]   -1    3   -1
[3,]    2    2    4
> (B <- matrix(c(3,2,-1,2,3,1,-1,1,3), ncol=3, nrow=3, byrow=TRUE))
 [,1] [,2] [,3]
[1,]    3    2   -1
[2,]    2    3    1
[3,]   -1    1    3
> (C <- matrix(c(2,0,-1,1,3,2), ncol=2, nrow=3, byrow=TRUE))
 [,1] [,2]
[1,]    2    0
[2,]   -1    1
[3,]    3    2
> A + B
[,1] [,2] [,3]
 [1,]    4    4   -2
 [2,]    1    6    0
 [3,]    1    3    7
> A - 2*B
[,1] [,2] [,3]
 [1,]   -5   -2    1
 [2,]   -5   -3   -3
 [3,]    4    0   -2
> t(A)+B # t(A) es la funcion para obtener la transpuesta de una matriz
[,1] [,2] [,3]
 [1,]    4    1    1
 [2,]    4    6    3
 [3,]   -2    0    7
> A + C
Error in A + C : non-conformable arrays
> t(A+B)
[,1] [,2] [,3]
 [1,]    4    1    1
 [2,]    4    6    3
 [3,]   -2    0    7
> t(3*t(A) - 2*B)
[,1] [,2] [,3]
 [1,]   -3    2   -1
 [2,]   -7    3   -5
 [3,]    8    4    6

1. Si

${\textbf {P}}={\begin{bmatrix}3&2&1\\2&5&-1\\1&-1&3\end{bmatrix}}\;\;\;\;{\textbf {Q}}={\begin{bmatrix}1&2&2&1\\1&1&1&4\\1&2&2&1\end{bmatrix}}\;\;\;\;{\textbf {R}}={\begin{bmatrix}1&2\\-5&2\\3&-1\\-2&2\end{bmatrix}}$
${\textbf {x}}={\begin{bmatrix}1\\0\\1\end{bmatrix}}\;\;\;\;{\textbf {y}}={\begin{bmatrix}2\\3\\2\end{bmatrix}}\;\;\;\;{\textbf {z}}={\begin{bmatrix}-1\\-2\\-3\\-4\end{bmatrix}}$

Calcúlense aquellos productos matriciales que se encuentren definidos:

1. 1. ${\textbf {PQ}}$
   2. ${\textbf {PQR}}$
   3. ${\textbf {QR}}^{\prime }$
   4. ${\textbf {yx}}^{\prime }$
   5. ${\textbf {x}}^{\prime }{\textbf {y}}$
   6. ${\textbf {x}}^{\prime }{\textbf {Px}}$
   7. ${\textbf {x}}^{\prime }{\textbf {Py}}$
   8. ${\textbf {P}}({\textbf {x+y}})$

> (P <- matrix(c(3,2,1,2,5,-1,1,-1,3), nrow=3, ncol=3, byrow=TRUE))
     [,1] [,2] [,3]
[1,]    3    2    1
[2,]    2    5   -1
[3,]    1   -1    3
> (Q <- matrix(c(1,2,2,1,1,1,1,4,1,2,2,1), nrow=3, ncol=4, byrow=TRUE))
     [,1] [,2] [,3] [,4]
[1,]    1    2    2    1
[2,]    1    1    1    4
[3,]    1    2    2    1
> (R <- matrix(c(1,2,-5,2,3,-1,-2,2), nrow=4, ncol=2, byrow=TRUE))
     [,1] [,2]
[1,]    1    2
[2,]   -5    2
[3,]    3   -1
[4,]   -2    2
> (x <- c(1,0,-1))
[1]  1  0 -1
> (y <- c(2,3,2))
[1] 2 3 2
> (z <- c(-1,-2,-3,-4))
[1] -1 -2 -3 -4
> P %*% Q # En R, el simbolo de multiplicacion de matrices es %*%
     [,1] [,2] [,3] [,4]
[1,]    6   10   10   12
[2,]    6    7    7   21
[3,]    3    7    7    0
> P %*% Q %*% R
     [,1] [,2]
[1,]  -38   46
[2,]  -50   61
[3,]  -11   13
> Q %*% t(R)
Error in Q %*% t(R) : non-conformable arguments
> y %*% t(x)
     [,1] [,2] [,3]
[1,]    2    0   -2
[2,]    3    0   -3
[3,]    2    0   -2
> t(x) %*% y
     [,1]
[1,]    0
> t(x) %*% P %*% x
     [,1]
[1,]    4
> t(x) %*% P %*% y
     [,1]
[1,]    9
> P %*% (x+y)
     [,1]
[1,]   16
[2,]   20
[3,]    3