Ver el código de «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 <math>E_{1}, ..., E_{k}</math>. 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 <math>p_{i}</math> entre cero y uno. Este valor numérico se le conoce como la ''probabilidad'' <math>P(E_{i})</math> de ese evento. Entonces <math>p_{i}</math> es la probabilidad de que ocurra el evento <math>E_{i}</math> 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 
 <math>\begin{matrix}P(E_{i}\cap E_{j})&=&  0\\ P(E_{i}\cup E_{j})&=& p_{i}+p_{j}\end{matrix}</math> donde el símbolo de intersección <math>\cap</math> denota el evento <math>E_{i}</math>
y <math>E_{j}</math> y el símbolo de unión <math>\cup</math> indica el evento <math>E_{i}</math> o <math>E_{j}</math>. 

Usando la propiedad aditiva de las probabilidades de eventos mutuamente exclusivos la probabilidad total del conjunto de eventos es
 <math>P(E_{1}\cup E_{2}\cdots\cup E_{k})=p_{1}+\cdots+p_{k}=1</math>
Ahora se le asigna el valor numérico <math>x_{i}</math> al i-ésimo resultado, donde, por conveniencia, los resultados se ordenan ascendentemente según sus valores <math>x_{i}</math>. Se define entonces la variable aleatoria discreta <math>X</math> como la cantidad que toma el valor <math>x_{i}</math> con probabilidad de <math>p_{i}</math> en cada corrida del experimento. Por ejemplo, si el experimento consiste en lanzamientos de una moneda, el resultado "cara" puede tener una cierta probabilidad <math>p_{1}</math> y el resultado "cruz" puede tener una probabilidad <math>p_{2}</math>. Esta variable aleatoria puede ser descrita entonces por una ''función de probabilidad'' <math>f(x_{i})</math> la cual especifica las probabilidades con las que la variable aleatoria <math>X</math> adquiere los valores de cara o cruz.

{|border="1" cellspacing="2" cellpadding="2" align="center"
|<math>X</math>  
| <math>f(x_{i})</math>
|-
|cara  
| <math>p_{1}</math>
|-
|cruz  
| <math>p_{2}</math>
|-
|}

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 <math>X</math> definida sobre el dominio de los números reales se caracteriza por su ''función de distribución''
 <math>F(x)=P(X\leq x) \;\; -\infty<x<\infty </math>

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

Si <math>F(x)</math> es una función absolutamente continua, la versión continua
de la función discreta de probabilidad es la ''función de densidad''
 <math>f(x)=\dfrac{dF(x)}{dx}</math>

Al mismo tiempo, por la propiedad de continuidad absoluta
 <math>F(x)=\int_{-\infty}^{x}f(u)\; du </math>
y de esta definición se genera el término equivalente ''función acumulativa de distribución'' de <math>F(x)</math>. 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 <math>f(x)</math> es igual a cero, y <math>F(x)</math> 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 <math>f(x)</math>, los valores en el intervalo se definen como la ''población'' o el ''universo'' de la variable aleatoria <math>X</math>.

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 <math>\frac{1}{2}(a+b)</math> de la densidad rectangular y <math>\mu</math> 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 <math>b-a</math> de la variable rectangular, el parámetro exponencial <math>\beta</math>, y el parámetro <math>\sigma</math> 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 <math>f(x;\alpha,\beta)</math> se puede expresar como
 <math>f\left(\dfrac{x-\alpha}{\beta}\right)</math>
entonces <math>\alpha</math> es un parámetro de ubicación y <math>\beta</math> es un parámetro de escala. Se hace notar que todas las funciones de densidad de esta forma contienen el factor <math>1/\beta</math> asociado con el elemento diferencial <math>dx</math>.

Las variables aleatorias y sus densidades tienen una interpretación física: considérese la función de densidad <math>f(x)</math> de la variable aleatoria <math>X</math> como la función que define la densidad de una barra sólida que se encuentra posicionada en el eje de las <math>x</math>. El k-ésimo momento alrededor del origen es
 <math>\mu_{k}^{'}=\int_{-\infty}^{\infty}x^{k}f(x)\; dx </math>
El primer momento es la ''media'' o el ''valor esperado'' 
 <math>E(X)=\int_{-\infty}^{\infty}xf(x)\; dx</math>
de la variable aleatoria y corresponde al concepto físico del centro de gravedad de la vara. El símbolo <math>E</math> denota la operación de calcular el valor esperado, y se conoce como el ''operador de esperanza''. Si <math>c</math> y <math>k</math> son cantidades no aleatorias, entonces podemos aseverar las siguientes propiedades de las esperanzas:
 <math> \begin{align}
E(c) & = c\\
E(cX) & = cE(X)\\
E(k+cX) & = k+cE(X)
\end{align}
</math>

Como ejemplo, considérese la función de densidad rectangular de la figura [[Archivo:Distribuciones Rectangular y Exponencial.jpg|thumb|Figura 1 - Funciones de Densidad Rectangular y Exponencial]]. Por definición:
 <math>\begin{align}
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{align}
</math>

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 <math>w=b-a</math> como el parámetro de escala. Entonces, la nueva densidad será:
 <math>
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}
</math>

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 <math>X</math> como
 <math>
\text{var}(X)
</math>
o también como
<math>
\sigma_{x}^{2}
</math>

Por definición
 <math>\begin{align}
\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{align} </math>

y el simbolo '''var''' aplicado a cualquier variable aleatoria denotará el cálculo de la varianza. Si <math>c</math> y <math>k</math> son constantes, la varianza posee las siguientes propiedades:
 <math>
\begin{align}
\text{var}(c)&= 0\\
\text{var}(cX)&= c^{2}\text{var}(X)\\
\text{var}(X+c)&= \text{var}(X)
\end{align}
</math>

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 <math>c</math> altera la varianza en <math>c^{2}</math>. 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 <math>X</math>. Se puede demostrar que las varianzas de las funciones de densidad rectangular, exponencial, y normal [[Archivo:TresFD.jpg|thumb|Figura 2 - Funciones de Densidad Rectangular, Exponencial, y Gaussiana]] son <math>\frac{1}{12}(b-a)^{2}</math>, <math>\beta^{2}</math>, y <math>\sigma^{2}</math> 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 <math>X</math> es la raiz cuadrada positiva de la varianza:

 <math>\sigma_{x}=+\sqrt{\text{var}(X)}</math>

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 <math>X_{1},\ldots,X_{p}</math> como

 <math>
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}
</math>
donde <math>f(u_{1},\ldots,u_{p})</math> es la densidad conjunta, entonces las variables se dice que son independientes ''si y solo si''
 <math>F(x_{1},\ldots,x_{p})=F_{1}(x_{1})\cdots F_{p}(x_{p})
</math>
donde <math>F_{i}(x_{i})</math> es la función de distribución de la variable sencilla <math>X_{i}</math>. Alternativamente, la independencia es válida si y sólo si la factorización de
 <math>f(x_{1},\ldots,x_{p})=f_{1}(x_{1})\cdots f_{p}(x_{p})
</math>
de la densidad conjunta en el producto de las funciones individuales de densidad es válida.

El momento de productos
 <math>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}
</math>
de las variables <math>X_{1},\ldots,X_{p}</math> se factoriza como el producto
 <math>
E(X_{1}^{k_1})\cdots E(X_{p}^{k_{p}})
</math>
de los <math>k</math> 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:
<math>
\text{var}(X_{1}+\cdots+X_{p})=\text{var}(X_{1})+\cdots+\text{var}(X_{p}).
</math>

== 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
 <math>
f(x)=\dfrac{1}{\sqrt{2\pi\sigma}}\exp\left[-\dfrac{1}{2}\left(\dfrac{x-\mu}{\sigma}\right)^{2}\right]
\;\; -\infty<x<\infty
</math>
y que su función de distribución se encuentra dada por la integral
 <math>
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)
</math>

donde <math>\Phi(z)</math> 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: <math>N(\mu,\sigma^{2})</math>;
la distribución estandar será indicada por <math>N(0,1)</math>.

Los valores de <math>\Phi(z)</math> 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 <math>z_{\alpha}</math> tal que
 <math>

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

</math>

