Box Cox

Descargar como ppt, pdf o txt
Descargar como ppt, pdf o txt
Está en la página 1de 14

Transformacin de Potencia Box-Cox

Modelos Estadsticos
Dra. Graciela Gonzlez Faras Jos Ramn Domnguez Molina 14/marzo/2003 Omar Posada Villarreal

Transformacin de potencia
s s

Simple Se requiere que la distribucin sea


Suave Continua X>0

Transformacin de potencia
X 1 0 T(X ) = Y = ln X =0

s s s

=2, Y=X2 =1/2, Y=X1/2 Se busca que la variable transformada se parezca a una distribucin normal

YX

( )

~ N ( , )
2

Ejemplo: X ~ Exp(1)
s s

Rango: [-2, 2] pasos de 0.5. La mejor fue = 0.5


Exponencial( 1 ) QQPlot normalizado. Desv. Tipica = 1.16
6

Tran. Box-Cox con SD min. QQPlot norm. (lambda = 0.5 , Desv. Tip. = 0.988 )
6 Y

-2

-2

-1

-2

-2

-1

Quantiles of Standard Normal

Quantiles of Standard Normal

Exponencial( 1 ) Histograma
60

Lambda vs. Desv. Tipica. (lambda = 0.5 , Desv. Tip. = 0.988 )

40

log(SD) 0 2 X 4 6

20

1 -2

5 10

50

-1

0 Lambda

Ejemplo: X ~ Exp(3)
s s

Rango: [0, 10] pasos de 0.05. La mejor fue = 3.05


Exponencial( 3 ) QQPlot normalizado. Desv. Tipica = 0.364
1.5

Tran. Box-Cox con SD min. QQPlot norm. (lambda = 3.05 , Desv. Tip. = 0.221 )
1.5 Y

1.0

0.5

0.0

-2

-1

0.0

0.5

1.0

-2

-1

Quantiles of Standard Normal

Quantiles of Standard Normal

Exponencial( 3 ) Histograma
40 1.0 log(SD) 0.0 0.5 X 1.0 1.5 0.3 0.5 0.7

Lambda vs. Desv. Tipica. (lambda = 3.05 , Desv. Tip. = 0.221 )

10

20

30

4 Lambda

10

Ejemplo: X ~ U(0.01, 1)
s s

Rango: [-10, 10] pasos de 0.5. La mejor fue >= 10


Uniforme( 0.01 , 1 ) QQPlot normalizado. Desv. Tipica = 0.286
1.0 0.8

Tran. Box-Cox con SD min. QQPlot norm. (lambda = 10 , Desv. Tip. = 0.0222 )
0.8 Y 1.0

0.6

0.4

0.2

0.0

-2

-1

0.0

0.2

0.4

0.6

-2

-1

Quantiles of Standard Normal

Quantiles of Standard Normal

Uniforme( 0.01 , 1 ) Histograma


10 12 14 10^7 10^10 log(SD) 0.0 0.2 0.4 X 0.6 0.8 1.0 10^-2 10^1 10^4

Lambda vs. Desv. Tipica. (lambda = 10 , Desv. Tip. = 0.0222 )

-10

-5

0 Lambda

10

Ejemplo: X ~ U(1, 5)
s s

Rango: [-10, 10] pasos de 1. La mejor fue <= -10


Uniforme( 1 , 5 ) QQPlot normalizado. Desv. Tipica = 1.23
5 4

Tran. Box-Cox con SD min. QQPlot norm. (lambda = -10 , Desv. Tip. = 0.00843 )
4 Y 5

-2

-1

-2

-1

Quantiles of Standard Normal

Quantiles of Standard Normal

Uniforme( 1 , 5 ) Histograma
10^4 log(SD) 1 2 3 X 4 5 10^-2 10^0 10^2 15

Lambda vs. Desv. Tipica. (lambda = -10 , Desv. Tip. = 0.00843 )

10

-10

-5

0 Lambda

10

Ejemplo: X ~ Beta(5, 2.5)


s s

Rango: [-10, 10] pasos de 1. La mejor fue >= 10


Beta( 5 , 2.5 ) QQPlot normalizado. Desv. Tipica = 0.165 Tran. Box-Cox con SD min. QQPlot norm. (lambda = 10 , Desv. Tip. = 0.0138 )

1.0

0.8

0.6

0.4

Y -2 -1 0 1 2

0.2

0.0

0.0

0.2

0.4

0.6

0.8

1.0

-2

-1

Quantiles of Standard Normal

Quantiles of Standard Normal

25

10

log(SD)

0.2

0.4 X

0.6

0.8

1.0

10^-2 -10

10^1

10^3

15

10^5

20

10^7

Beta( 5 , 2.5 ) Histograma

Lambda vs. Desv. Tipica. (lambda = 10 , Desv. Tip. = 0.0138 )

-5

0 Lambda

10

Listado S-Plus (1)


s s s s s s s s s

