180 likes | 464 Views
Transformación de Potencia Box-Cox. Modelos Estadísticos Dra. Graciela González Farías José Ramón Domínguez Molina 14/marzo/2003 Omar Posada Villarreal. Transformación de potencia. Simple Se requiere que la distribución sea Suave Continua X>0. Transformación de potencia. =2, Y=X 2
E N D
Transformación de Potencia Box-Cox Modelos Estadísticos Dra. Graciela González Farías José Ramón Domínguez Molina 14/marzo/2003 Omar Posada Villarreal
Transformación de potencia • Simple • Se requiere que la distribución sea • Suave • Continua • X>0
Transformación de potencia • =2, Y=X2 • =1/2, Y=X1/2 • Se busca que la variable transformada se parezca a una distribución normal
Ejemplo: X ~ Exp(1) • Rango: [-2, 2] pasos de 0.5. • La mejor fue = 0.5
Ejemplo: X ~ Exp(3) • Rango: [0, 10] pasos de 0.05. • La mejor fue = 3.05
Ejemplo: X ~ U(0.01, 1) • Rango: [-10, 10] pasos de 0.5. • La mejor fue >= 10
Ejemplo: X ~ U(1, 5) • Rango: [-10, 10] pasos de 1. • La mejor fue <= -10
Ejemplo: X ~ Beta(5, 2.5) • Rango: [-10, 10] pasos de 1. • La mejor fue >= 10
Listado S-Plus (1) • # Realiza una transformación 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)
Listado S-Plus (2) • # 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) • 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) • 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) • 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) • # 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 # Parámetro 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 # Parámetro para exp • sTitle = paste("Exponencial(", lambda1, ")") • leftLambda = 0 • rightLambda = 10 • eachLambda = 0.05 • cXOrig = rexp(n, lambda1)
Listado S-Plus (5) • } else if (example == 3) { • print("Unif") • alfa = 0.01 # Parámetro para Unif • beta = 1 # Parámetro para Unif • sTitle = paste("Uniforme(", alfa, ", ", beta, ")") • leftLambda = -10 • rightLambda = 10 • eachLambda = 0.5 • cXOrig = runif(n, min=alfa, max=beta) • } else if (example == 4) { • print("Unif") • alfa = 1 # Parámetro para Unif • beta = 5 # Parámetro para Unif • sTitle = paste("Uniforme(", alfa, ", ", beta, ")") • leftLambda = -10 • rightLambda = 10 • eachLambda = 1 • cXOrig = runif(n, min=alfa, max=beta) • } else if (example == 5) { • print("Beta") • alfa = 5 # Parámetro para Unif • beta = 2.5 # Parámetro para Unif • sTitle = paste("Beta(", alfa, ", ", beta, ")") • leftLambda = -10 • rightLambda = 10 • eachLambda = 1 • cXOrig = rbeta(n, alfa, beta) • } • # 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)
Listado S-Plus (6) • # 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)