Estas relaciones entre porcentajes y sus probabilidades se ilustran en la figura 3 [[Archivo:Normal dist density cumulative 0001.jpg|thumb|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 <math>X_{1},\ldots,X_{N}</math> se encuentran distribuídas independientemente de acuerdo a alguna función común de distribución con media <math>E(X_{i})=\mu</math> y varianza finita <math>\text{var}(X_{i})=\sigma^{2}</math>, entonces a medida que el número de variantes <math>N</math> aumenta sin cotas, la variante
 <math>
Z=\dfrac{\sum_{i=1}^{N}X_{i}-N_{\mu}}{\sqrt{N\sigma^{2}}}
</math>
: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 <math>\sqrt{N}(\bar{X}-\mu)</math> de la media de las muestras de una secuencia de variables aleatorias independientes
 <math>
\bar{X}=\frac{1}{N}(X_{i}+\ldots+X_{N})
</math>
cuya distribución común tiene media <math>\mu</math> y varianza <math>\sigma^{2}</math> tiende a distribuírse como una variante de tipo <math>N(0,\sigma^{2})</math> a medida que <math>N</math> se vuelve más grande. Una y otra vez estaremos
invocando una propiedad muy importante de la distribución normal: Si <math>X_{1},\ldots,X_{n}</math> son variantes independientes con distribuciones <math>N(\mu_{1},\sigma_{1}^{2}),\ldots,N(\mu_{n},\sigma_{n}^{2})</math> entonces la combinación lineal
 <math>Y=a_{1}X_{1}+\ldots+a_{n}X_{n} </math>
también se encuentra normalmente distribuída con media <math>\sum_{i=1}^{n}a_{i}\mu_{i}</math> y varianza <math>\sum_{i=1}^{n}a_{i}^{2}\sigma_{i}^{2}</math>. 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 <math>\left( \Gamma(x) \right) </math>. 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 <math>X_{1},\ldots,X_{n}</math> son independientes y se encuentran normalmente distribuídas con media cero y varianza igual a uno, entonces 
 <math>
  \chi^{2}=X_{1}^{2}+\ldots+X_{n}^{2}
</math>
tiene la función de densidad
 <math>
  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
</math>
y se dice que es una variante ''ji-cuadrada'' con <math>n</math> 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 <math>n</math>. Aún más, si las <math>X_{i}</math> se encuentran independientemente distribuídas como variantes <math>N(\mu_{i},\sigma_{i}^{2})</math>, entonces la cantidad
 <math>\sum_{i=1}^{n}\dfrac{(X_{i}-\mu_{i})^{2}}{\sigma_{i}^{2}}
</math>
tiene una distribucion ji-cuadrada con <math>n</math> 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 <math>\textsf{R}</math> 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 <math>n</math> grados de libertad como:
 <math>
\chi_{\alpha;n}^{2}
</math>
donde <math>\alpha=P(\chi^{2}>\chi_{\alpha;n}^2)</math>.

'''La distribución t'''

La variable aleatoria <math>t</math> con <math>n</math> grados de libertad se define como el cociente 
 <math>
t=\dfrac{z}{\sqrt{\chi^{2}/n}}
</math>

de una variante normal estándar <math>z</math> y la raíz cuadrada de una variante ji cuadrada independiente dividida esta por sus <math>n</math> grados de libertad. La variable <math>t</math> es adimensional y su función de densidad
 <math>
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
</math>
depende solamente del parámetro de grados de libertad <math>n</math>. Veremos más adelante como generar valores de una distribución <math>t</math> con <math>n</math> grados de libertad usando el lenguaje <math>\textsf{R}</math>.

Se expresara el porcentaje 100α como <math>t_{\alpha;n}</math> o
 <math>
 \alpha = P(t > t_{\alpha;n})
</math>

'''La distribución F'''

El cociente
 <math>
 F = \dfrac{\chi_1^2/m}{\chi_2^2/n}
</math>
de variantes ji-cuadradas independientes divididas por sus respectivos grados de libertad, tiene la distribución <math>F</math> o ''cociente de varianzas'' con función de densidad
 <math>
 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
</math>
Se denotará el porcentaje superior 100α de la distribución <math>F</math> con <math>m,n</math> grados de libertad por medio de <math>F_{\alpha; m,n}</math>:
 <math>
 \alpha =  P( F > F_{\alpha; m,n})
</math>
De la definición de la variante <math>F</math> se deduce que los puntos porcentuales inferiores se pueden obtener de los recíprocos de los valores superiores con los grados de libertad invertidos:
 <math>
 F_{1-\alpha;m,n} = \dfrac{1}{ F_{\alpha; m,n}}
</math>

 '''Una Nota Breve Sobre la Funcion <math>\Gamma(n)</math>:''' La función <math>\Gamma(n)</math> 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:
 <math>
\Gamma(n) = (n-1)!

</math>
 La función gamma se puede definir como una integral definida en <math>R[z]>0</math> (La forma integral de Euler):
 <math>
\Gamma(z) \equiv \int_0^\infty t^{z-1}e^{-t}\;dt</math>
 o también
 <math>
\Gamma(z) = \int_0^\infty e^{-t^2}t^{2z-1}\; dt.</math>

== 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 <math>F(x)</math> 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 <ref name="Her">[Heurística: Técnicas empíricas usadas para resolver problemas, para el aprendizaje y el
descubrimiento. Cuando no es posible realizar una búsqueda exhaustiva, se usan métodos heurísticos para encontrar
una solucion satisfactoria.]</ref> el cual lleva a obtener estimaciones con algunas propiedades requeridas.

Empecemos suponiendo que los valores de la variable aleatoria continua <math>X</math> 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
<math>f(x; \theta_1, \ldots, \theta_k)</math> de <math>X</math> se sabe a partir de algunas
suposiciones anteriores aunque no se sepan los valores de los parámetros <math>\theta_i</math>.

Definamos ahora el ''espacio de parámetros'' de una densidad. Supongamos
que la densidad depende de un sólo parámetro real <math>\theta</math>, como en el caso de
una exponencial descendente. Entonces el espacio de parámetros de <math>\theta</math> es
la porción del eje de los reales que contiene todos los valores permitidos de
<math>\theta</math>. 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
<math>k</math> parámetros sería alguna región del espacio euclideano en <math>k</math>
dimensiones. Por ejemplo, el de una distribución normal sería la mitad
superior del plano cartesiano con ejes <math>(\mu, \sigma^2)</math> (Fig.4).
[[Archivo:EspacioDeParametros.jpg|thumb|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 <math>X</math>. 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 <math>X</math>. 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 <math>N</math> unidades aleatoriamente. Sus
valores observados de <math>X</math> se denotaran como <math>x_1, \ldots, x_N</math>. El primer
problema es estimar los parametros <math>\theta_i</math> a traves de una función
apropiada de las observaciones. Tales estimaciones serán denotadas por un
parámetro con capucha: <math>\hat{\theta_i}</math>. Estas cantidades serán llamadas
''estimaciones puntuales'', ya que son valores asignados de
<math>\theta</math>, 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. <math>E(\hat{\theta})=\theta</math>.

* ''Consistencia''. A medida que aumenta el tamaño de la muestra, <math>\hat{\theta}</math> debe converger a <math>\theta</math> 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 <math>\hat{\theta}</math> es suficiente si contiene toda la información en las observaciones para estimar <math>\theta</math>. Esto es, el conocimiento directo de los valores <math>x_1, \ldots, x_N</math> no proporcionará más información que la que contiene <math>\hat{\theta}</math>. En términos matemáticos, la distribución conjunta de las <math>x_i</math> para un valor fijo de <math>\hat{\theta}</math> no depende de <math>\theta</math>.

'''Algunas Observaciones Sobre el Sesgo:''' <pre>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.
</pre>
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 <math>x_i, \ldots, x_N</math>:
 <math>
 L(\theta_1, \ldots, \theta_k) = \prod_{i=1}^{N} f(x_i; \theta_1, \ldots,
\theta_k)
</math>

La función de posibilidad es una medida relativa de la posibilidad de la
muestra en particular <math>x_i, \ldots, x_N</math>. El método de máxima posibilidad
demanda que los estimadores de <math>\theta_i</math> 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:
 <math>
  l(\theta_i, \ldots, \theta_k) = \ln L(\theta_1, \ldots, \theta_k)
</math>
Entonces, si la función de posibilidad tiene un máximo relativo, las
estimaciones de <math>\theta_i</math> se pueden encontrar resolviendo el sistema de <math>k</math>
ecuaciones simultáneas
 <math>
  \dfrac{\partial l(\theta_1, \ldots, \theta_k)}{\partial \theta_j} = 0; \ \ j=1,
\ldots,k
 </math>
para encontrar el valor de los estimadores
<math>  \hat{\theta_1}, \ldots, \hat{\theta_k}
</math>
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 <math>a</math> y <math>b</math> para
la densidad rectangular. La función de posibilidad de la muestra <math>x_1, \ldots, x_N</math>, es
 <math>
L(a,b) = \dfrac{1}{(b-a)^N}
</math>
Es evidente que <math>a</math> no puede ser menor que la observación más pequeña, y que
<math>b</math> no puede ser mayor que la observación más grande. Si se denotan las
observaciones ordenadas como <math>x_{(1)}\leq\cdots\leq x_{(N)}</math>,
 <math>
a \leq x_{(1)}\leq\cdots\leq x_{(N)} \leq b
</math>
La función de posibilidad adquiere su valor máximo cuando <math>(b-a)</math> adquiere su valor m+inimo lo que es consistente con
el segundo grupo de desigualdades, y las estimaciones que minimizan el intervalo son
 <math>
\hat{a}=x_{(1)} \; \;\hat{b} = x_{(N)}
</math>

'''Ejemplo:'''
La función de posibilidad de la muestra <math>x_1, \ldots, x_N</math> de <math>N</math> variables
aleatorias normales independientes es
 <math>
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]
</math>
tomando el logaritmo
 <math>
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
</math>
Las derivadas parciales con respecto a <math>\mu</math> y <math>\sigma^2</math> son
 {|
  | <math>\dfrac{\partial l(\mu,\sigma^2)}{\partial \mu} = \dfrac{1}{\sigma^2}\sum_{i=1}^{N}(x_i-\mu)</math>
  |-
  | <math>\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</math>
 |}

Igualando a cero y cancelando factores, queda el siguiente sistema de
ecuaciones simultáneas
 {|
  | <math>  \sum_{i=1}^N x_i -N\mu =  0 </math>
  |-
  |<math>\sum_{i=1}^N (x_i -\mu)^2 - N\sigma^2 =  0</math>
 |}

Despejar el estimador de <math>\mu</math> primero, y luego usar ese valor para obtener la
solución del estimador de <math>\sigma^2</math>:
{|
 | <math>  \hat{\mu} = \dfrac{1}{N}\sum_{i=1}^N x_i  =  \bar{x}</math>
 |-
 |<math>\hat{\sigma}^2 = \dfrac{1}{N}\sum_{i=1}^N (x_i - \bar{x})^2</math>
 |-
 |<math> =  \dfrac{1}{N}\left[ \dfrac{1}{N}\sum_{i=1}^N x_i^2 - \dfrac{1}{N}(\sum
x_i)^2 \right]</math>
|}
Estas estimaciones son intuitivamente plausibles, ya que la media de la
población se estima a través de la media de la muestra <math>\bar{x}</math>, y la
estimación de la varianza es el promedio de la desviación cuadrada de la
media.

Determinemos ahora si <math>\hat{\mu}</math> y <math>\hat{\sigma}^2</math> son estimadores
insesgados. Reemplacemos las observaciones por las variables aleatorias <math>X</math> y
tomemos las esperanzas:
{|
 |<math> E(\hat{\mu}) </math>
 |<math>= \dfrac{1}{N}\sum_{i=1}^N E(X_i) </math>
 |<math>= \dfrac{1}{N}\sum_{i=1}^N \mu = \mu </math>
|}
Esto significa que la media muestral es un estimador insesgado de <math>\mu</math>. El
calcular la esperanza de <math>\hat{\sigma}^2</math> involucra un poco más de cómputo:
{|
  |<math>E(\hat{\sigma}^2)</math>
  |-
  |<math>= \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]</math>
  |-
  |<math>=  \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]</math>
  |-
  |<math>=  \dfrac{1}{N}\left[NE(X_i^2) - E(X_i^2) - \dfrac{1}{N}\sum\sum E(X_iX_j)\right] </math>
  |-
  |<math>=  \dfrac{1}{N} \left\lbrace (N-1) E(X_i^2) - \dfrac{1}{N} N(N-1) [E(X_i)]^2 \right\rbrace </math>
  |-
  |<math> = \dfrac{N-1}{N} \sigma^2 </math>
|}
Entonces <math>\hat{\sigma}^2</math> no es un estimador insesgado de la varianza. Sin
embargo, si reemplazamos el divisor <math>N</math> en la fórmula original por <math>N-1</math>, se
eliminará el sesgo y se obtendrá la conocida formula de la varianza de la muestra
 <math>
s^2 = \dfrac{1}{N-1} \left[\sum_{i=1}^N x_i^2 - \dfrac{1}{N}\left( \sum x_i \right) ^2 \right].
</math>

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:
{| cellpadding="5" cellspacing="0" align="center"
 | '''Hipótesis original:'''
 | <math>\mu=2.43</math>
 | <math>\sigma^2=0.04</math>
 |-
 | '''Hipotesis alternativa:'''
 | <math>\mu>2.43</math>
 | <math>\sigma^2<0.04</math>
 |-
|}'
Le llamaremos a la descripción original de la población como la '''hipótesis nula''; ésta se denota normalmente como <math>H_0</math>. La
'''hipótesis alternativa''' será denotada como <math>H_1</math>. Se supondrá también que la población se encuentra normalmente distribuída, pero por lo
general se asume por omisión. <math>H_0</math> 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, <math>H_1</math> 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 <math>X</math> se
encuentra normalmente distribuída con varianza desconocida, las hipótesis

 <math>H_0: \; \mu = \mu_0</math>
y
 <math>H_1: \; \mu = \mu_1 > \mu_0</math>

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 <math>\mu</math> en <math>\mu_0</math> y en <math>\mu_1</math>.
Habría entonces que desarrollar pruebas estadísticas para tal hipótesis de
tal forma que no se afectaran por el valor desconocido de <math>\sigma^2</math>.
[[Archivo:Fig1 4.jpg|thumb|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 <math>X</math> una variable aleatoria continua con los números reales como su dominio, y su función de densidad
depende de los parámetros <math>\theta_1,\ldots,\theta_k</math>. El ''espacio muestral'' <math>W</math> de todos los resultados posibles
de <math>N</math> observaciones en <math>X</math> ser el espacio N-dimensional euclideano. Una muestra particular de observaciones se
expresará en forma de vector como <math>\left[x_1, \ldots, x_N\right]</math> y denotará un punto en el espacio muestral. Sean
<math>\omega_1</math> y <math>\omega_2</math> dos regiones disjuntas en el espacio de parámetros <math>\Omega</math>, y se proponen las siguientes
hipótesis sobre los parámetros de la distribución de <math>X</math>.

 <math>
  H_0:  [\theta_1, \ldots, \theta_k]\quad \text{se encuentra contenida en }
\omega_1 </math>
 <math>
  H_1: [\theta_1, \ldots, \theta_k]\quad \text{se encuentra contenida en }
\omega_2
</math>

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 <math>H_0</math> sobre <math>H_1</math> y otros
conjuntos apoyarían la alternativa, la regla de decisión debe tener la
siguiente forma:
 <math>
 \text{Aceptar } H_0\text{ si } \left[x_1, \ldots, x_N\right] \text{cae en } W-w
</math>
 <math> \text{Aceptar } H_1\text{ si } \left[x_1, \ldots, x_N\right] \text{cae en } w 
 </math>
donde <math>w</math> es una parte en particular del espacio muestral llamada '''región crítica''' o '''región de rechazo''' de <math>H_0</math>. Si el verdadero estado de la
naturaleza se encuentra descrito por <math>H_0</math> o <math>H_1</math>, 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 <math>H_1</math> como
verdadera cuando en realidad <math>H_0</math> es la correcta. Un error del segundo tipo, o '''error Tipo II''' se comete cuando se toma <math>H_0</math> como verdadera cuando es
<math>H_1</math> la que describe el verdadero estado de la naturaleza. Resumimos los
resultados de las acciones anteriores en la siguiente tabla:

{| border="1" cellpadding="5" cellspacing="0" align="center"
|
!colspan="2"|Estado de la naturaleza
|-
|'''Decisión'''
|<math>H_0</math> verdadera 
|<math>H_1</math> verdadera
|-
|Se acepta <math>H_0</math>
|Correcto
|Error Tipo II
|-
| Se acepta <math>H_1</math>
|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 <math>\alpha</math> de un error
de Tipo I es la probabilidad de que las observaciones muestrales se encuentren
en la región crítica <math>w</math> cuando <math>H_0</math> es verdadera, o lo mismo que
 <math>
  \alpha = P( \left[x_1, \ldots, x_N\right] \in w|H_0 \text{ es verdadera})
</math>
donde el símbolo <math>\in</math> denota la inclusión en el conjunto <math>w</math> y la barra
vertical indica que el enunciado de probabilidad esta condicionado a la certeza
de la hipótesis nula. El símbolo <math>\alpha</math> se conoce como el '''tamaño'''
de la prueba. La probabilidad de un error del Tipo II es
 <math>
  \beta = P( \left[x_1, \ldots, x_N\right] \in W-w|H_1 \text{ es verdadera})
</math>
El complemento <math>1-\beta</math> 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 <math>\alpha</math> también aumentarán <math>\beta</math>. Al mismo tiempo,
un aumento de <math>\alpha</math> resultará en una reducción de <math>\beta</math>. El enfoque
clásico de las pruebas de hipótesis requiere una prueba con <math>\alpha</math> fijo
cuya región de rechazo se elige de tal manera que minimize <math>\beta</math>, o que
maximize la potencia. Si ambas hipótesis son sencillas, esto es, tienen la
forma
{| cellpadding="5" cellspacing="0" align="center"
| 
| <math>\theta_1 = \theta_{10}</math>
|
|
|<math>\theta_1 = \theta_{11}</math>
|-
|<math>H_0:</math>
|<math>\cdots</math>
|
|<math>H_1:</math>
|<math>\cdots</math>
|-
|
|<math>\theta_k = \theta_{k0}</math>
|
|
|<math>\theta_k = \theta_{k1}</math>
|}

entonces el lema de Neyman-Pearson afirma que la prueba más poderosa con
tamaño <math>\alpha</math> tendrá una región crítica definida por la siguiente regla
de decisión:
 <math>  \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
</math>
y
 <math>
  \text{Aceptar } H_1 \text{ si }\  \lambda < c
</math>
donde <math>c</math> es una constante tal que
 <math>  P(\lambda < c) = \alpha
</math>
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 <math>\lambda</math> es menor que <math>c</math>.

Nótese que el lema de Neyman-Pearson requiere que ambas hipótesis sean
simples, esto es, que las regiones <math>\omega_1</math> y <math>\omega_2</math> 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''':
 <math>  \lambda = \dfrac{L(\hat{\omega})}{L(\hat{\Omega})}
</math>
donde
 <math>  L(\hat{\Omega}) = f(x_1; \hat{\theta}_{10}, \ldots,\hat{\theta}_{k0}) \cdots
f(x_N;\hat{\theta}_{10}, \ldots,
\hat{\theta}_{k0})
</math>
es la función de posibilidad maximizada bajo la suposición de que
 <math>
  H_0: \left[\theta_1, \ldots, \theta_k \right] \in \omega
</math>
es cierta, y que <math>L(\hat{\Omega})</math> es la posibilidad maximizada con <math>\theta_1,
\ldots, \theta_k</math> la cual se permite tomar valores a través del todo el
espacio de parámetros <math>\Omega</math>. Aceptamos <math>H_0</math> si <math>\lambda > c</math> y si no, se
acepta la alternativa
 <math>  H_1: \left[\theta_1, \ldots, \theta_k \right] \in \Omega-\omega
</math>
La constante <math>c</math> se escoge de tal forma que <math>P(\lambda < c |H_0 \text{ es verdadera})
\le \alpha</math>. Tal como en el caso de las hipótesis compuestas, cuando el
verdadero tamaño de la prueba es realmente menor o igual a <math>\alpha</math>, entonces
decimos que la prueba es de nivel <math>\alpha</math>. Se puede demostrar que si <math>H_0</math> es
verdadera, la variante
 <math>  \chi^2 = -2 \ln \lambda
</math>
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 <math>\Omega</math> y el
subespacio de la hipótesis nula <math>\omega</math>, o lo que es mismo, el número de
parámetros determinados por <math>H_0</math>.

=== Pruebas de la Media de una Variante Normal con Varianza Conocida ===
Sea <math>x_1, \ldots, x_N</math> una muestra de observaciones independientes de la
variable aleatoria con distribución <math>N(\mu, \sigma^2)</math>. La varianza
<math>\sigma^2</math> es conocida aunque <math>\mu</math> no. Basándonos sobre las observaciones
muestrales deseamos probar la hipótesis
 <math> H_0: \mu = \mu_0
</math>
de que la media poblacional tiene un cierto valor determinado <math>\mu_0</math> contra la
hipótesis alternativa:
 <math>  H_1: \mu = \mu_1 > \mu_0
</math>
que la media es de un valor mayor <math>\mu_1</math>.
Usando el lema de Neyman-Pearson la prueba mas poderosa de tamaño <math>\alpha</math> se
basa en la estadística de prueba:
 <math>  z = \dfrac{(\bar{x} - \mu_0)\sqrt{N}}{\sigma}
</math>
cuya region crítica es
 <math>  z > z_\alpha</math>
donde <math>z_\alpha</math> es el punto que divide el <math>100\alpha\%</math> superior del resto de
la distribucion normal unitaria <math>N(0,1)</math>. En términos de la media de las
observaciones originales debemos rechazar <math>H_0</math> en favor de <math>H_1</math> si
 <math>
  \bar{x} > \mu_0 + \dfrac{\sigma}{\sqrt{N}}z_\alpha
 </math>
y si no, entonces aceptar la hipótesis nula. Si la hipotesis alternativa
hubiera sido
 <math>  H_1': \mu = \mu_1 < \mu_0
</math>
entonces la misma estadística de prueba hubiera sido empleada, aunque la
región crítica hubiera sido
 <math>  z < -z_\alpha
</math>
o, equivalentemente
 <math>  \bar{x} < \mu_0 - \dfrac{\sigma}{\sqrt{N}}z_\alpha
</math>

Las pruebas anteriores y las hipótesis alternativas <math>H_1, H'_1</math> se conocen
como hipótesis '''unilaterales''', debido a que se indica claramente la
dirección del cambio de <math>\mu_0</math> a <math>\mu_1</math>. 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:
 <math>
  H_1^{''}: \mu = \mu_1 \neq \mu_0
</math>
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 <math>\alpha</math> es de
 <math>
  |z| > z_{\alpha/2}
</math>
donde <math>z_{\alpha/2}</math> es el punto que divide el <math>50\alpha\%</math> superior del resto
de la distribución normal unitaria. Equivalentemente, se rechaza <math>H_0</math> si
 <math>   \bar{x} > \mu_0 + \dfrac{\sigma}{\sqrt{N}}z_{\alpha/2} \;\;\;\; \text{ó}  \;\;\;\; \bar{x} < \mu_0 - \dfrac{\sigma}{\sqrt{N}}z_{\alpha/2}
</math>
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 <math>\textsf{R}</math>.
En la práctica uno selecciona un valor apropiado de <math>\alpha</math> y una o más
probabilidades <math>\beta</math> que sean tolerables. El tamaño de la muestra <math>N</math> se escoge entonces para
lograr que la <math>\beta</math> 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 
 <math>
  H_0: \mu = \mu_0 \; \; \; \; H_1: \mu > \mu_0
</math>
o
 <math>   H_0: \mu < \mu_0 \; \; \; \; H_1: \mu > \mu_0
</math>
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 <math>N</math> observaciones independientes
<math>x_1, \ldots, x_N</math> con media <math>\bar{x}</math> y varianza <math>s^2</math> lleva a la estadística
de prueba
 <math>  t = \dfrac{(\bar{x} - \mu_0)\sqrt{N}}{s}
</math>
Si <math>H_0</math> es verdadera, <math>t</math> tiene la distribución Student-Fisher con <math>N-1</math>
grados de libertad y rechazamos la hipótesis nula en una prueba con tamaño <math>\alpha</math> si
 <math>
  t > t_{\alpha;N-1}
</math>
donde <math>t_{\alpha; N-1}</math> es el punto que separa el <math>100\alpha\%</math> superior de la
distribución <math>t</math> definida en la sección anterior. De la misma forma, si la
hipótesis alternativa hubiera sido <math>H_1^{'}</math>, la región de rechazo
estaria definida por
 <math>
  t < -t_{\alpha;N-1}
</math>
y para la alternativa bilateral, \ref{H1biprime} la hipótesis nula será
rechazada si
 <math>  |t| > t_{\alpha/2; N-1}
</math>
La potencia de las pruebas que involucran la estadística <math>t</math> se puede calcular
de la rutina Pearson-Hartley con <math>\nu_1=1</math> y <math>\nu_2 = N-1</math> grados de libertad y
parametro de descentralidad
 <math>  \phi = \dfrac{|\mu - \mu_0|}{\sigma}\left(\dfrac{N}{2}\right)^{1/2}
</math>
Entonces se puede determinar, dado un <math>\alpha</math> 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 <math>\sigma^2</math> pero con
medias posiblemente distintas <math>\mu_1</math> y <math>\mu_2</math>. Denótense las observaciones
de las muestras como <math>x_1, \ldots, x_{N_1}, y_1, \ldots, y_{N_2}.</math> Entonces el 
criterio generalizado del cociente de posibilidades para la prueba de hipótesis
 <math>  H_0: \mu_1 = \mu_2
</math>
de medias poblacionales iguales contra la alternativa
 <math>  H_1: \mu_1 > \mu_2
</math> 
lleva a la estadística de prueba
 <math>  t =
\dfrac{\bar{x}-\bar{y}}{\sqrt{\dfrac{\sum^{N_1}_{i=1}(x_i-\bar{x})^2 +
\sum^{N_2}_{i=1}(y_i-\bar{y})^2}{N_1 + N_2 -2 }
\left(\dfrac{1}{N_1}+\dfrac{1}{N_2} \right) } }
</math>
Para una prueba con tamaño <math>\alpha</math> se rechaza <math>H_0</math> en favor de <math>H_1</math> si
 <math>  t > t_{\alpha; N_1+N_2-2}
</math>
Las regiones de rechazo para las otras hipótesis alternativas <math>H'_1: \mu_1 <
\mu_2</math> y <math>H''_1: \mu_1 \neq \mu_2</math> 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 <math>N</math> observaciones ha sido tomada de
alguna población con densidad continua <math>f(x;\theta)</math>. El intervalo de
confianza del <math>100(1-\alpha)\%</math> para <math>\theta</math> es aquel conjunto de
valores
 <math>  t_1 < \theta < t_2 </math>
cuyos límites fueron calculados a partir de los puntos porcentuales de la
función de distribución de la estimación <math> \hat{\theta}</math> de <math>\theta</math> tal que
 <math>  P(t_1 \leq \theta \leq t_2) = 1 - \alpha </math>
donde <math>\alpha</math> generalmente toma el valor 0.05 o alguna probabilidad más
pequeña. <math>1-\alpha</math> 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 <math>t_1, t_2</math> contenga a
<math>\theta</math> es de <math>1-\alpha</math>" 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 <math>t_1</math> y <math>t_2</math> de tal manera que la longitud de <math>t_2 - t_1</math> 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 <math>N</math> observaciones
independientes con media <math>\bar{x}</math> de la variante <math>N(\mu,\sigma^2)</math> con
varianza conocida, <math>H_0: \mu = \mu_0</math> se
rechaza en favor de
 <math>H_1: \mu = \mu_1 > \mu_0</math>
usando la prueba de tamaño <math>\alpha_1</math> cuando el valor de <math>\mu_0</math> está dado por
 <math>  \mu_0 = \bar{x} - \dfrac{\sigma}{\sqrt{N}}z_{\alpha_1}</math>
De la misma forma, el valor más grande de <math>\mu_0</math> para el cual, usando una
prueba con tamaño <math>\alpha - \alpha_1</math>, <math>H_0</math> es rechazada en favor de <math>H'_1</math>
cuando <math>H'_1: \mu = \mu_1 < \mu_0</math>, entonces ese valor de <math>\mu_0</math> será
 <math>
  \mu_0 = \bar{x} + \dfrac{\sigma}{\sqrt{N}}z_{\alpha -\alpha_1}
</math>
Donde <math>z_{\alpha_1}</math>y <math>z_{\alpha -\alpha_1}</math> son los porcentajes superiores
<math>100\alpha_1</math> y <math>100(\alpha - \alpha_1)</math> de la distribución normal unitaria.
Entonces el intervalo de confianza del <math>100(1-\alpha)</math> por ciento para <math>\mu</math> es
 <math>
 \bar{x} - \dfrac{\sigma}{\sqrt{N}}z_{\alpha_1} \leq \mu \leq \bar{x} +
 \dfrac{\sigma}{\sqrt{N}}z_{\alpha -\alpha_1}
 </math>
y su longitud es de <math>(\sigma/\sqrt{N})(z_{\alpha_1} + z_{\alpha -\alpha_1})</math>.
Pero se puede demostrar que la longitud mínima se obtiene cuando 
 <math>z_{\alpha_1}
 = z_{\alpha -\alpha_1}</math> 
o si
 <math>\alpha_1 = \frac{1}{2}\alpha</math>. 
Entonces el intervalo más corto con coeficiente <math>1-\alpha</math> tiene la forma simétrica
 <math>
 \bar{x} - \dfrac{\sigma}{\sqrt{N}}z_{\alpha/2} \leq \mu \leq \bar{x} +
 \dfrac{\sigma}{\sqrt{N}}z_{\alpha/2}
</math>

De la misma forma, si <math>\sigma^2</math> no es conocida, el intervalo de confianza del
<math>100(1-\alpha)</math> por ciento para <math>\mu</math> está dado por
 <math>
 \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}
</math>
donde <math>s</math> es la desviación estándar de la muestra, y <math>t_{\alpha/2; N-1}</math> es
el <math>50\alpha\%</math> superior de la distribución <math>t</math> con <math>N-1</math> 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 <math>100(1-\alpha)</math>
por ciento para el cambio atribuible a la condición experimental es
 <math>
  \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}
