R help es - Aug 2012 - Prueba de Thom y series temporales homogéneas de prueba

Hola a todos:

Hace unos meses estuve preguntando sobre la prueba de Thom para 
comprobar la homogeneidad de series temporales sobre datos de radiación 
e insolación (después he visto que también se usa con otras variables 
meteorológicas) [0].

El caso es que por distintas razones no había podido avanzar en el 
asunto y ahora vuelvo sobre él.

Creé una función para realizar la prueba de Thom pero el caso es que 
series temporales de radiación que «deberían» pasar la prueba no lo 
hacen. (Digo deberían porque he encontrado múltiples artículos donde 
mencionan los mismos datos que yo uso (AEMET) y pasan la prueba).

Como de alguna forma quiero comprobar que la función esta bien 
construida he intentado buscar alguna documentación para crear series 
temporales homogéneas pero no he conseguido dar con nada para hacerlo 
(puede deberse también a mi torpeza, claro).

Expongo el código a continuación:

--------------- código R -------------
contar.pasos<-function(vector,valor.cruce) {
     # Función que cuenta los cruces de un vector por un
     # valor de referencia (valor.cruce)

     paso<-0
     i<-1

     while ( i <= (length(vector)-1) ) {

       cruce.corto<-(vector[i] < valor.cruce & vector[(i+1)] > 
valor.cruce ) |
	  (vector[i] > valor.cruce & vector[(i+1)] < valor.cruce )

       if ( vector[i] == valor.cruce & i!=1 ) {
       # La función esta en el valor de cruce y no es el primer valor,
       # se presentan varios casos

	baja<-( (vector[i-1] > vector[i]) & (vector[i] > vector[i+1]) )
	sube<-( (vector[i-1] < vector[i]) & (vector[i] < vector[i+1]) )

	if (sube | baja ) {pcaso<-TRUE} else {pcaso<-FALSE}
	# La función cruza el valor de referencia.

	if (vector[(i-1)] == vector[i]) {
	  # El valor anterior de la función también esta en el valor de cruce,
	  # entonces comprueba como era la función después del tramo horizontal
	  # y como evoluciona después
	  z<-i-1
	  while (vector[z] == vector[i] & z>1) {z<-z-1}
	  # Retrocedemos hasta el último valor distinto sin pasar
	  # del primer valor de la función.
	  menor<-FALSE
	  mayor<-FALSE
	  if (vector[z] < vector[i]) { menor<-TRUE } else {mayor<-TRUE}
	  if (
	    ( menor & (vector[i] < vector[(i+1)]) ) |
	    ( mayor & (vector[i] > vector[(i+1)]) )
	  ) {scaso<-TRUE}
	} else {scaso<-FALSE}
	
	cruce.largo<-(pcaso | scaso)

       } else {cruce.largo<-FALSE}

       excepcion<-if (
	i==1 | (vector[i] == valor.cruce & !cruce.largo)
	) {TRUE} else {FALSE}
	
       if ( (cruce.corto | cruce.largo) & !excepcion ) {paso<-paso+1}
	# Si la función cruza y no hay excepción

       i<-i+1
     }

     return(paso)
}

pdThom<-function(vector,na.rm=F) {
     # Prueba de Thom, de las rachas o de las alternancias. (1.966)
     # Prueba de homogeneidad absoluta aconsejada para distribuciones no
     # conocidas o no normales.
     # Depende de la función contar.pasos

     longitud<-length(vector)
     mediana<-median(vector,na.rm=na.rm)
     rachas<-contar.pasos(vector,mediana)

     resultado<-(rachas-((longitud+2)/2))/
       sqrt((longitud*(longitud-2))/(4*(longitud-1)))

     return(abs(resultado))
}

library(solaR)

SIAR <- read.csv('http://solar.r-forge.r-project.org/data/SIAR.csv')
#
Estaciones.Almeria<-subset(SIAR,Provincia == "Almería")
#
# Datos estación de La Mojonera
# Almería, provincia = 4, estación = 1
Estacion.La.Mojonera<-readSIAR(4,1,"01/01/1970","31/12/2011")
#
# Datos estación de Adra
# Estación = 10
Estacion.Adra<-readSIAR(4,10,"01/01/1970","31/12/2011")
#
# Datos estación Almería
# Estación = 2
Estacion.Almeria<-readSIAR(4,2,"01/01/1970","31/12/2011")
#
# Radiación diaria en Adra.
Adra.radiacion.diaria<-getG0(Estacion.Adra)

