Box Cox
Box Cox
Box Cox
Modelos Estadsticos
Dra. Graciela Gonzlez Faras Jos Ramn Domnguez Molina 14/marzo/2003 Omar Posada Villarreal
Transformacin de potencia
s s
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
Tran. Box-Cox con SD min. QQPlot norm. (lambda = 0.5 , Desv. Tip. = 0.988 )
6 Y
-2
-2
-1
-2
-2
-1
Exponencial( 1 ) Histograma
60
40
log(SD) 0 2 X 4 6
20
1 -2
5 10
50
-1
0 Lambda
Ejemplo: X ~ Exp(3)
s s
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
Exponencial( 3 ) Histograma
40 1.0 log(SD) 0.0 0.5 X 1.0 1.5 0.3 0.5 0.7
10
20
30
4 Lambda
10
Ejemplo: X ~ U(0.01, 1)
s s
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
-10
-5
0 Lambda
10
Ejemplo: X ~ U(1, 5)
s s
Tran. Box-Cox con SD min. QQPlot norm. (lambda = -10 , Desv. Tip. = 0.00843 )
4 Y 5
-2
-1
-2
-1
Uniforme( 1 , 5 ) Histograma
10^4 log(SD) 1 2 3 X 4 5 10^-2 10^0 10^2 15
10
-10
-5
0 Lambda
10
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
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
-5
0 Lambda
10
# 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
# 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) }
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')
# 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
} 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)
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
# 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