</math>
donde <math>\bar{x}</math> y <math>\bar{y}</math> son las medias de la primera y segunda muestras
compuestas de <math>N_1</math> y <math>N_2</math> observaciones y
 <math>
  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)}
</math>
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 <math>x_1, \ldots, x_N</math> constituyen
una muestra aleatoria de <math>N(\mu,\sigma^2)</math>, entonces la cantidad
 <math>
  \dfrac{(N-1)s^2}{\sigma^2} = \dfrac{\sum^{N}_{i=1}(x_i - \bar{x})^2}{\sigma^2}
</math> 
se encuentra distribuida como una variante ji-cuadrada con <math>N-1</math> grados de
libertad. El criterio generalizado del cociente de posibilidad para la prueba de
la hipótesis
 <math>
  H_0: \sigma^2 = \sigma_0^2
</math> 
contra la hipótesis alternativa
 <math>
  H_1: \sigma^2 > \sigma_0^2
</math> 
especifica que la región de rechazo para una prueba de tamaño <math>\alpha</math> es
 <math>
  \dfrac{(N-1)s^2}{\sigma_0^2} > \chi^2_{\alpha;N-1}
</math> 
donde <math>\chi^2_{\alpha;N-1}</math> es el <math>100\alpha</math> porciento superior de la
distribución ji-cuadrada con <math>N-1</math> grados de libertad. De la misma forma, las
regiones de rechazo para probar la hipótesis nula contra las alternativas
<math>H'_1: \sigma^2<\sigma_0^2</math> y <math>H''_1:\sigma^2 \neq \sigma_0^2</math> sería
 <math>
   \dfrac{(N-1)s^2}{\sigma_0^2} < \chi^2_{1-\alpha;N-1}
</math> 
y
 <math>
   \dfrac{(N-1)s^2}{\sigma_0^2} < \chi^2_{1-\alpha;N-1} \; \; o \; \;
\dfrac{(N-1)s^2}{\sigma_0^2} > \chi^2_{\alpha_2;N-1}
</math> 
respectivamente, donde en el segundo caso <math>\alpha_1 + \alpha_2 = \alpha</math>.
Estrictamente hablando, <math>\alpha_1</math> y <math>\alpha_2</math> 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 <math>\alpha</math>, pero para muchas aplicaciones con muestras grandes es
válido partir <math>\alpha</math> en dos partes iguales. Los intervalos de confianza para
<math>\sigma^2</math> se pueden obtener directamente de la región de rechazo \ref{RejChi}
como se demostrará más adelante.

Si las muestras independientes de <math>N_1</math> y <math>N_2</math> observaciones se han tomado
aleatoriamente de las poblaciones <math>N(\mu_1,\sigma_1^2)</math> y
<math>N(\mu_2,\sigma_2^2)</math>, entonces la hipótesis
 <math>
  H_0: \sigma^2 = \sigma_0^2
</math> 
se puede probar contra la hipótesis alternativa
 <math>
  H_1: \sigma^2 > \sigma_0^2
</math> 
por medio de la estadística
 <math>
  F = \dfrac{s_1^2}{s_2^2}
</math> 
Si la hipótesis nula es verdadera, la estadística tiene una distribución <math>F</math>
con <math>N_1-1, N_2-1</math> grados de libertad. La región de rechazo para una prueba de
tamano <math>\alpha</math> es
 <math>
  F > F_{\alpha; N_1-1, N_2-1}
</math> 
Por otro lado, <math>H_0</math> se puede probar contra la otra hipótesis alternativa
<math>H'_1: \sigma_1^2 < \sigma_2^2</math> por medio de la estadística
 <math>
  F = \dfrac{s^2_2}{s_1^2}
</math> 
cuya region crítica es
 <math>
  F > F_{\alpha_1; N_1-1, N_2-1} \ \  o' \ \  F < \dfrac{1}{ F_{\alpha_2; N_1-1,
N_2-1}}
</math> 
donde la repartición <math>\alpha_1 = \alpha_2 = \frac{1}{2}\alpha</math> 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 <math>\sigma^2</math> en cada cepa. Sea <math>\mu_j</math> la media poblacional
para la j-ésima cepa. Es posible entonces construir una prueba de
 <math>
  H_0: \mu_1 = \cdots = \mu_5
</math> 
contra la hipótesis alternativa de que algunas medias son distintas usando el
criterio generalizado del cociente de posibilidad. Aun mas, si <math>H_0</math> 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 <math>k</math> 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 <math>k</math> grupos) se puede expresar por
medio del modelo matemático
 <math>
  x_{ij} =  \mu + \tau_j + e_{ij} \ \ i=1, \ldots, N_j, j=1, \ldots, k.
</math> 
donde
 <math>
 \begin{align}
  \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{align}
 </math>
Los términos variantes <math>e_{ij}</math> 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 <math>\mu</math> en la cantidad <math>\tau_j</math>, 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, <math>\mu_j
\equiv \mu + \tau_j</math>. La hipótesis de los efectos equivalentes del tratamiento
se puede expresar como
 <math>
  H_0: \tau_1 = \cdots = \tau_k
</math> 
y se tomará como la alternativa a <math>H_0</math> el modelo general \ref{GenModel} para
las observaciones.

Este modelo se dice que tiene '''efectos fijos''', ya que los <math>\mu_j</math> son
parámetros y cualquier inferencia a partir de las observaciones sólo se puede
hacer con respecto a los <math>k</math> 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
{| border="0" cellpadding="5" cellspacing="0" align="center"
 |<math>x_{11}</math>
 |<math>\cdots</math>
 |<math>x_{1k}</math>
 |-
 |<math>\cdots</math>
 |<math>\cdots</math>
 |<math>\cdots</math>
 |-
 |<math>x_{N_11}</math>
 |<math>\cdots</math>
 |<math>x_{N_kk}</math>
|}
Denótese el total de las observaciones para el j-ësimo tratamiento como
 <math>
  T_j = \sum^{N_j}_{i=1} x_{ij}
</math> 
y se denotará la media de ese tratamiento como <math>\bar{x}_j = T_j/N_j</math>. La suma
de todas las observaciones será denotada por
 <math>
  G = T_1 + \cdots + T_k
</math> 
el número total de unidades experimentales por
 <math>
  N = N_1 + \cdots + N_k
</math> 
y la media general por <math>\bar{x} = G/N</math>. Expresamos entonces la suma total de
los cuadrados
 <math>
 \begin{align}
 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{align}
 </math>
como la suma dos componentes independientes
 <math>S = SST + SSE</math>
donde
 <math>
 \begin{align}
 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{align}
 </math>
y
 <math>
 \begin{align}
 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{align}
 </math>

{| border="1" cellpadding="5" cellspacing="0" align="center"
 |Fuente
 |Suma de cuadrados
 |Grados de Libertad
 |Cuadrado Medio
 |-
 |Tratamiento
 |SST
 |<math>k-1</math>
 |<math>\dfrac{SST}{k-1}</math>
 |-
 |Intra-tratamiento
 |SSE
 |<math>N-k</math>
 |<math>\dfrac{SSE}{N-k}</math>
 |-
 |Total
 |S
 |<math>N-1</math>
 |
|}

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 <math>H_0</math> es
 <math>
  F = \dfrac{N-k}{k-1}\dfrac{SSE}{SST}
</math> 
y, cuando <math>H_0</math> es verdadera, <math>SST/(k-1)\sigma^2</math>, <math>SSE/(N-k)\sigma^2</math> son
variantes independientes con distribucion ji-cuadrada con <math>k-1</math> y <math>N-k</math>
grados de libertad, respectivamente. Entonces la estadística de prueba tiene
una distribución <math>F</math> con <math>k-1</math> y <math>N-k</math> grados de libertad. Rechazamos <math>H_0</math>
con una prueba de nivel <math>\alpha</math> si
 <math>
  F > F_{\alpha; k-1, N-k}
</math> 
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 <math>\nu_1 = k-1</math> y <math>\nu_2 = N-k</math>, 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
 <math>
  \phi = \dfrac{\sqrt{\sum_{j=1}^k N_j \left(\tau_j - \dfrac{\sum
N_j \tau_j}{N}\right)^2}}{\sigma \sqrt{k}}
</math> 
o, como en el caso de muestras con tratamientos equivalentes que lo que
generalmente se encuentra en problemas de diseño experimental
 <math>
  \phi = \dfrac{\sqrt{N_j \sum_{j=1}^k \tau_j^2}}{\sigma \sqrt{k}}
</math> 
En esta última expresión, el término de corrección ha desaparecido de la
habitual restricción <math>\tau_1 + \ldots +  \tau_k =0</math> que se impone sobre los
efectos del tratamiento. El parámetro <math>\phi</math> 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 <math>\alpha</math>. 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 <math>\tau_j</math> del modelo de una vía como cualquier función
lineal
 <math>
  \sum_{j=1}^{k} c_j\tau_j
</math> 
cuyos coeficientes tienen la propiedad de
 <math>
  \sum_{j=1}^{k} c_j = 0
</math> 
Por lo tanto <math>\tau_1 - \tau_2</math> y <math>3\tau_1 +(\tau_2+\tau_3+\tau_4)</math> son
contrastes mientras que <math>\tau_2 - (\tau_3 + \tau_4)</math> 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 <math>1-\alpha</math> para todos los contrastes <math>\tau_j</math> tienen la forma

 <math>
  \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}}
</math> 

donde todas las sumas son sobre los <math>k</math> tratamientos. <math>\sum c_j \bar{x}_j</math> es la
estimación muestral del contraste <math>\sum c_j \tau_j</math> y
 <math>
  s = \sqrt{\dfrac{SSE}{N-k}}
</math> 
Nótese que <math>s^2 +\sum c_j^2N_j</math> es la estimación de la varianza del contraste
estimado. Aceptamos la hipótesis nula
 <math>
  H_0: \sum_{j=1}^k c_j \tau_j =0
</math> 
al nivel <math>\alpha</math> si el intervalo simultáneo de confianza para ese contraste
incluye el valor cero. Si, por otro lado
 <math>
   \sum c_j\bar{x}_j > s\sqrt{(k-1)F_{\alpha;k-1,N-k}\sum \dfrac{c_j^2}{N_j}}
</math> 
o
 <math>
  \sum c_j\bar{x}_j < -s\sqrt{(k-1)F_{\alpha;k-1,N-k}\sum \dfrac{c_j^2}{N_j}}
</math> 
rechazamos <math>H_0</math> en favor de la hipótesis alternativa unilateral <math>H_1: \sum c_j \tau_j > 0</math> o <math>H'_1: \sum c_j \tau_j<0</math>. El nivel conjunto de todas las pruebas es de <math>\alpha</math>.

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 <math>m</math>
declaraciones de intervalos de confianza <math>H_1,\ldots, H_m</math> de los parámetros del modelo lineal. La probabilidad de que
<math>H_i</math> sea una declaración verdadera (esto es, que <math>H_i</math> contenga el valor de la i-ésima función paramétrica) es
igual a <math>P(H_i)</math> y la probabilidad de que todas las declaraciones sean correctas simultáneamente es de <math>P(H_1 \cap
\cdots \cap H_m)</math>, donde la notación de intersección <math>H_i \cap H_j</math> denota el evento donde ambos <math>H_i</math> y <math>H_j</math> 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 <math>1-\alpha</math>, donde <math>\alpha</math> 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 <math>H_i</math>. En lugar de calcular esta probabiliad, debemos contentarnos con la
siguiente cota inferior:
 <math>
 P(H_1 \cap \cdots H_m) \geq 1 - \sum_{i=1}^m P(\bar{H}_i).