# Radiación diaria en La Mojonera
Mojonera.radiacion.diaria<-getG0(Estacion.La.Mojonera)

# Radiación diaria en Almería
Almeria.radiacion.diaria<-getG0(Estacion.Almeria)

# Temperatura diaria en Almeria
Almeria.temp<-Estacion.Almeria en data$TempMedia

## Pruebas de homogeneidad
which(is.na(as.vector(Almeria.temp)))
# Resultado: integer(0), no hay valores perdidos

pdThom(as.vector(Almeria.temp))
# Resultado: [1] 59.20887, bastante alejado de 2,58

which(is.na(as.vector(Adra.radiacion.diaria)))
#  [1]    1   11  825 1435 1449 1450 1451 1452 1453 1454 1756 1868 1895 
# 1896 1897
# [16] 1903 1947 1948 1949 1950 1951 1952 1954 3465 3614

pdThom(as.vector(Adra.radiacion.diaria)[12:824])
# [1] 22.63508

which(is.na(as.vector(Almeria.radiacion.diaria)))
# [1]  886 1838 3221 3222 3223 3224 3227 3402 3470 3592 3693

pdThom(as.vector(Almeria.radiacion.diaria)[1839:3220])
# [1] 28.201

library(TSA)
data(tempdub)
pdThom(tempdub)
# [1] 8.195372
data(larain)
pdThom(larain)
# [1] 0.4683109
------------------- Fin del código R -----------

Como veis la única serie que «pasa» la prueba de Thom es la de lluvia en 
Los Angeles del paquete TSA, pero estoy desorientado, no me lo creo.

¿Que opináis?

Gracias por adelantado.

Un saludo.

R. Gómez

[0] https://stat.ethz.ch/pipermail/r-help-es/2012-February/003236.html
-- 
Libertad es poder elegir en cualquier momento. Ahora yo elijo GNU/Linux,
para no atar mis manos con las cadenas del soft propietario.
---------
Desde El Ejido, en Almería, usuario registrado Linux #294013
http://www.counter.li.org

Buenos días.

En lo que se refiere a la función runs.test del paquete tserie es la misma
prueba.
En el help de runs.test hace referencia a

S. Siegel and N. J. Castellan (1988): *Nonparametric Statistics for the
Behavioural Sciences*, 2nd edn, McGraw-Hill, New York.

en su página 62 da la fórmula que emplea par el caso de muestras grandes.

mu_R = (2mn/N)+1 = {bajo la hipotesis nula , m = n = N/2} = (N+2)/2