# Realiza una transformacin que se ajuste a la normal # @param fX Datos # @param leftlambda Limite inferior para probar lambda # @param rightLambda Limite superior para probar lambda # @param eachLambda Intervalo entre marcas boxCox = function(fX, leftLambda, rightLambda, eachLambda) { cX = data.matrix(fX) dimX = dim(cX) n = dimX[1] origSD = stdev(cX) #Equivale a cXLambda1 = (cX ^ 1) - 1 #origSD = stdev(cXLambda1) # Checa que xi>0 for (i in 1:n) { if (cX[i] <= 0) { stop("Debe ser: x[i]>0") } } # Inicializar # Rango de lambdas a probar minLambda = rightLambda rLambda = seq(leftLambda, rightLambda, by=eachLambda) nLambda = length(rLambda) minSD = 1E100 rSD = vector(mode="numeric", length=nLambda) cY = vector(mode="numeric", length=n)

s s s

s s s s s s s s s s s s s s

Listado S-Plus (2)


s s s s s s s s s s s s s s s s s s s s s

# Para cada lambda for (j in 1:nLambda) { # Transformacion Box-Cox # print(paste("- i=", i , " j=", j)) if (rLambda[j] != 0) { cY = (cX ^ rLambda[j] - 1) / rLambda[j] } else { cY = log(cX) } # Recuerda el vector con min stdev rSD[j] = stdev(cY) if (rSD[j] < minSD) { cMinY = cY minLambda = rLambda[j] minSD = rSD[j] } } return (cX, origSD, rLambda, rSD, cMinY, minLambda, minSD) }

Listado S-Plus (3)


s s s s s s s s s

plotBoxCox = function(sTitle, cX, origSD, rLambda, rSD, cMinY, minLambda, minSD) {

print("Original")
# En una pagina par(mfrow = c(2,2)) options(digits=3) # Conserva la mayor escala de los datos orig y tran en el eje Y minY = min(cX, cMinY) maxY = max(cX, cMinY) # Grafica qqplot normalizado de los datos originales # Muestra la varianza actual. sTitle2 = paste(sTitle, "\nQQPlot normalizado. Desv. Tipica = ", format(origSD)) qqnorm(cX, main=sTitle2, ylab="X", ylim=c(minY, maxY)) qqline(cX)

s s s s s

s s s

print("Transformada")
# Grafica transformacion con Desv. Tipica sTitle2 = paste("Tran. Box-Cox con SD min. QQPlot norm.\n(lambda = ", format(minLambda), ", Desv. Tip. = ", format(minSD), ")") qqnorm(cMinY, main=sTitle2, ylab="Y", ylim=c(minY, maxY)) qqline(cMinY)

s s

s s s s s s s

print("Histograma")
sTitle2 = paste(sTitle, "\nHistograma") hist(cX, main=sTitle2, xlab="X")

print("Lambda")
# Grafica lambda vs. Desv. Tipica sTitle2 = paste("Lambda vs. Desv. Tipica.\n(lambda = ", format(minLambda), ", Desv. Tip. = ", format(minSD), ")") plot(rLambda, rSD, main=sTitle2, xlab="Lambda", ylab="log(SD)", log='y')

Listado S-Plus (4)


s s s s

# PARAMETROS DEL PROGRAMA # Inicializar archivo example = 5 n = 100 # Tamano de muestra # Parametros de los ejemplos # El dominio debe ser X>0 if (example == 1) { print("Exp") lambda1 = 1 # Parmetro para exp sTitle = paste("Exponencial(", lambda1, ")") leftLambda = -2 rightLambda = 2 eachLambda = 0.05 cXOrig = rexp(n, lambda1) } else if (example == 2) { print("Exp") lambda1 = 3 # Parmetro para exp sTitle = paste("Exponencial(", lambda1, ")") leftLambda = 0 rightLambda = 10 eachLambda = 0.05 cXOrig = rexp(n, lambda1)

s s s s s s

s s s s

s s s s

s s s s

Listado S-Plus (5)


s s s s s s s s s s s s s s s s s s s s s s s s s

} else if (example == 3) { print("Unif") alfa = 0.01 # Parmetro para Unif beta = 1 # Parmetro para Unif sTitle = paste("Uniforme(", alfa, ", ", beta, ")") leftLambda = -10 rightLambda = 10 eachLambda = 0.5 cXOrig = runif(n, min=alfa, } else if (example == 4) { print("Unif") alfa = 1 # Parmetro para beta = 5 # Parmetro para sTitle = paste("Uniforme(",

max=beta)

Unif Unif alfa, ", ", beta, ")")

leftLambda = -10 rightLambda = 10 eachLambda = 1 cXOrig = runif(n, min=alfa, max=beta) } else if (example == 5) { print("Beta") alfa = 5 # Parmetro para Unif beta = 2.5 # Parmetro para Unif sTitle = paste("Beta(", alfa, ", ", beta, ")") leftLambda = -10 rightLambda = 10 eachLambda = 1 cXOrig = rbeta(n, alfa, beta)

s s s s s

Listado S-Plus (6)


s s s s

# Escribe en archivo una muestra aleatoria con distribucion exponencial cXOrig = t(cXOrig) cXOrig = t(cXOrig) # Dos veces para transponer renglon a columna (?) exportData(cXOrig, "D:\\Posada\\ModEst\\ModEst4\\ExpSample.txt", type="ASCII") fX = importData("D:\\Posada\\ModEst\\ModEst4\\ExpSample.txt", type="ASCII") res = boxCox(fX, leftLambda, rightLambda, eachLambda) plotBoxCox(sTitle, res$cX, res$origSD, res$rLambda, res$rSD, res$cMinY, res$minLambda, res$minSD)

s s s

También podría gustarte