</math> 
donde <math>P(\bar{H}_i) = 1 - P(H_i)</math> 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 <math>m</math> intervalos de confianza
simultáneos, por lo general se le asigna a cada declaración una proporción de error de <math>\alpha</math>, de tal manera que
el coeficiente para toda la familia es de al menos <math>1-\alpha</math>.

Los intervalos de confianza de Bonferroni sobre <math>m</math> contrastes dados
 <math>
 \psi_i = \sum_{j=1}^k c_{ij}\tau_j \; \; \;  i=1,  \ldots, m
</math> 
de los parámetros del modelo lineal de una vía son
 <math>
 \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}}
</math> 
La regla de decisión para probar <math>H_0: \Psi_i=0</math> es equivalente a aquella para el metodo de Scheffe: si
(\ref{Bonf01}) no incluye el valor cero, se rechaza <math>H_0</math>. Para todas las comparaciones apareadas de los tratamientos,
<math>m = \tfrac{1}{2}k(k-1)</math>, mientras que para comparaciones sucesivas (suponiendo que hubo algun ordenamiento a priori),
<math>m=k-1</math>. Si <math>m</math> 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 <math>1-\alpha</math>.
Para valores de <math>m</math> 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:
{| border="1" cellpadding="5" cellspacing="0" align="center"
|
 |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
 <math>
T_1=21 \; \; \; T_2=9 \; \; \; T_3 =6 \; \; \; G=36 \; \; \; \sum_{j=1}^3 \sum_{i=1}^6 x_{ij}^2 = 106
</math>
La suma de cuadrados debida a los medicamentos es
 <math>
SST = \frac{1}{6}(21^2 + 9^2 + 6^2) - \dfrac{36^2}{18} = 21
</math>
y la suma total de cuadrados es <math>S =   106 -72 = 34</math>. Los valores se resumen en la tabla de análisis de varianza:
{| border="1" cellpadding="5" cellspacing="0" align="center"
  |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 <math>F_{0.01;2,15} = 6.36</math> 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 <math>F</math> observamos que la <math>F</math> excede el valor crítico
para <math>\alpha =0.001</math>

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:
 <math>
\sqrt{(k-1)F_{0.01;2.15}} = 3.567
</math>
Primero compárense los medicamentos B y C calculando el intervalo de confianza para <math>\tau_2 - \tau_3</math>. Aqui <math>c_1=0,c_2=1, c_3=-1</math> 
y la desviación estándar estimada de la población <math>\bar{x_2} -\bar{x_3}</math> es de
 <math>
\sqrt{0.867(\frac{1}{6}+\frac{1}{6})} =0.5376
</math>
El intervalo de confianza es
<math>
-1.42 \leq \tau_1 - \tau_2 \leq 2.42
</math>
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
 <math>
0.08 \leq \tau_1 - \tau_2 \leq 3.92
</math>
y es posible aceptar la hipótesis alternativa
 <math>
H_1: \tau_1 > \tau_2
</math>
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 <math>c_1 =1, \ c_2 = -\frac{1}{2}, \ c_3 = -\frac{1}{2}</math> y la desviación estándar estimada de
ese contraste es
 <math>
\sqrt{0.867(\frac{1}{6}+\frac{1}{6}+\frac{1}{24})} = 0.4653
</math>
La estimación del contraste es de 2.25, y el intervalo de confianza es de
 <math>
0.59 \leq \tau_1 - \frac{1}{2}(\tau_2 +\tau_3) \leq 3.91
</math>
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
{| cellpadding="5" cellspacing="0" align="center"
 |0.13 
| <math>\leq \tau_1 - \tau_2 \leq</math>
 | 3.87
 |-
 | -1.37
| <math> \leq \tau_2 - \tau_3 \leq</math>
| 2.37
 |-
 |0.63
| <math> \leq \tau_1- \tau_3 \leq </math> 
| 4.37
|}

= Introducción al lenguaje R =
== ¿Qué es R? ==
El sistema <math>\textsf{R}</math> para cómputos estadísticos es un ambiente para análisis de datos y gráficos. El origen de
<math>\textsf{R}</math> es el lenguaje <math>\textsf{S}</math> el cual fue desarrollado por John Chambers y sus colegas en los laboratorios Bell por la década de los 1960s. El lenguaje <math>\textsf{S}</math> 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 <math>\textsf{R}</math> 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 <math>\textsf{R}</math> se encuentra localizada en el siguiente URL:
 [http://www.R-project.org http://www.R-project.org]
Todos los recursos de <math>\textsf{R}</math> 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 <math>\textsf{R}</math> 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 <math>\textsf{R}</math> 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 <math>\textsf{R}</math> Archive Network" (CRAN), el cual se encuentra en
 [[The R Project|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, <math>\textsf{R}</math> 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:

<pre>
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.

>
</pre>

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 <math>\textsf{R}</math> 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 <math>\textsf{R}</math> 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
<math>\textsf{R}</math> 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 <math>\textsf{R}</math> que tiene el nombre de variable <tt>x</tt>. El operador de asignación <tt><-</tt>
afianza el valor a su lado derecho con el nombre de la variable a su lado izquierdo. El valor del objeto <tt>x</tt>
puede ser inspeccionado escribiendo

  R> x
    [1] 7

esto llama implícitamente al método <tt>print</tt>

  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 <math>\textsf{R}</math>.
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 <math>\textsf{R}</math>: 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 <math>\textsf{R}</math>. 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:
<pre>
  > 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"
</pre>

== 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 <math>\textsf{R}</math> 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 <math>\textsf{R}</math>
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 <math>\textsf{R}</math>.

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
<pre>
  , , 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

</pre>
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
<pre>
   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"
</pre>
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: <math>4 \times 4 \times 2</math>. 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 <math>\textsf{R}</math> usa para almacenar datos. Escribamos el siguiente mando:
  R> str(USArrests)
El sistema responde con:
<pre>
  '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 ...
</pre>
<math>\textsf{R}</math> 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". <ref name="dataframe">[Voy a dejar
por un lado el purismo del idioma castellano y me estaré refiriendo, cuando sea necesario, al marco de datos como
"dataframe".]</ref> Los objetos tipo dataframe son los objetos más importantes dentro de <math>\textsf{R}</math>. 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''':
<pre>
  R> dim(USArrests)
  [1] 50  4
</pre>
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''':
<pre>
  R> nrow(USArrests)
  [1] 50
  R> ncol(USArrests)
  [1] 4
</pre>
Los nombres de las variables que describen las observaciones, esto es, los nombres de las columnas del dataframe se
extraen usando
<pre>
  R> names(USArrests)
[1] "Murder"   "Assault"  "UrbanPop" "Rape"
</pre>
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
<pre>
  R> length(USArrests[,"Murder"])
  [1] 50
</pre>
Un vector es la estructura más elemental para manejar datos en <math>\textsf{R}</math> 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:
<pre>
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
</pre>
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:
<pre>
  R> class(USArrests[,"Murder"])
  [1] "numeric"
</pre>
El tamaño del vector donde se encuentran almacenados
<pre>
  R> length(USArrests[,"Murder"])
  [1] 50
</pre>
Y el primer elemento del vector es
<pre>
  R> USArrests[,"Murder"][1]
  [1] 13.2
</pre>
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:
<pre>
  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 ...
</pre>
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.
<pre>
  R> nlevels(InsectSprays[,"spray"])
  [1] 6
  R> levels(InsectSprays[,"spray"])
  [1] "A" "B" "C" "D" "E" "F"
</pre>
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
<pre>
  R> table(InsectSprays[,"spray"])

 A  B  C  D  E  F
12 12 12 12 12 12
</pre>
<references />

== 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 <math>\textsf{R}</math>. 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:
<pre>
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"
</pre>
Usaremos entonces la función '''read.table''' para leer el archivo:
<pre>
R> users <- read.table("users.csv", header=TRUE, sep=",")
</pre>
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 <math>\textsf{R}</math>:
<pre>
  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 ...
</pre>
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 <math>\textsf{R}</math> 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:
<pre>
  R> write.table(users, file="USERS.csv", sep=",", col.names=NA)
</pre>
Si se desea guardar los datos para que solo sean procesados por <math>\textsf{R}</math>, cualquier objeto en <math>\textsf{R}</math> se puede
salvar en un archivo binario externo a través de:
<pre>
  R> save(users, file="USERS.rda")
</pre>
donde la extensión '''.rda''' es estándar. Si deseamos cargar los objetos de un archivo tipo '''rda''' a
<math>\textsf{R}</math>.
  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 <math>\textsf{R}</math>. Internamente, un dataframe es una lista de vectores con una longitud común <math>n</math>, donde <math>n</math> 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''':
<pre>
  R> occ <- users[,"Occupation"]
  R> occ[1]
  [1] technician
</pre>
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:
<pre>
  R> occ[125:128]
  [1] lawyer    lawyer    none      marketing
  21 Levels: administrator artist doctor educator engineer ... writer
</pre>
<math>\textsf{R}</math> 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 <math>\textsf{R}</math> 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:
<pre>
  R> occ[-(6:943)]
  [1] technician other      writer     technician other
</pre>
Y en forma convencional con índices positivos:
<pre>
  R> occ[1:5]
  [1] technician other      writer     technician other
</pre>
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:
<pre>
  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
</pre>
Con esta instrucción, le estamos diciendo a <math>\textsf{R}</math> 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:
<pre>
  R> occ <- users$Occupation
  R> occ[1:5]
  [1] technician other      writer     technician other
  21 Levels: administrator artist doctor educator engineer ... writer
</pre>
La instrucción mostrada arriba es equivalente a la instrucción
<pre>
  R> occ <- users[,"Occupation"]
</pre>

'''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'''.

<pre>
  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
</pre>
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:
<pre>
  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>
</pre>
donde la expresion '''users$Age <15''' indica un vector lógico de longitud 943 donde
<pre>
R> table(users$Age<15)

FALSE  TRUE
  932    11
</pre>
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
<pre>
 R> na_occ <- is.na(users$Occupation)
 R> table(na_occ)
na_occ
FALSE
  943
</pre>
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
<pre>
 R> mayores50 <- subset(users, Age> 50)
 R> dim(mayores50)
[1] 105   5
</pre>
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 <math>\textsf{R}</math>: '''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.
<pre>
  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
</pre>
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 <math>\textsf{R}</math>.
Sin embargo, el verdadero poder de <math>\textsf{R}</math> 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 <math>\textsf{R}</math> 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 <math>\textsf{R}</math>

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:
<pre>
R> iqr <- function(x) {
       q <- quantile(x, prob = c(0.25, 0.75), names = FALSE)
	return(diff(q))
}
</pre>
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
<pre>
R> xdata <- rnorm(100)
R> iqr(xdata)
[1] 0.9977416
R> IQR(xdata)
[1] 0.9977416
</pre>
Sin embargo, cuando el vector numérico contiene valores faltantes, la función falla como se muestra en el siguiente
ejemplo:
<pre>
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

</pre>
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 
<pre>
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
</pre>
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. [[Archivo:Fig2 1.jpg|thumb|Histogramas de edades]] 

Las relaciones bivariantes de dos variables continuas se grafican generalmente como una gráfica de puntos. En <math>\textsf{R}</math>, 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
<pre>
  R> fm <- ingresos ~ ventas
  R> class(fm)
   [1] "formula"
</pre>
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 <math>\textsf{R}</math> para obtener una gráfica de cajas de las edades es:
<pre>
R> tmp <- subset(users, Gender %in% c("M", "F"))
R> plot(Age ~ Gender, data = tmp, ylab="Edades", xlab="Genero", varwidth=TRUE)
</pre>

Podemos observar el resultado en la figura 7 [[Archivo:Fig2 2.jpg|thumb|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
<math>\textsf{R}</math>, 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 <math>\textsf{R}</math> usando el mando '''source''':
<pre>
  R> source("analisis.R", echo = TRUE)
</pre>
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
<math>
[X_1, \ldots, X_p]
</math>
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 <math>p</math> 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 <math>[X_1, X_2, X_3, X_4]</math> con alguna distribución en el espacio de 4 dimensiones. Si esos pesos fueron registrados con una muestra de <math>N</math> ratones, las observaciones se pueden resumir por medio del arreglo
 <math>
 \begin{bmatrix}
  x_{11} & \cdots & x_{14}  \\
 \cdots  & \cdots & \cdots \\
  x_{N1} & \cdots & x_{N4}
 \end{bmatrix}
 </math>
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
 <math>
A =
\begin{bmatrix}
	a_{11} & \cdots & a_{1c} \\
	\cdots & \cdots & \cdots \\
	a_{r1} & \cdots & a_{rc}
\end{bmatrix}

 </math>
como un arreglo rectangular ordenado de los elementos. El término general de la matriz se expresará como <math>a_{ij}</math>
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 <math>r</math> filas y <math>c</math> columnas será referida como una
matriz de <math>r \times c</math> 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 <math>1 \times 1</math>. 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 <math>n</math>
componentes como
 <math>
	x =
	\begin{bmatrix}
		x_1\\
		\vdots \\
		x_n
		\end{bmatrix}

 </math>
De la misma forma, el ''vector fila'' lo expresaremos como
 <math>
	\textbf{x}^{\prime} = [ x_1, \ldots, x_n]

 </math>
y consiste de una sola fila de <math>n</math> elementos. Cada vector especifica las coordenadas de un punto en el espacio
euclideano <math>n</math>-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 <math>\textbf{x}^\prime</math> es la ''transpuesta'' del
vector
columna '''x'''. En general, si '''A''' es cualquier matriz de <math>r \times c</math> dimensiones, la transpuesta de
'''A''' es la matriz <math>\textbf{A}^\prime</math> de <math>c \times r</math> dimensiones formada al intercambiar las filas y las
columnas:
 <math>
	\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}

 </math>
Si una matriz es cuadrada e igual a su transpuesta, se dice que es ''simétrica''. Entonces <math>a_{ij} = a_{ji}</math> para
todos los pares de <math>i</math> y <math>j</math>. Por ejemplo
 <math>
	\textbf{A} =
	\begin{bmatrix}
		3 & 0 & -1\\
		0 & 1 & 2\\
		-1 & 2 & -4
	\end{bmatrix}

 </math>
es simétrica, mientras que
 <math>
	\textbf{B} =
	\begin{bmatrix}
		2 & 3 \\
		1 & 2
	\end{bmatrix}

 </math>
no lo es. Por brevedad, se omitirán con frecuencia los elementos inferiores duplicados. Los elementos <math>a_{ii}</math> 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
 <math>
	\text{tr}\; \textbf{A} = \sum_{i=1}^n a_{ii}

 </math>
Ciertas matrices son particularmente importantes. La matriz ''identidad''
 <math>
	\textbf{I} =
	\begin{bmatrix}
		1 & \cdots & 0 \\
		\cdots & \cdots & \cdots \\
		0 & \cdots & 1
	\end{bmatrix}

 </math>
es una matriz cuadrada con unos en la diagonal principal y ceros en las posiciones restantes. La ''matriz diagonal''
de <math>p \times p</math>
 <math>
	\textbf{D}(a_i) =
	\begin{bmatrix}
		a_1 & \cdots & 0 \\
		\cdots & \cdots & \cdots \\
		0 & \cdots & a_p
	\end{bmatrix}

 </math>
tiene los elementos <math>a_1, \ldots, a_p</math> en las posiciones de la diagonal principal y ceros en las posiciones restantes.
Algunos de los elementos <math>a_i</math> pueden ser cero. De aquí en adelante el término diag('''A''') = '''D'''(<math>a_{ii}</math>)
denotará la matriz diagonal formada por la matriz cuadrada '''A'''.
Una ''matriz triangular'' de <math>p \times p</math> tiene el patrón
 <math>
	\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}

 </math>
La ''matriz nula'' de <math>r \times c</math>
 <math>
	\begin{bmatrix}
		0 & \cdots & 0 \\
		\cdots & \cdots & \cdots \\
		0 & \cdots & 0
	\end{bmatrix}

 </math>
tiene un cero en cada una de sus posiciones. A veces necesitaremos el vector
 <math>
	\textbf{j}^\prime = [ 1, \ldots, 1 ]

 </math>
y la matriz
 <math>
	\textbf{E} =
	\begin{bmatrix}
		1 & \cdots & 1 \\
		\cdots & \cdots & \cdots \\
		1 & \cdots & 1
	\end{bmatrix}

 </math>
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 <math>r \times c</math> '''A''' y '''B''' se dice que son iguales si y sólo si
 <math>
	a_{ij} = b_{ij}
 </math>
''Suma''. La suma de dos matrices de iguales dimensiones es la matriz de la suma de los elementos correspondientes.
Si
 <math>
	\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}
 </math>
entonces
 <math>
	\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}
 </math>
Una matriz se puede sustraer de otra con iguales dimensiones formando la matriz de diferencias de los elementos
individuales. Entonces
 <math>
\begin{align}
	\mathbf{A + B} &= \mathbf{B + A}\\
	\mathbf{A +(B +C)} &= \mathbf{(A +B) + C} \\
	\mathbf{A - (B -C)} &= \mathbf{A - B + C}
\end{align}
 </math>
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 <math>c</math> multiplicando cada elemento
de '''A''' por <math>c</math>:
 <math>
	c\textbf{A} =
	\begin{bmatrix}
		ca_{11} & \cdots & ca_{1k} \\
		\cdots & \cdots & \cdots \\
		ca_{r1} & \cdots & ca_{rk}
	\end{bmatrix}
 </math>
''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 <math>p \times r</math>, y
'''B''' tiene dimensiones de <math>r \times q</math>, entonces el <math>ij</math>-ésimo elemento del producto <math>\textbf{C = AB}</math> se calcula
como
 <math>
c_{ij} = \sum^r_{k=1}a_{ik}b_{kj}
 </math>
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 <math>p \times q</math>. Por ejemplo, si
 <math>
	\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}
 </math>
entonces
 <math>
	\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}
 </math>
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:
 <math>