sigma_R = (2mn(2mn-N))/(N^2 (N-1) =  {bajo la hipotesis nula , m = n = N/2}
= N(N-2) / 4(N-1)

que son las expresiones para la media y la varianza que se dan en el pdf,

http://www.ub.edu/gc/**Documentos/Sevilla_SunDu.pdf<http://www.ub.edu/gc/Documentos/Sevilla_SunDu.pdf>


Respecto a la otra función runs.test, del paquete lawstat, no la he mirado,
pero si te da los mismos resultados parece que es otra implementación del
mismo test.

Puedes encontrar  en,

http://bellman.ciencias.uniovi.es/estadistica2/estadistica2_archivos/aleatoriedad.pdf

una explicación más detallada del test.

Por último, efectivamente la eliminación de la copia a la lista fue
accidental.

Un saludo


El 20 de agosto de 2012 01:37, Rubén Gómez Antolí
<lobo@mucharuina.com>escribió:
> Hola José Luis:
>
> El 19/08/12 22:24, Jose Luis Brita escribió:
>
>  Buenas noches.
>>
>> Respecto a la prueba de Thom, si no estoy confundido, es el clásico
test
>> de rachas.
>>
>
> No soy estadístico y no estoy seguro.
>
> En la página 2 de este artículo esta descrita la prueba de Thom:
>
>
http://www.ub.edu/gc/**Documentos/Sevilla_SunDu.pdf<http://www.ub.edu/gc/Documentos/Sevilla_SunDu.pdf>
>
> que, me reitero, no se si es exactamente este test de rachas.
>
>
>  Si es así, puedes utilizar la función runs.test del paquete
>> tseries para chequear tu función. Un ejemplo de como se usa la función
>> runs.test lo tienes en la página 99 del libro Estadística Básica con R
y
>> R–Commander, de la Universidad de Cadiz que es libre.
>>
>
> Gracias por el aporte. Carlos Ortega en su día también me recomendó «otra»
> función runs.test, del paquete lawstat.
>
> Como te decía, ante la duda de si me valía o no, opté por implementar la
> prueba de Thom ya que, la formula en si, no entrañaba dificultad (la
> dificultad se encontraba en contar las rachas -función contar.pasos-).
>
> En concreto, con las dos pruebas runs.test obtengo el mismo resultado
> aunque su utilización es diferente:
>
> ---------------- consola de R ------------
>
> >
tseries::runs.test(as.factor(**Almeria.temp>median(Almeria.**temp)))
>
>         Runs Test
>
> data:  as.factor(Almeria.temp > median(Almeria.temp))
> Standard Normal = -59.178, p-value < 2.2e-16
> alternative hypothesis: two.sided
>
> > lawstat::runs.test(as.vector(**Almeria.temp))
>
>         Runs Test - Two sided
>
> data:  as.vector(Almeria.temp)
> Standardized Runs Statistic = -59.178, p-value < 2.2e-16
>
> ----------- fin de consola --------
>
> Por otro lado, voy a tener que repasar mi estadística para ver si este
> resultado valida mis datos o no.
>
>
>  Un saludo
>> --
>> Jose Luis
>>
>
> Gracias por tu interés.
>
> Salud y Revolución.
>
> Lobo.
>
> Pd: Curioso, mi prueba de Thom me esta devolviendo el valor «standardized
> runs statistic»:
>
> > pdThom(as.vector(Almeria.temp)**)
> [1] 59.20887
>
> > runs.test(tempdub)
>
>         Runs Test - Two sided
>
> data:  tempdub
> Standardized Runs Statistic = -8.0281, p-value = 9.898e-16
>
> > pdThom(tempdub)
> [1] 8.195372
>
>
> así que, voy a reconocerlo, no tengo ni idea de lo que estoy haciendo.
> Ahora si que voy a ponerme a repasar estadística. :^(
>
> Pd2: la eliminación de la copia a la lista, ¿es accidental?
>
> --
> Libertad es poder elegir en cualquier momento. Ahora yo elijo GNU/Linux,
> para no atar mis manos con las cadenas del soft propietario.
> ---------
> Desde El Ejido, en Almería, usuario registrado Linux #294013
> http://www.counter.li.org
>


-- 
Jose Luis

	[[alternative HTML version deleted]]

Hola:

El 20/08/12 09:55, Jose Luis Brita escribió:
> Buenos días.
>
> En lo que se refiere a la función runs.test del paquete tserie es la
> misma prueba.
> En el help de runs.test hace referencia a
>
> S. Siegel and N. J. Castellan (1988): /Nonparametric Statistics for the
> Behavioural Sciences/, 2nd edn, McGraw-Hill, New York.
>
> en su página 62 da la fórmula que emplea par el caso de muestras grandes.
>
> mu_R = (2mn/N)+1 = {bajo la hipotesis nula , m = n = N/2} = (N+2)/2
>
> sigma_R = (2mn(2mn-N))/(N^2 (N-1) = {bajo la hipotesis nula , m = n >
N/2} = N(N-2) / 4(N-1)
>
> que son las expresiones para la media y la varianza que se dan en el pdf,
>
> http://www.ub.edu/gc/__Documentos/Sevilla_SunDu.pdf
> <http://www.ub.edu/gc/Documentos/Sevilla_SunDu.pdf>
>
>
> Respecto a la otra función runs.test, del paquete lawstat, no la he
> mirado, pero si te da los mismos resultados parece que es otra
> implementación del mismo test.
Si, tiene toda la pinta, aunque a una se le pasan factores y a la otra 
vectores.
> Puedes encontrar  en,
>
>
http://bellman.ciencias.uniovi.es/estadistica2/estadistica2_archivos/aleatoriedad.pdf
>
>
> una explicación más detallada del test.
Gracias por ambas referencias, me pondré con ellas a ver si me aclaro.
> Por último, efectivamente la eliminación de la copia a la lista fue
> accidental.
Eso me pareció.
> Un saludo
Gracias por todo.

Salud y Revolución.

Lobo.
-- 
Libertad es poder elegir en cualquier momento. Ahora yo elijo GNU/Linux,
para no atar mis manos con las cadenas del soft propietario.
---------
Desde El Ejido, en Almería, usuario registrado Linux #294013
http://www.counter.li.org