\begin{align}
	\mathbf{A(B+C)} &= \mathbf{AB + AC} \\
	\mathbf{A(BC)} &= \mathbf{AB(C)}
\end{align}
 </math>
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
 <math>
	\mathbf{A} =
	\begin{bmatrix}
		1 & 2\\
		2 & 3
	\end{bmatrix} \;\;
	\mathbf{B} =
	\begin{bmatrix}
		3 & -1 \\
		-1 & 1
	\end{bmatrix}
 </math>
Entonces
 <math>
	\mathbf{AB} =
	\begin{bmatrix}
		1 & 1 \\
		3 & 1
	\end{bmatrix}
	\;\;
	\mathbf{BA} =
	\begin{bmatrix}
		1 & 3\\
		1 & 1
	\end{bmatrix}
 </math>
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 <math>d_1, \ldots, d_r</math> posee el efecto de multiplicar cada elemento en la
i-ésima fila por <math>d_i</math>:
 <math>
	\mathbf{D}(d_i)\mathbf{A} =
	\begin{bmatrix}
		d_1a_{11} & \cdots & d_1a_{1c} \\
		\cdots & \cdots & \cdots \\
		d_ra_{r1} & \cdots & d_ra_{rc}
	\end{bmatrix}
 </math>
La posmultiplicación por una matrix diagonal similar de <math>c \times c</math> multiplica cada elemento en la j-ésima columna
por <math>d_j</math>.

Una matriz puede considerarse como una transformación lineal de los vectores en un espacio a vectores en otro espacio.
Si '''x''' tiene <math>m</math> componentes y '''y''' tiene <math>n</math> componentes, es posible expresar una transformación del
sistema de coordenadas en <math>m</math> dimensiones de los elementos de '''x''' al espacio <math>n</math>-dimensional de
'''y''' en forma matricial como
 <math>
	\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}
 </math>
La transformación a un tercer conjunto de variables especificado por el vector <math>\textbf{z}</math> de <math>p</math> componentes puede
representarse como
 <math>
\begin{align}
	\mathbf{z} &= \mathbf{By} \\
	&= \mathbf{BAx}
\end{align}
 </math>
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 <math>x_{ij}</math> es la medición del <math>i</math>-ésimo individuo en
<math>ji</math>-ésima característica, y <math>\bar{x}_j</math> y <math>s_j</math> 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'':
<math>
z_{ij} = \dfrac{z_{ij} - \bar{x}_j}{s_j}
</math>
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 <math>N \times p</math>
 <math>
	\textbf{X} =
	\begin{bmatrix}
		x_{11} & \cdots & x_{1p} \\
		\cdots & \cdots & \cdots \\
		x_{N1} & \cdots & x_{Np}
	\end{bmatrix}
 </math>
y formamos la matriz diagonal
 <math>
	\textbf{D}=
	\begin{bmatrix}
		\dfrac{1}{s_1} & \cdots & 0 \\
		\cdots & \cdots & \cdots \\
		0 & \cdots & \dfrac{1}{s_p}
	\end{bmatrix}
 </math>
a partir de las desviaciones estándard de la muestra. Si introducimos la matriz de <math>N \times N</math>
 <math>
	\textbf{E} =
	\begin{bmatrix}
		1 & \cdots & 1\\
		\cdots & \cdots & \cdots \\
		1 & \cdots & 1
	\end{bmatrix}
 </math>
con un uno en cada posición, entonces la matriz '''Z''' de <math>N \times p</math> dimensiones correspondiente a las mediciones
normalizadas se puede calcular como
 <math>
	\textbf{Z} = \left( \textbf{I}- \frac{1}{N}\textbf{E}\right) \textbf{XD}
 </math>
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:
 <math>
	\textbf{x}^\prime \textbf{y}= [x_1, \ldots, x_p]
	\begin{bmatrix}
		y_1\\
		\vdots\\
		y_p
	\end{bmatrix}
 </math>
Ya que el producto interno es escalar, entonces <math>\textbf{x}^\prime \textbf{y} = \textbf{y}^\prime \textbf{x}</math>. 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 <math>p</math> dimensiones. En general, la distancia
entre los puntos con coordenadas dadas por '''x''' y '''y''' es
 <math>
	d = \left[ \sum^p_{i=1} (x_i -y_i)^2 \right]^{\frac{1}{2}}
 </math>
El coseno <math>\theta</math> del ángulo entre los vectores '''x''' y '''y''' es
 <math>
	\cos \theta = \frac{\textbf{x}^\prime \textbf{y}}{(\textbf{x}^\prime \textbf{x})^{\frac{1}{2}}(\textbf{y}^\prime
	\textbf{y})^{\frac{1}{2}}}
 </math>
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
 <math>
\begin{align}
	\textbf{x}^\prime &= [1,0,0] \\
	\textbf{y}^\prime &= [0,1,0] \\
	\textbf{z}^\prime &= [0,0,1]
\end{align}
 </math>
se pueden interpretar como si fueran los ejes coordenados. El vector
 <math>
\textbf{e}^\prime = [1,1,1]
 </math>
hace un ángulo de <math>\sqrt{3}/3</math> con <math>\textbf{x}^\prime, \textbf{y}^\prime</math> y <math>\textbf{z}^\prime</math> y caracteriza la
''línea equiangular'' en el espacio tridimensional. Las cantidades <math>\sqrt{3}/3</math> se conocen como los ''cosenos
direccionales} de la línea. Similarmente, los vectores
 <math>
\begin{align}
	\textbf{u}^\prime &= [1,1,0] \\
	\textbf{v}^\prime &= [1, -1, 0]
\end{align}
 </math>
forman un ángulo de 90 grados en el plano <math>xy</math> del espacio. Los vectores <math>\textbf{u}^\prime</math> y <math>\textbf{v}^\prime</math>
forman cada uno un ángulo de 45 grados con <math>\textbf{x}^\prime</math>, y sus ángulos con <math>\textbf{y}^\prime</math> 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:
 <math>
	(\textbf{AB})^\prime = \textbf{B}^\prime \textbf{A}^\prime
 </math>
En general, si <math>\textbf{A}_1, \ldots, \textbf{A}_k</math> son conformables,
 <math>
	(\textbf{A}_1, \cdots \textbf{A}_k)^\prime = \textbf{A}_k^\prime \cdots \textbf{A}_1^\prime
 </math>
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 <math>A</math> tiene el valor
 <math>
	\det \textbf{A} \equiv |\textbf{A}| = \sum^k_{i=1} a_{ij}C_{ij}
 </math>
donde se supone que no hay suma sobre el índice <math>j</math> y donde <math>C_{ij}</math> es el ''cofactor'' de <math>a_{ij}</math> definido por
 <math>
	C_{ij} \equiv (-1)^{i+j}M_{ij}
 </math>
y <math>M_{ij}</math> es el <math>\emph{menor}</math>de la matriz <math>\textbf{A}</math> formado al eliminar la fila
<math>i</math> y la columna <math>j</math> de <math>\textbf{A}</math>. 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 <math>(1,
\ldots, n) </math> tomando cada permutación como los subíndices de las letras <math>a, b, \ldots</math> y sumando con los signos dados
por <math>\epsilon_p = (-1)^{i(p)}</math> donde <math>i(p)</math> es el número de ''inversiones de la permutación'' en la permutación
<math>p</math>, y <math>\epsilon_{n_1 n_2}</math> es el símbolo de permutación. Por ejemplo, para <math>n=3</math>, 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 <math>3 \times 3</math> es
 <math>
	\begin{vmatrix}
		a_1 & a_2 & a_3\\
		b_1 & b_2 & b_3\\
		c_1 & c_2 & c_3\\
	\end{vmatrix} = a_1b_2c_3 - a_1b_3c_2 -a_2b_1c_3+a_2b_3c_1+a_3b_1c_2 -a_3b_2c_1.
 </math>
Para determinantes mayores de <math>3 \times 3</math> 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 <math>\textsf{R}</math> como calcular
determinantes.

Ciertas propiedades de los determinantes nos serán útiles en los capítulos subsiguientes:
# 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.
# Si los elementos de una sola fila o columna de la matriz \textbf{A} de <math> n \times n</math> son multiplicados por el escalar <math>c</math>, el determinante de la nueva matriz es igual a <math>c|\textbf{A}|</math>. Si cada elemento es multiplicado por <math>c</math>, entonces <math>|c\textbf{A}| = c^n|\textbf{A}</math>.
# Si dos filas o columnas de una matriz cuadrada son intercambiadas, entonces el signo del determinante se invierte.
# 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.
# El determinante de una matriz no cambia agregando un múltiplo de alguna columna a otra columna. Lo mismo vale para las filas.
# Si todos los elementos de una fila o una columna son iguales a cero, entonces el determinante es cero.
# Si <math>\textbf{A}</math> y <math>\textbf{B}</math> son matrices de <math>n \times n</math>, el determinante de <math>\textbf{AB}</math> es igual al producto de los determinantes individuales.
# 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 <math>\textbf{A}</math> es la matriz única <math>\textbf{A}^{-1}</math> con elementos tales que
 <math>
	\textbf{AA}^{-1} = \textbf{A}^{-1}\textbf{A} = \textbf{I}
 </math>
Es posible que <math>\textbf{A}^{-1}</math> no exista, en la misma forma que no es posible dividir entre cero con escalares. Entonces
se dice que <math>\textbf{A}</math> es <math>\emph{singular}</math>. Las matrices cuyas inversas existen se conocen como 
<math>\emph{no singulares}</math>.

Los elementos de <math>\textbf{A}^{-1}</math> pueden ser calculadas a partir de dos resultados de la sección anterior. Se forma la
matriz de cofactores <math>C_{ij}</math>:
 <math>
	\textbf{C} =
	\begin{bmatrix}
		C_{11} & \cdots & C_{1n}\\
		\cdots & \cdots & \cdots \\
		C_{n1} & \cdots & C_{nn}
	\end{bmatrix}
 </math>
conocida como la <math>\emph{adjunta}</math> de <math>\textbf{A}</math> Entonces el producto interno de la <math>i</math>-ésima fila de <math>\textbf{C}</math> y la
<math>h</math>-ésima fila de <math>\textbf{A}</math> es igual a <math>|\textbf{A}|</math> si <math>i=j</math> y es igual a cero si <math>i \neq h</math>. Si tomamos la
transpuesta de <math>\textbf{C}</math> y dividimos cada elemento entre <math> |\textbf{A}|</math>, tendremos entonces la inversa deseada:
 <math>
	\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}
 </math>
Es evidente de esta última expresión que la matriz inversa existe si y sólo si el determinante de <math>\textbf{A}</math> es
distinto de cero. El cálculo de <math>\textbf{A}^{-1}</math> 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:

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

'''Ejemplo: ''' La inversa de
 <math>
	\textbf{A} =
	\begin{bmatrix}
		2 & 0& 0& 0\\
		0 & 4 & 0 & 0\\
		0 & 0 & 3 & 2\\
		0 & 0 & 1 & 2
	\end{bmatrix}
 </math>
es
 <math>
	\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}
 </math>

== El Rango de una Matriz ==
Dos vectores de <math>p</math> 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
<math>
\textbf{x}^\prime = [1, 0, -1] \;\;\; \textbf{y}^\prime =[4,0,-4]
</math>
son linealmente dependientes, mientras que
<math>
\textbf{u}^\prime = [2,-1,0,7] \;\;\; \textbf{v}^\prime = [6,2,0,0]
</math>
son linealmente independientes. Un conjunto de <math>k</math> 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 <math>\textbf{x}_1, \ldots, \textbf{x}_k</math> forman un conjunto linealmente independiente si no existen escalares
<math>c_1, \ldots, c_k</math> tales que para algun vector <math>\textbf{x}_i</math> en el conjunto se da que
<math>
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\;.
</math>
Los vectores
<math>
\textbf{x}^\prime = [1,-1,2] \;\;\; \textbf{y}^\prime =[2,0,-1] \;\;\; \textbf{z}^\prime =[0,-2,5]
</math>
constituyen un conjunto linealmente dependiente ya que <math>\textbf{z}^\prime = 2\textbf{x}^\prime  - \textbf{y}^\prime</math>,
mientras que los vectores unitarios
<math>
\textbf{t}^\prime = [1,0,0] \;\;\; \textbf{u}^\prime = [0,1,0] \;\;\; \textbf{v}^\prime = [0,0,1]
</math>
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 <math>N</math> 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 <math>k</math> de vectores no excede la dimensionalidad <math>p</math> de los
vectores. Entonces el rango de
 <math>
	\textbf{X}=
	\begin{bmatrix}
		\textbf{x}^\prime \\
		\vdots\\
		\textbf{x}^\prime _k
	\end{bmatrix}
 </math>
es el número de vectores fila linealmente independientes de la matriz. El rango puede variar desde cero para una matriz
nula, hasta <math>k</math> para una matriz de ''rango completo''. Si, por otro lado, el número de filas de <math>\textbf{X}</math> 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 <math>\textbf{A}</math> ''tiene un rango igual a r si contiene al menos un menor distinto de cero de'' <math>r
\times r</math> ''dimensiones y ningún menor distinto de cero con dimensionalidad mayor que r''.

'''Ejemplo: ''' La matriz
 <math>
\begin{bmatrix}
	1 & 2 & 3 &4\\
	3 & 6 & 9 & 12\\
	4 & 3 & 2 & 1\\
	-1 & 3 & 7 & 11\\
	8 & 6 & 4 & 2
\end{bmatrix}
 </math>
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:
# El rango de <math>\textbf{A}^\prime</math> es igual al rango de <math>\textbf{A}</math>. Esto viene de la equivalencia de las definiciones de rango para filas y columnas.
# El rango de <math>\textbf{A}^\prime\textbf{A}</math> es igual al de <math>\textbf{A}</math>. Similarmente, el rango de <math>\textbf{A}\textbf{A}^\prime</math> es equivalente al rango de <math>\textbf{A}</math>.
# El rango de <math>\textbf{A}</math> no cambia si se pre- o posmultiplica <math>\textbf{A}</math> 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:
# Intercambio de dos filas o dos columnas.
# La multiplicación de todos los elementos de una fila o una columna por una constante escalar.
# 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
 <math>
\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}
 </math>
se pueden usar para llevar a cabo transformaciones elementales sobre las primeras dos filas de
 <math>
\textbf{A}=
\begin{bmatrix}
	5 & 1\\
	-2 & 3\\
	3 & 2
\end{bmatrix}
 </math>
dando así
 <math>
\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}
 </math>

Usemos ahora las operaciones elementales de fila para calcular el rango de la matriz
 <math>
\begin{bmatrix}
	1 & 0 & 2 \\
	-3 & 2 & 1\\
	-1 & 2 & 5\\
	3 & 0 & 6
\end{bmatrix}
 </math>
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 <math>A</math> en la siguiente forma:
 <math>
\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}
 </math>
Ya que las filas restantes son linealmente independientes, se concluye que <math>A</math> es de rango igual a dos.

Ambas operaciones de fila y columna se pueden aplicar simultáneamente a cualquier matriz de <math>m \times n</math> para
reducirla a su forma canónica
 <math>
	\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}
 </math>
la cual contiene un uno en las primeras <math>r</math> posiciones diagonales y cero en todas las demás. Esto quiere decir que para
cada matriz <math>A</math> existen matrices elementales de transformación de filas y columnas llamadas <math>\textbf{R}_1
\ldots, \textbf{R}_p, \textbf{C}_1, \ldots, \textbf{C}_q</math> tales que
 <math>
\begin{align}
	\textbf{R}_p \cdots \textbf{R}_1\textbf{A}\textbf{C}_1\cdots \textbf{C}_q &= \textbf{PAQ}\\
	&= \textbf{F}
\end{align}
 </math>
es la forma canónica de <math>A</math>. Pero si <math>\textbf{A}</math> es cuadrada y de rango completo <math>n</math>, entonces las operaciones de
fila producirán
 <math>
	\textbf{R}_p \cdots \textbf{R}_1\textbf{A} = \textbf{I}
 </math>
esto significa que
 <math>
	\textbf{A}^{-1} = \textbf{R}_p \cdots \textbf{R}_1\textbf{I}
 </math>
Por lo tanto, la misma secuencia de operaciones elementales de fila que reduce <math>A</math> a la matriz identidad
transforma la matriz identidad de <math>n \times n</math> a la inversa de <math>A</math>. Esta propiedad de las transformaciones
elementales nos proporciona con una herramienta poderosa para invertir matrices: simplemente se expresa la matriz de <math>n
\times 2n</math>
 <math>
	\begin{bmatrix}
		\textbf{A}& & \textbf{I}
	\end{bmatrix}
 </math>
y se aplican operaciones elementales a cada fila con el fin de transformar las primeras <math>n</math> columnas en la matriz
identidad. Cuando esto ha sido logrado, aparecerá <math>\textbf{A}^{-1}</math> en las últimas <math>n</math> columnas.


'''Ejemplo:''' La transformaciones sucesivas en la inversión de
 <math>
\textbf{A}=
\begin{bmatrix}
	1 &1 &2\\
	2 & 3 & 3\\
	3 & -1 & 8
\end{bmatrix}
 </math>
se muestran a continuación:
 <math>
\begin{bmatrix}
	1 & 1 & 2 & 1 & 0 & 0\\
	2 & 3 & 3 & 0 & 1 & 0\\
	3 & -1 & 8 & 0 & 0 & 1
\end{bmatrix}
 </math>
 <math>
\begin{bmatrix}
	1 & 1 & 2 & 1 & 0 & 0\\
	0 & 1 &-1 & -2 & 1 & 0\\
	0 & -4 & 2 & -3 & 0 & 1
\end{bmatrix}
 </math>
 <math>
\begin{bmatrix}
	1 & 0 & 3 & 3 & -1 & 0\\
	0 & 1 & -1 & -2 & 1 & 0\\
	0 & 0 & -2 & -11 & 4 & 1
\end{bmatrix}
 </math>
 <math>
\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}
 </math>
y
 <math>
\textbf{A}^{-1} =
\begin{bmatrix}
	-13.5 & 5 & 1.5 \\
	3.5 & -1 & -0.5 \\
	5.5 & -2 & -0.5
\end{bmatrix}
 </math>
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
 <math>
\textbf{A}=
\begin{bmatrix}
	10 & 5 & -2 & 0 \\
	6 & 3 & 2 & 1 \\
	4 & 5& 12 & 3 \\
	5 & 1 & 3 & 8
\end{bmatrix}
 </math>
Empezamos restando múltiplos apropriados de la primera fila de las filas restantes:
 <math>
\begin{bmatrix}
	10 & 5 & -2 & 0 \\
	0 & 0 & 3.2 & 1 \\
	0 & 3 & 12.8 & 3 \\
	0 & -1.5 & 4 & 8
\end{bmatrix}
 </math>
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:
 <math>
\begin{bmatrix}
	10 & 5 &-2 & 0 \\
	0 & -3 & -12.8 & -3 \\
	0 & 0 & 3.2 & 1 \\
	0 & -1.5 & 4 & 8
\end{bmatrix}
 </math>
Una operación de fila pivoteando en la segunda columna transforma la matriz anterior en:
 <math>
\begin{bmatrix}
	10 & 5 & -2 & 0 \\
	0 & -3 & -12.8 & -3 \\
	0 & 0 & 3.2 & 1 \\
	0 & 0 & 10.4 & 9.5
\end{bmatrix}
 </math>
Finalmente tenemos
 <math>
\begin{bmatrix}
	10 & 5 & -2 & 0 \\
	0 & -3 & -12.8 & -3 \\
	0 & 0 & 3.2 & 1 \\
	0 & 0 & 0 & 6.25
\end{bmatrix}
 </math>
y el determinante es:
 <math>
\begin{align}
|\textbf{A}| &= (10)(-3)(3.2)(6.25)\\
 &=-600
\end{align}
 </math>

== 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 <math>\textbf{G}</math> la cual
satisface
 <math>
	\textbf{AGA} = \textbf{A}
 </math>
donde <math>\textbf{A}</math> es de <math>p \times q</math> dimensiones y de rango <math>r</math>, mientras que <math>\textbf{G}</math> tiene que ser necesariamente de
<math>q \times p</math> dimensiones y del mismo rango que <math>\textbf{A}</math>. Existen otros tipos de inversas generalizadas, las que
satisfacen las siguientes condiciones:
 <math>
	\textbf{GAG} = \textbf{G} \;\;\; (\textbf{GA})^\prime = \textbf{GA} \;\;\; (\textbf{AG})^\prime = \textbf{AG}
 </math>
Nuestra matriz <math>\textbf{G}</math> se puede calcular a partir de la reducción canónica de <math>\textbf{A}</math>. Expresamos
 <math>
	\textbf{PAQ} =
	\begin{bmatrix}
		\textbf{D}& \textbf{0}\\
		\textbf{0}& \textbf{0}
	\end{bmatrix}
 </math>
donde <math>\textbf{D}</math> es una matriz diagonal de <math>r \times r</math> (que no es necesariamente la identidad) y las matrices nulas
tienen las dimensiones apropiadas. Las matrices <math>\textbf{P}</math> y <math>\textbf{Q}</math> son los productos
respectivos definidos anteriormente.
Presentamos ahora la matriz <math>\textbf{F}^-</math> de <math>q \times p</math>:
 <math>
	\textbf{F}^- =
	\begin{bmatrix}
		\textbf{D}^{-1}& \textbf{0}\\
		\textbf{0}& \textbf{0}
	\end{bmatrix}
 </math>
Entonces <math>\textbf{G}= \textbf{QF}^-\textbf{P}</math> es una inversa generalizada de la matriz <math>\textbf{A}</math> ya que uno puede
verificarlo reemplazando <math>\textbf{A}</math> en la expresión por su representación
<math>\textbf{P}^{-1}\textbf{FQ}^{-1}</math>. La matriz <math>\textbf{G}</math> no es única. a menos que <math>\textbf{A}</math> sea cuadrada y no singular.

<math>\textbf{Ejemplo:}</math> Vamos a calcular la inversa generalizada de la matriz <math>\textbf{A}</math> vista en un ejemplo anterior. Las
matrices de operadores de fila que reducen la tercera y la cuarta fila a los vectores nulos son:
 <math>
\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}
 </math>
Entonces
 <math>
\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}
 </math>
Las transformaciones de columna que reducen a <math>\textbf{PA}</math> a la forma diagonal canónica son
 <math>
\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}
 </math>
y
 <math>
\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}
 </math>
Entonces
 <math>
\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}
 </math>
y
 <math>
\textbf{G} =
\begin{bmatrix}
	1 & 0 & 0 & 0\\
	\frac{3}{2} & \frac{1}{2} & 0 & 0\\
	0 & 0 & 0 & 0
\end{bmatrix}
 </math>
es una inversa generalizada de <math>\textbf{A}</math>.

== Ecuaciones Lineales Simultáneas ==
El conjunto de ecuaciones con las incógnitas <math>x_1, \ldots, x_n</math>
 <math>
 \begin{align}
	a_{11}x_1 + \cdots + a_{1n}x_n &= c_1\\
	\vdots \\
	a_{n1}x_1 + \cdots + a_{nn}x_n &= c_n
 \end{align}
 </math>
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
 <math>
	\textbf{Ax} = \textbf{c}
 </math>
donde <math>\textbf{A}</math> es la matriz de coeficientes de <math>m \times n</math>, <math>\textbf{x}^\prime = [x_1, \ldots, x_n]</math>, y
<math>\textbf{c}^\prime = [c_1, \ldots, c_m]</math> El sistema general <math>\textbf{Ax=c}</math> de <math>m</math> ecuaciones en <math>n</math> incógnitas tiene
una solución si y solo si la ''matriz aumentada'' de <math>m \times n</math>
	<math>
	\left[ \textbf{A}\;\;\textbf{c}\right]
	</math>
tiene el mismo rango <math>r</math> que <math>\textbf{A}</math>. Si no es así, entonces se dice que el sistema es inconsistente. Es esencial
distinguir entre sistemas homogéneos, para los cuales <math>\textbf{c}=0</math>, 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 <math>\textbf{A}^{-1}</math> existe, la solución única del sistema es
 <math>
	\textbf{x} = \textbf{A}^{-1}\textbf{c}
 </math>
Esto viene de premultiplicar ambos lados de la ecuación por <math>\textbf{A}^{-1}</math>.

<math>\textbf{Ejemplo: }</math> La matriz inversa del sistema
 <math>
\begin{align}
	x_1 + x_2 - x_3 & =1 \\
	-x_1+x_2+x_3 &= -1\\
	x_1 -x_2+x_3 &= 1
\end{align}
 </math>
es
 <math>
\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}
 </math>
La solución es
 <math>
\begin{align}
	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{align}
 </math>

=== Sistemas No Homogéneos: A <math>m \times n</math> y Rango <math>r</math> ===
Si se satisface la condición sobre los rangos de las matrices aumentadas y de coeficientes, se puede encontrar una
solución seleccionando cualesquier <math>r</math> ecuaciones linealmente independientes y obteniendo <math>r</math> de las <math>n</math> incógnitas en
términos de las constantes <math>c_i</math> y de las <math>n-r</math> variables restantes.

<math>\textbf{Ejemplo:}</math> Las ecuaciones
 <math>
	\begin{align}
		x_1 + x_2 &= 2 \\
		x-1 + x_2 &= 1
	\end{align}
 </math>
se pueden interpretar geométricamente como dos líneas paralelas en el plano <math>x_1x_2</math>. Debido a que no intersectan, el
sistema no puede tener solución. El rango de la matriz de coeficientes
 <math>
\begin{bmatrix}
	1 & 1 \\
	1 & 1
\end{bmatrix}
 </math>
es uno, mientras que el de la matriz aumentada
 <math>
	\begin{bmatrix}
		1 & 1 & 2\\
		1 & 1 & 1
	\end{bmatrix}
 </math>
es dos.

<math>\textbf{Ejemplo:}</math> El sistema
 <math>
	\begin{align}
		2x_1 + 3x_2 - x_3 &= 1 \\
		x_2 + x_3 &= 2
	\end{align}
 </math>
tiene matrices aumentadas y de coeficientes con rango igual a dos. Se expresan las ecuaciones como
 <math>
	\begin{align}
		2x_1 + 3x_2 - x_3 &= 1 \\
		 x_2 &= 2-x_3
	\end{align}
 </math>
y se obtiene <math>x_1</math> por sustitución. La solución general está dada por el vector
 <math>
	[\frac{1}{2}(-5+4x_3),2-x_3,x_3]
 </math>
donde <math>x_3</math> 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 <math>[\textbf{A} \;
\textbf{0}]</math> es el mismo que el de <math>\textbf{A.}</math> Todos los sistemas homogéneos se ven satisfechos por la solución
trivial <math>\textbf{x}^\prime = [0, \ldots, 0]</math>, y si el rango de la matriz de coeficientes es igual al número de
incógnitas <math>n</math>, é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 <math>r</math> de <math>\textbf{A}</math> es estrictamente menor que <math>n</math>, y siempre será posible encontrar
<math>n-r</math> 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.

<math>\textbf{Ejemplo:}</math> El sistema
 <math>
	\begin{align}
		3x_1 - 2x_2 &= 0\\
		5x_1 + x_2 &= 0
	\end{align}
 </math>
define dos líneas en el plano <math>x_1x_2</math> que intersectan en el origen. Por lo tanto la única solución es la trivial:
<math>[0,0]</math>.

<math>\textbf{Ejemplo:}</math> En el sistema
 <math>
	\begin{align}
		x_1 - x_2 + 2x_3 &= 0\\
		x_1 + 3x_2 - 2x_3 &=0\\
		3x_1 + x_2 + 2x_3 &= 0
	\end{align}
 </math>
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 <math>x_3</math> resolviendo el sistema
 <math>
	\begin{align}
		x_1 - x_2 &= -2x_3\\
		x_1 + 3x_2 &= 2x_3
	\end{align}
 </math>
La solución general del conjunto original es, en unidades de <math>x_3</math>
 <math>
\begin{align}
	x_1 &= -x_3 \\
	x_2 &= x_3
\end{align}
	</math>

<math>\textbf{Ejemplo:}</math> El sistema
 <math>
\begin{align}
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{align}
 </math>
también es de rango dos. La primera y tercera ecuaciones se pueden resolver convenientemente en términos de <math>x_3</math> y
<math>x_4</math> para dar el vector solución
 <math>
[ \frac{1}{2}(-x_3 +x_4), \frac{1}{2}(3x_3+x_4), x_3, x_4]
 </math>
Si hacemos <math>x_3</math> igual a cero y <math>x_4 </math>igual a uno, y alternativamente hacemos <math>x_3</math> igual a uno y <math>x_4</math> igual a cero,
tendremos los siguientes vectores solución linealmente independientes:
 <math>
[\frac{1}{2}, \frac{1}{2},0,1]
 </math>
 <math>
[-\frac{1}{2},\frac{3}{2},1,0]
 </math>
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
 <math>
	\begin{align}
		x_1 &= \frac{1}{2}(a-b)\\
		x_2 &= \frac{1}{2}(a+3b)\\
		x_3 &= b\\
		x_4 &= a
	\end{align}
 </math>
para <math>a</math> y <math>b</math> arbitrarios.

=== Solución por medio de inversas generalizadas ===
Se ha demostrado que el sistema <math>\textbf{Ax}=\textbf{c}</math> de <math>m</math> ecuaciones consistentes en <math>n</math> incógnitas tiene la
solución
 <math>
\textbf{x}^{*} = \textbf{Gc} + ( \textbf{GA - I})\textbf{z}
 </math>
donde <math>\textbf{G}</math> es la inversa generalizada de <math>\textbf{A}</math> definida por , <math>\textbf{I}</math> es la matriz identidad
de <math>n \times n</math>, y <math>\textbf{z}</math> es cualquier vector de <math>n \times 1</math> constantes arbitrarias. Hacemos notar que el término
que contine a <math>\textbf{z}</math> desaparece cuando <math>\textbf{G} = \textbf{A}^{-1}</math>.

<math>\textbf{Ejemplo:}</math> Obtendremos soluciones a los sistemas de los ejemplos X y Y por el método de las inversas
generalizadas. Para el primer sistema, <math>\textbf{P} = \textbf{I}</math> y
 <math>
\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}
 </math>
 <math>
\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}
 </math>
La solución más general es
 <math>
\textbf{x}^\prime = [-1/2(5+z_3), 2+z_3, -z_3]
 </math>
Para el ejemplo Y
 <math>
\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}
 </math>
 <math>
\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}
 </math>
Si hacemos que <math>z_3</math> sea igual a <math>-x_3 </math>la solución general es <math>\textbf{x}^\prime = [-x_3, x_3, x_3]</math>.

=== 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
 <math>
	\cos \theta = \frac{\textbf{x}^\prime \textbf{y}}{(\textbf{x}^\prime \textbf{x})^{\frac{1}{2}}(\textbf{y}^\prime
	\textbf{y})^{\frac{1}{2}}}

 </math>
para el coseno del ángulo entre los vectores <math>\textbf{x}</math> y <math>\textbf{y}</math> 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'' <math>\textbf{T}</math> es una
matriz cuadrada cuyas files son un conjunto de vectores ortonormales. Por lo tanto
 <math>
 \begin{align}
	\textbf{TT}^\prime &= \textbf{T}^\prime\textbf{T}\\
	&= \textbf{I}
\end{align}
 </math>
y la inversa de <math>\textbf{T}</math> es simplemente su transpuesta <math>\textbf{T}^\prime</math>. Una matriz ortogonal sencilla de <math>2 \times
2</math> es
 <math>
\begin{bmatrix}
	0 & 1\\
	1 & 0
\end{bmatrix}
 </math>
Esta es la matriz de transformación de una reflexión de los puntos en el plano <math>xy</math> alrededor de la línea de 45
grados. Una matriz ortogonal más general de esa dimensión es
 <math>
	\textbf{T}=
	\begin{bmatrix}
		\cos \theta & \sin \theta\\
		-\sin \theta & \cos \theta
	\end{bmatrix}
 </math>
y se interpreta como la matriz de transformación para la rotación de los ejes coordenados <math>xy</math> en un ángulo <math>\theta</math>.
Esto es, si <math>x</math> y <math>y</math> fueran las coordenadas de un punto bajo los ejes originales, el punto después de haber girado los
ejes, tendría las coordenadas
 <math>
\begin{align}
	u &= x \cos \theta + y \sin \theta \\
	v &= -x \sin \theta + y \cos \theta
\end{align}
 </math>
donde los nuevos ejes coordenados han sido girados un ángulo de <math>\theta</math> 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
 <math>
\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}
 </math>
es una matriz ortogonal de <math>3 \times 3</math>.

La matriz ortogonal <math>\textbf{T}</math> de <math>3 \times 3</math> 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 <math>n</math> ejes coordenados alrededor del origen. Las
distancias en este espacio permanecen invariantes ante esta rotación, ya que si ejecutamos la transformación
 <math>
\textbf{y} = \textbf{Tx}
 </math>
entonces
 <math>
\begin{align}
\textbf{x}^\prime \textbf{x} &= \textbf{y}^\prime \textbf{T}^\prime \textbf{Ty}\\
&= \textbf{y}^\prime\textbf{Iy}
\end{align}
 </math>
y las suma de los cuadrados será invariante.
Las matrices ortogonales tienen estas útiles propiedades:
# Las columnas de una matriz ortogonal son ortogonales.
# El determinante de una matriz ortogonal siempre será 1 ó -1. Ya que <math>|\textbf{TT}^\prime| =
		|\textbf{I}| = 1</math> y <math>|\textbf{T}|= |\textbf{T}^\prime|</math>, <math>|\textbf{T}| =
		\pm 1\,</math>.
# 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 <math>x_1, \ldots, x_n</math> es una expresión del tipo
 <math>
	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_1x_2 + \cdots + 2a_{2n}x_1x_n + \cdots 
	+ 2a_{n-1,n}x_{n-1}x_n = \sum_{i=1}^n\sum_{j=1}^n a_{ij}x_ix_j
 </math>
donde <math>a_{ij} = a_{ji}</math>. Algunas de las <math>a_{ij}</math> pueden ser iguales a cero. Hacemos notar de inmediato que la forma cuadrática se puede expresar en notación matricial como
<math>
\textbf{x}^\prime \textbf{Ax}
</math>
donde <math>\textbf{x}^\prime = [x_1, \ldots, x_n] </math>y <math>\textbf{A}</math> es la matriz simétrica de <math>n \times n</math> de coeficientes. La forma cuadrática más sencilla es simplemente <math>f(x) = a_{11}x^2</math>, la ecuación de una parábola en la variable <math>x</math>. Las formas cuadráticas en donde las <math>x_i</math> 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
 <math>
	\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
 </math>
se puede expresar como una forma cuadrática en términos de las observaciones <math>x_i</math> usando la matriz
 <math>
	\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}
 </math>

Una matriz simétrica <math>\textbf{A}</math> y su forma simétrica asociada se llaman ''definidas positivas'' si
<math>\textbf{x}^\prime \textbf{Ax} > 0</math> para toda <math>\textbf{x}</math> no nula. Si <math>\textbf{x}^\prime \textbf{Ax}\leq 0</math> 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 <math>n</math> variables a la forma
 <math>
	d_1y_1^2 + \cdots + d_n y_n^2
 </math>
la cual solo contiene cuadrados de las nuevas variables <math>y_i</math> y con coeficientes <math>d_i > 0</math>. En forma similar, una forma positiva semidefinida cuadrática se puede reducir a
 <math>
	d_1y_1^2 + \cdots + d_ry_r^2
 </math>
donde todos los coeficientes son positivos, y <math>r \leq n</math> 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
 <math>
	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}|
 </math>
Si <math>\textbf{A}</math> es de rango <math>r</math> se dice que es regular si <math>p_r \neq 0</math> y no hay <math>p_i</math> 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 <math>\textbf{A}</math> es una matriz arreglada regularmente,
# Una condición necesaria y suficiente para la ''definitez'' positiva es que <math>p_i >0</math> para todo <math>i=1, \ldots ,n</math> .
# Una condición necesaria y suficiente para la ''semidefinitez'' positiva es que <math>p_1 >0, \ldots, p_r>0</math> y los <math>n-r\,\,p_i</math> restantes sean iguales a cero, donde <math>r</math> puede ser igual a <math>n</math>.


<math>\textbf{Ejemplo: }</math> La matriz
 <math>
\begin{bmatrix}
 4 & 6 & 0\\
 6 & 9 & 0 \\
 0 & 0 & 0
\end{bmatrix}
 </math>
tiene la secuencia de determinantes principales menores <math>p_0=1</math>, <math>p_1=4</math>, <math>p_2=0</math>,<math>p_3=0</math>, y por eso no se encuentra
regularmente arreglada. Podemos reordenar las filas y las columnas para obtener la nueva matriz:
 <math>
\begin{bmatrix}
	2 & 0 & 0\\
	0 & 4 & 6\\
	0 & 6 & 9
\end{bmatrix}
 </math>
la cual tiene la secuencia de determinantes <math>p_0=1</math>, <math>p_1=2</math>, <math>p_2=8</math>, <math>p_3=0</math>. La matriz y su forma cuadrática son positivas semidefinidas con rango igual a dos.


<math>\textbf{Ejemplo: }</math> 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 <math>a</math> y elementos <math>b</math> que se encuentran fuera de la
diagonal. El determinante de una matriz así de <math>i \times i</math> dimensiones es
 <math>
(a-b)^{i-1}[a+(i-1)b]
 </math>
de tal forma que los determinantes menores principales de son
 <math>
p_i = \frac{1}{N}[N-1-(i-1)]
 </math>
Los primeros <math>N-1</math> determinantes son positivos, mientras que el <math>N</math>-ésimo siempre será cero. Entonces la suma de las desviaciones cuadradas es una forma cuadrática positiva semidefinida de rango <math>N-1</math>.

== Las raíces y vectores característicos de una matriz ==
Las ''raíces características'' de la matriz <math>\textbf{A}</math> de <math>p \times p</math> son las soluciones de la ecuación
determinante
 <math>
	| \textbf{A} - \lambda \textbf{I} | = 0
 </math>
El determinante es un polinomio de <math>p</math>-ésimo grado en la variable <math>\lambda</math>, y por eso <math>\textbf{A}</math> tiene <math>p</math> raíces
características. El desarrollo de Laplace del ''determinante característico'' nos permite expresar el
''polinomio característico'' como
 <math>
	| \textbf{A} - \lambda \textbf{I} | = (-\lambda)^p + S_1(-\lambda)^{p-1} + S_2(-\lambda)^{p-2} + \cdots
	-S_{p-1}\lambda + |\textbf{A}|
 </math>
donde <math>S_i</math> es la suma de todos los <math>i \times i</math> determinantes menores. <math>S_1</math> es simplemente la suma de los elementos
diagonales de <math>\textbf{A}</math>, o <math>\text{tr} \textbf{A}</math>. De la teoría de ecuaciones polinomiales se puede obtener que:
# El producto de las raíces características de <math>\textbf{A}</math> es igual a <math>| \textbf{A} |</math>.
# La suma de las raíces características de <math>\textbf{A}</math> es igual a la ''traza'' de <math>\textbf{A}</math>.



<math>\textbf{Ejemplo: }</math> La matriz
 <math>
\textbf{A} =
\begin{bmatrix}
	2 & 1 & 1\\
	1 & 2 & 1 \\
	1 & 1 & 2
\end{bmatrix}
 </math>
tiene una traza de 6 y tres determinantes menores principales de <math>2 \times 2</math> cada uno igual a
 <math>
\begin{vmatrix}
	2 & 1\\
	1 & 2
\end{vmatrix}
= 3
 </math>
mientras que <math>|\textbf{A}|=4</math>. La ecuación característica de <math>\textbf{A}</math> es
 <math>
\lambda^3 - 6 \lambda^2 + 9\lambda - 4 = 0
 </math>
y sus raíces son 1, 1, y 4.


De aquí en adelante usaremos frecuentemente las siguientes propiedades de las raíces características:
# Las raíces características de una matriz simétrica con elementos reales son todas reales.
# Las raíces características de una matriz positiva definida son todas positivas.
# Si una matriz de simétrica de <math>n \times n</math> es positiva semidefinida con rango <math>r</math>, contiene exactamente <math>r</math> raíces características positivas y <math>n-r</math> raíces cero.
# Las raíces características distintas de cero del producto <math>\textbf{AB}</math> son iguales a las raíces distintas de cero del producto <math>\textbf{BA}</math>. Por lo tanto las trazas de <math>\textbf{AB}</math> y <math>\textbf{BA}</math> son iguales.
# Las raíces características de una matriz diagonal son los elementos diagonales mismos.

Con cada raíz característica <math>\lambda_i</math> de la matriz cuadrada <math>\textbf{A}</math> se encuentra asociado un ''vector característico'' <math>\textbf{x}_i</math> cuyos elementos satisfacen el sistema homogéneo de ecuaciones
 <math>
	\left[ \textbf{A} - \lambda_i \textbf{I}\right]\textbf{x}_i = \textbf{0}
 </math>
Por definición de la raíz característica el determinante del sistema desaparece, y siempre existe una solución no
trivial <math>\textbf{x}_i</math>. 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:
# Si <math>\lambda_i</math> y <math>\lambda_j</math> son raíces características distintas de la matriz simétrica <math>\textbf{A}</math>, sus vectores asociados <math>\textbf{x}_i</math> y <math>\textbf{x}_j</math> serán ortogonales. Esto es evidente si premultiplicamos las definiciones de los vectores
 <math>
		\textbf{Ax}_i = \lambda_i\textbf{x}_i \;\;\;\;\; 	\textbf{Ax}_j = \lambda_j\textbf{x}_j
</math>
por <math>\textbf{x}_j ^\prime</math> y <math>\textbf{x}_i^\prime</math> respectivamente. Pero esto implica que
 <math>
 \lambda_i \textbf{x}_j ^\prime\textbf{x}_i = \lambda_j \textbf{x}_i ^\prime \textbf{x}_j
 </math>
y ya que <math>\lambda_i \neq \lambda_j</math> podemos concluír que los vectores son ortogonales.
# Para cada matriz simétrica real <math>\textbf{A}</math> existe una matriz ortogonal <math>\textbf{P}</math> tal que
 <math>
			\textbf{P}^\prime\textbf{AP} = \textbf{D}
 </math>
donde <math>\textbf{D}</math> es la matriz diagonal de las raíces características de <math>\textbf{A}</math>. Los vectores característicos normalizados de <math>\textbf{A}</math> pueden tomarse como las columnas de <math>\textbf{P}</math>. 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
 <math>
	\textbf{x} = \textbf{Py}
 </math>
a las <math>p</math> variables en la forma cuadrática <math>\textbf{x}^\prime\textbf{Ax}</math>, entonces la forma se convierte en
 <math>
\begin{align}
	\textbf{x}^\prime \textbf{Ax} &= \textbf{y}^\prime\textbf{P}^\prime\textbf{APy}	\\
	&= \textbf{y}^\prime\textbf{Dy}\\
	&= \lambda_1y_i^2 + \cdots + \lambda_ry_r^2
\end{align}
 </math>
donde la <math>\lambda_i</math> son las raíces características de la matriz de coeficientes, y <math>r</math> 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.


<math>\textbf{Ejemplo: }</math> El vector característico de la raíz más grande del ejemplo anterior debe satisfacer las ecuaciones
 <math>
\begin{align}
	-2x_{11} + x_{12} + x_{13} &= 0 \\
	x_{11} - 2x_{12} + x_{13} &= 0 \\
	x_{11} + x_{12} -2x_{13} &=0
\end{align}
 </math>
Si fijamos <math>x_{13}</math> arbitrariamente igual a uno y resolvemos las primeras dos ecuaciones no homogéneas, el vector
característico será
 <math>
\textbf{x}^\prime = [1,1,1]
 </math>
y por supuesto cualquier vector no nulo <math>[a,a,a]</math> es un vector característico para <math>\lambda_3=4</math>. El vector asociado
con la doble raíz <math>\lambda_1 = \lambda_2 =1</math> debe satisfacer el sistema con la ecuación linealmente independiente
 <math>
x_{21} + x_{22} + x_{23} = 0
 </math>
La solución más general a esto es el vector
 <math>
[a, b, -a-b]
 </math>
y seleccionando adecuadamente <math>a</math> y <math>b</math> se pueden generar dos vectores característicos linealmente independientes para
la raíz doble. Nótese que ''cualquier'' elección de <math>a</math> y <math>b</math> dará vectores ortogonales al primer vector
<math>[1,1,1]</math>.
Podemos elegir <math>a</math>y <math>b</math> para continuar esta ortogonalidad. Por ejemplo,
 <math>
\begin{align}
	\textbf{x}_1^\prime &= [0,1,-1] \\
	\textbf{x}_2^\prime &= [-2, 1,1]
\end{align}
</math>

<math>\textbf{Ejemplo: }</math> Las matrices cuadradas con la propiedad de
 <math>
\textbf{AA} = \textbf{A}
 </math>
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 <math>r</math> términos cuadrados. Es fácil verificar que la matriz es idempotente y
que es posible realizar la transformación
 <math>
\sum_{i=1}^N (x_i - \bar{x})^2 = y_1^2 + \cdots + y_{N-1}^2
 </math>
a una suma de los cuadrados de <math>N-1</math> 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 <math>\textsf{R}</math>.

== 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
 <math>
	\textbf{A} =
	\begin{bmatrix}
		\textbf{A}_{11} & \cdots & \textbf{A}_{1n} \\
		\cdots & \cdots & \cdots \\
		\textbf{A}_{m1} & \cdots & \textbf{A}_{mn}
	\end{bmatrix}
 </math>
de submatrices <math>\textbf{A}_{ij}</math> que contienen <math>r_i</math> filas y <math>c_j</math> columnas, donde es obvio que todas las submatrices en
una fila dada de <math>\textbf{A}</math> 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
<math>\textbf{A}</math> y <math>\textbf{B}</math> cuyas submatrices tienen dimensiones similares es
 <math>
	\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}
 </math>
El producto de las matrices partidas <math>\textbf{A}</math> y <math>\textbf{B}</math> es
 <math>
	\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}
 </math>
Las dimensiones dentro de cada producto de las submatrices debe conformarse; si las submatrices de <math>\textbf{A}</math> tienen el
número de columnas <math>c_1, \ldots, c_n</math>, el número de filas de <math>\textbf{B}</math> debe ser <math>c_1, \ldots, c_n</math>.

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 <math>f(\textbf{x})</math> una función continua de los elementos del vector
<math>\textbf{x}^\prime = [x_1, \ldots, x_p]</math> cuya primera y segunda derivadas 
 <math>
\frac{\partial f(\textbf{x})}{\partial x_i} \;\;\;\; \frac{\partial^2f(\textbf{x})}{\partial x_i \partial x_j}
 </math>
existen para todos los puntos \textbf{x} en alguna región del espacio euclideanto de <math>p</math> dimensiones. Se define el
operador vectorial de derivada parcial como 
 <math>
	\frac{\partial}{\partial\textbf{x}} = 
	\begin{bmatrix}
		\frac{\partial}{\partial x_i} \\
		\vdots \\
		\frac{\partial}{\partial x_p}
	\end{bmatrix}
 </math>
Aplicando este operador a <math>f(\textbf{x})</math> se obtiene el vector de derivadas parciales
 <math>
	\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}
 </math>
Las siguientes funciones y sus derivadas son especialmente importantes:
# Si <math>f(\textbf{x})</math> es constante para toda \textbf{x},
 <math>
			\frac{\partial f(\textbf{x})}{\partial\textbf{x}} = 
			\begin{bmatrix}
				0 \\
				\vdots\\
				0
			\end{bmatrix}
</math>

# <math>f(\textbf{x}) = \textbf{a}^\prime \textbf{x}</math>:
		 <math>
			\frac{\partial f(\textbf{x})}{\partial\textbf{x}} = 
			\begin{bmatrix}
				a_1 \\
				\vdots\\
				a_p
			\end{bmatrix}
		 </math>
		Nótese que la forma columnar del vector de derivadas no cambia si expresamos <math>f(\textbf{x}) = \textbf{x}^\prime\textbf{a}</math>.
# Para calcular el vector de derivadas de la forma cuadrática <math>\textbf{x}^\prime\textbf{AX}</math> se expresa como
	<math>
		\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
	</math>
y se deriva con respecto a <math>x_i</math>. Si hacemos que <math>\textbf{a}_i^\prime</math> sea la <math>i</math>-ésima fila de la matriz simétrica <math>\textbf{A}</math>, la derivada parcial de la forma cuadrática con respecto a <math>x_i</math> será
 <math>
		2\textbf{a}^\prime_i\textbf{x} = 2 \sum_{j=1}^p a_{ij}x_j
 </math>
y el vector de derivadas parciales es
 <math>
			\frac{\partial\textbf{x}^\prime\textbf{Ax}}{\partial\textbf{x}}= 2\textbf{Ax}
 </math>
En particular, el vector de derivadas de la suma de cuadrados <math>\textbf{x}^\prime\textbf{x}</math> es simplemente
<math>2\textbf{x}</math>.
# En la función cuadrática más general
 <math>
			h(\textbf{x}) = (\textbf{a}- \textbf{Cx})^\prime\textbf{K}(\textbf{a} - \textbf{Cx})
</math>
sea \textbf{K} una matriz simétrica de <math>N \times  N</math>, <math>\textbf{C} = [\textbf{C}_1, \ldots, \textbf{C}_p]</math> es una
matriz de constantes <math>N  \times  p</math>, y <math>\textbf{a}</math> es un vector de <math>N \times 1</math>. Si hacemos que <math>\textbf{u} =
		\textbf{a} - \textbf{Cx}</math>, las derivadas de <math>h(\textbf{x})</math> se pueden calcular por medio de la regla de la cadena:
 <math>
		\begin{align}
			\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{align}
 </math>
y
 <math>
			\frac{\partial h(\textbf{x})}{\partial\textbf{x}}= -2\textbf{C}^\prime\textbf{K}(\textbf{a} - \textbf{Cx)}
</math>
# La matriz de derivadas parciales de segundo orden de una función de <math>p</math> variables se conoce como la matriz ''hesiana'' :
 <math>
			\textbf{H} = \frac{\partial^2f(\textbf{x})}{\partial\textbf{x}^\prime\partial\textbf{x}} = 
			\begin{bmatrix}
				\frac{\partial^2f(\textbf{x})}{\partial x_i^2} & \cdots &
				\frac{\partial^2f(\textbf{x})}{\partial x_1\partial x_p}\\
				\cdots & \cdots & \cdots \\
				\frac{\partial^2f(\textbf{x})}{\partial x_1\partial x_p} & \cdots & \frac{\partial^2f(\textbf{x})}{\partial x_p^2}
			\end{bmatrix}
</math>
Por ejemplo, la matriz hesiana de <math>f(\textbf{x}) = \textbf{x}^\prime\textbf{Ax}</math> es simplemente <math>2\textbf{A}</math>. La
matriz hesiana es necesariamente simétrica si <math>f(\textbf{x})</math> 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 <math>f(\textbf{x})</math> en <math>\textbf{x}=\textbf{x}_0</math> es que 
 <math>
	\frac{\partial f(\textbf{x})}{\partial\textbf{x}} = \textbf{0}
 </math>
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 <math>\textbf{x}</math> 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 <math>\textbf{x}_0</math> donde se satisface la ecuación X es que la matriz hesiana evaluada en <math>\textbf{x}_0</math> 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 <math>f(\textbf{x})</math> en
serie de Taylor en la vecindad del punto estacionario. Si la hesiana es una matriz indefinida de <math>\textbf{x}_0</math>, 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 <math>h(\textbf{x})</math> 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 
 <math>
	\textbf{C}^\prime\textbf{KCx} = \textbf{C}^\prime\textbf{KA}
 </math>
y si <math>\textbf{C}^\prime\textbf{KC}</math> es de rango completo <math>p</math>, entonces
 <math>
	\textbf{x} = (\textbf{C}^\prime\textbf{KC})^{-1}\textbf{C}^\prime\textbf{Ka}
 </math>
Si <math>\textbf{K}</math> es positiva definida y <math>\textbf{C}</math> es de rango <math>p</math>, la hesiana <math>2\textbf{C}^\prime\textbf{KC}</math> será
positiva definida y la solución () minimizará <math>h(\textbf{x})</math>.

=== Maximización sujeta a restricciones ===
Con frecuencia será necesario maximizar o minimizar una función <math>f(\textbf{x})</math> sujeta a la restricción
<math>g(\textbf{x})=c</math> sobre los valores de <math>\textbf{x}</math>. 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
 <math>
	h(\textbf{x},\lambda) = f(\textbf{x}) - \lambda[g(\textbf{x)} - c]
 </math>
Para un valor estacionario restringido se deben cumplir estas condiciones:
 <math>
\begin{align}
	\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{align}
 </math>
La segunda condición es simplemente la restricción original. En la práctica uno debe resolver las ecuaciones
() con respecto a la variable <math>\textbf{x}</math> después de haber eliminado <math>\lambda</math> 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 <math>f(\textbf{x}) = \textbf{x}^\prime\textbf{Ax}</math> 
donde <math>\textbf{A}</math> es una matriz positiva definida, sujeta a la restricción
<math>\textbf{x}^\prime\textbf{x}=1</math>. Entonces
 <math>
h(\textbf{x},\lambda) = \textbf{x}^\prime\textbf{Ax} - \lambda(\textbf{x}^\prime\textbf{x} -1)
</math>
y () es el sistema de ecuaciones lineales
 <math>
	[\textbf{A} - \lambda\textbf{I}]\textbf{x}=\textbf{0}
 </math>
que no es nada más sino la definición de los valores característicos de <math>\textbf{A}</math>. Premultiplicando por
<math>\textbf{x}^\prime</math> y usando la restricción nos proporciona <math>\lambda = \textbf{x}^\prime\textbf{Ax}</math>. Obviamente, si la
forma cuadrática va a ser un máximo, entonces <math>\lambda</math> debe ser la raíz característica más grande de
<math>\textbf{A}</math>, y <math>\textbf{x}</math> 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 <math>f(\textbf{x})</math> 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 <math>\textbf{A}</math> de <math>n  \times n</math> con respecto a su elemento <math>a_{ij}</math> se puede
encontrar a partir de la expansión de <math>\textbf{A}</math> en los cofactores de la <math>i</math>-ésima fila o <math>j</math>-ésima columna. Usando
la primera expansión
 <math>
\begin{align}
	\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{align}
</math>
Si <math>\textbf{A}</math> es simétrica, 
 <math>
	\frac{\partial |\textbf{A}|}{\partial a_{ii}} = A_{ii}
 </math>
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
 <math>
	\frac{\partial| \textbf{A}|}{\partial a_{ij}} = A_{ij}
 </math>
Más adelante será necesario derivar formas cuadráticas y otros escalares con respecto a los elementos de sus matrices
constituyentes. Si <math>\textbf{X}</math> es una matriz de <math>m \times n</math> con un elemento general <math>x_{ij}</math>, su derivada con respecto a ese elemento es 
 <math>
	\frac{\partial\textbf{X}}{\partial x_{ij}} = \textbf{J}_{ij}
 </math>
donde <math>\textbf{J}_{ij}</math> es una matriz de <math>m  \times n</math> con uno en la posición <math>ij</math> y ceros en las demás. Si
<math>\textbf{X}</math> es simétrica
 <math>
	\frac{\partial\textbf{X}}{\partial x_{ij}}= J_{ij} + J_{ji} \;\;\; i \neq j
 </math>

La regla para derivar un producto de matrices es parecida a la de los escalares. Supóngase que los elementos de las
matrices conformables <math>\textbf{X}</math> y <math>\textbf{Y}</math> son funciones <math>x_{ij}(z), y_{ij}(z)</math>de alguna variable <math>z</math>. Entonces
 <math>
	\frac{\partial\textbf{XY}}{\partial z}=\frac{\partial\textbf{X}}{\partial x}\textbf{Y} +
	\textbf{X}\frac{\partial\textbf{Y}}{\partial z}
 </math>
Esta fórmula lleva a una forma de derivar la inversa de la matriz cuadrada no singular <math>\textbf{X}</math>. Expresamos
 <math>
\textbf{I} = \textbf{XX}^{-1}
</math>
Entonces
 <math>
\begin{align}
\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{align}
</math>
y
 <math>
	\frac{\partial\textbf{X}^{-1}}{\partial x_{ij}}= -\textbf{X}^{-1}\textbf{J}_{ij}\textbf{X}^{-1}
 </math>
Si <math>\textbf{X}</math> es simétrica
 <math>
	\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}
 </math>

== Avances ==
El siguiente capítulo tratará la implementación de rutinas para el manejo de matrices y vectores usando el lenguaje
<math>\textsf{R}</math>. Esto concluye el presente capítulo.

= Vectores y Matrices en R =
=== Vectores ===
Como ya hemos visto, <math>\textsf{R}</math> 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 <math>\textsf{R}</math>. Es necesario saber como crear, unir, y obtener subconjuntos de estos objetos. Dentro de
<math>\textsf{R}</math>, 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
<math>\texttt{NA}</math> que quiere decir ''Not Available'' (No se encuentra disponible). El tipo <math>\texttt{NA}</math> 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 <math>\texttt{vector()}</math> produce un vector de tipo <math>\texttt{logical}</math> por
omisión. Si el vector tiene una longitud, todos los elementos se inicializan con el valor de <math>\texttt{FALSE}</math>. 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 <math>\texttt{c()}</math> regresa con un vector.
Cuando se le asigna un objeto al resultado de <math>\texttt{c()}</math>, 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 <math>\texttt{c()}</math>,
<math>\textsf{R}</math> 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 <math>\texttt{x}</math>. En la línea 3
se crea un vector numérico de tipo <math>\texttt{double}</math> con los valores 1, 2, y 3. Al encadenar los vectores
<math>\texttt{x}</math> y <math>\texttt{y}</math>, <math>\textsf{R}</math> convierte todo a modo carácter. Por eso terminan encerrados con comillas. Esto
puede parecer confuso, pero al tomar en cuenta que <math>\textsf{R}</math> 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:

<math>\texttt{length()}</math>: Para obtener la longitud de un vector
 > (t.f <- c(TRUE, TRUE, FALSE, TRUE)); length(t.f)
 [1] TRUE TRUE FALSE TRUE
 [1] 4

<math>\texttt{sum()}</math>: 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 <math>\texttt{x}</math> son alterados individualmente por la operación especificada.

Cuando dos vectores son de la misma longitud y del mismo tipo, el resultado de <math>\texttt{x/y}</math> es un vector cuyo primer elemento es <math>\texttt{x[1]/y[1]}</math>, el segundo es <math>\texttt{x[2]/y[2]}</math>, 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
<math>\texttt{x > y}</math> y se usa aritmética binaria, <math>\textsf{R}</math> iterará sobre los elementos de <math>\texttt{y}</math> hasta que su longitud sea igual a la de <math>\texttt{y}</math>. Sin embargo, <math>\textsf{R}</math> 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, <math>\textsf{R}</math> procederá sin problemas. Por ejemplo, si se divide un vector <math>\texttt{x}</math> de longitud 10 entre un vector <math>\texttt{y}</math> de longitud 5, <math>\textsf{R}</math> usará los elementos del vector
<math>\texttt{y}</math> 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, <math>\textsf{R}</math> 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 <math>\texttt{x}</math> no es un múltiplo entero de <math>\texttt{y}</math>.

Esta misma regla se aplica con funciones que implementan aritmética vectorial. Por ejemplo, considérese la función de raíz cuadrada: <math>\texttt{sqrt()}</math>. <math>\textsf{R}</math> 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. <math>\textsf{R}</math> 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 <math>\texttt{s}</math> usamos la función
<math>\texttt{paste()}</math>:
 > paste(s[1], s[2])
 [1] "Our friend O'Brien the engineer"
Por omisión, <math>\texttt{paste()}</math> separa sus argumentos con un espacio. Si se desea un carácter diferente para espaciar
los elementos, entonces se usa el argumento <math>\texttt{sep}</math>:
 > paste(s[1], s[2], sep = '-')
 [1] "Our friend O'Brien-the engineer"
Si se desea más información sobre la función <math>\texttt{paste()}</math> úsese <math>\texttt{help(paste)}</math>.

==== 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 <math>\texttt{NA}</math>, el resultado va a ser <math>\texttt{NA}</math>. De modo que
necesitamos alguna forma de extraer del vector los datos que no son <math>\texttt{NA}</math>. Vamos a crear el siguiente vector:
 > (i <- is.na(x))
 [1] FALSE FALSE  TRUE FALSE  TRUE FALSE
La función <math>\texttt{is.na(x)}</math> examina cada elemento del vector <math>\texttt{x}</math> y si es <math>\texttt{NA}</math> le asigna al
correspondiente elemento del vector <math>\texttt{i}</math> el valor de <math>\texttt{TRUE}</math>; de lo contrario, le asigna el valor
<math>\texttt{FALSE}</math>. Ahora vamos a usar al vector <math>\texttt{i}</math> como un vector índice para extraer los valores que no son
<math>\texttt{NA}</math> del vector <math>\texttt{x}</math>. Tomando en cuenta que si <math>\texttt{x}</math> tiene un valor numérico, el elemento
correspondiente del vector <math>\texttt{i}</math> tiene el valor <math>\texttt{FALSE}</math>. Si <math>\texttt{x}</math> no tiene valor, entonces el elemento
correspondiente de <math>\texttt{i}</math> es <math>\texttt{TRUE}</math>. Usamos entonces <math>\texttt{i}</math> y la negación de <math>\texttt{i}</math> 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 <math>\texttt{x}</math>, 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 <math>\texttt{mean}</math> es una función intrínseca de <math>\textsf{R}</math> que calcula el promedio de un vector. El primer
argumento de <math>\texttt{mean}</math> es el vector de datos, y el segundo es un argumento lógico que si se le asigna el valor de
<math>\texttt{TRUE}</math> calcula el promedio sin tomar en cuenta los valores <math>\texttt{NA}</math>. El valor de <math>\texttt{na.rm}</math> por omisión
es <math>\texttt{FALSE}</math>.

'''Ejemplo: Vector índice de enteros positivos.''' Esta es otra forma para excluír datos tipo <math>\texttt{NA}</math>:
 > (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 <math>\texttt{x}</math> 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 <math>\texttt{20:30}</math> crea un vector de enteros que corren del 20 al 30 y se usa para extraer los
correspondientes elementos de <math>\texttt{x}</math>, mismos que se redondean a tres decimales usando la función <math>\texttt{round}</math>.
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 <math>\texttt{x}</math> <math>\textbf{excepto}</math> los
elementos del 20 al 30, se expresará como: <math>\texttt{x[-(20:30)]}</math>. Obsérvese como encerramos a <math>\texttt{20:30}</math> entre
paréntesis, ya que si no lo hubiéramos hecho, <math>\textsf{R}</math> 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 <math>\texttt{x}</math> y la altura
<math>\texttt{y}</math> 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 <math>\texttt{names(x)}</math> les da un nombre a los elementos de <math>\texttt{x}</math>. Para ver los datos los asociamos por
columna con la función <math>\texttt{cbind()}</math>:
 > 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: <math>\texttt{v[i]}</math>. Las
matrices son arreglos rectangulares en dos dimensiones. Los elementos en las matrices se accede a ellos por medio de la
notación <math>\texttt{m[i,j]}</math>. Los arreglos tienen <math>k</math> dimensiones. Cada elemento de un arreglo es accedido por medio de <math>k</math>
índices: <math>\texttt{a[i1, i2, \dots, ik]}</math>.

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 <math>\texttt{dim}</math>. Observemos aquí algunas
instrucciones en <math>\textsf{R}</math> 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 <math>\texttt{dim}</math>. 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 <math>\texttt{matrix()}</math>:
 > (m <- matrix(0, ncol=3, nrow=2))
 [,1] [,2] [,3]
 [1,]    0    0    0
 [2,]    0    0    0
Ahora, verificamos si <math>\texttt{m}</math> es un arreglo con la función <math>\texttt{is.array()}</math>:
 > c('matrix?' = is.matrix(m), árray?' = is.array(m))
 matrix?  array?
 TRUE    TRUE
Como nos muestra la función <math>\texttt{is.matrix()}</math>, <math>\texttt{m}</math> 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 <math>\texttt{byrow}</math> en la función <math>\texttt{matrix()}</math>.

=== Problemas Resueltos ===
# Sea
 <math>\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}
 </math>
Realícense las siguientes operaciones que se encuentren definidas:
## <math>\textbf{A + B}</math>
## <math>\textbf{A} - 2\textbf{B}</math>
## <math>\textbf{A}^\prime + \textbf{B}</math>
## <math>\textbf{A}+\textbf{C}</math>
## (<math>\textbf{A + B)}^\prime</math>
## (<math>3\textbf{A}^\prime - 2\textbf{B})^\prime</math>
 > (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

# Si
 <math>\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}
 </math>
 <math>\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}
 </math>
Calcúlense aquellos productos matriciales que se encuentren definidos:
## <math>\textbf{PQ}</math>
## <math>\textbf{PQR}</math>
## <math>\textbf{QR}^\prime</math>
## <math>\textbf{yx}^\prime</math>
## <math>\textbf{x}^\prime\textbf{y}</math>
## <math>\textbf{x}^\prime\textbf{Px}</math>
## <math>\textbf{x}^\prime\textbf{Py}</math>
## <math>\textbf{P}(\textbf{x+y})</math>
 > (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