potencianormal

La distribución potencia normal en R

(The Power normal distribution with R)

Considere la familia de distribuciones, propuesta por Lehmann (1953) y extendida por Durrans (1992), con función de distribución acumulada y función de densidad de probabilidad dadas, respectivamente, por:

donde F es en sí una función de distribución.

Una extensión, de localización escala, de esta familia de distribuciones con distribución generadora normal se encuentra en Pewsey, Gómez y Bolfarine (2012) y se denomina la distribución potencia normal con función de densidad de probabilidad dada por:

A continuación se encuentran funciones R para la fdp, la fda, la función cuantil y la generación aleatoria de la distribución potencia normal con parametros location, scale y shape.

dpn = function(x,location=0,scale=1,shape){

if(shape<=0 | scale<=0) stop("the shape and scale parameters must be greater than zero")

if(shape==1) d = dnorm(x,location,scale)

if(shape!=1) d = shape*dnorm((x-location)/scale)*pnorm((x-location)/scale)^(shape-1)/scale

return(d)

}

pap = function(x,location=0,scale=1,shape){

if(shape<=0 | scale<=0) stop("the shape and scale parameters must be greater than zero")

if(shape==1) p = pnorm((x-location)/scale)

if(shape!=1) p = pnorm((x-location)/scale)^shape

return(p)

}

qpn = function(p,location=0,scale=1,shape){

if(shape<=0 | scale<=0) stop("the shape and scale parameters must be greater than zero")

if(shape==1) q = qnorm(p,location,scale)

if(shape!=1) q = location+scale*qnorm(p^(1/shape))

return(q)

}

rpn = function(n,location=0,scale=1,shape){

if(shape<=0 | scale<=0) stop("the shape and scale parameters must be greater than zero")

if(shape==1) r = rnorm(n,location,scale)

if(shape!=1) r = location+scale*qnorm(runif(n)^(1/shape))

return(r)

}

Función de log-verosimilitud, score y de estimación de los parámetros:

lpn = function(theta,x){

require(sn)

location = theta[1] ; scale = theta[2] ; shape = theta[3]

n = length(x) ; z = (x-location)/scale

n*(log(shape)-log(scale)-0.5*log(2*pi))-sum((z^2)/2)+(shape-1)*sum(zeta(0,z)-log(2))

}

score = function(theta,x){

require(sn)

location = theta[1] ; scale = theta[2] ; shape = theta[3] ; n = length(x)

z = (x-location)/scale ; w = zeta(1,z) ; u = zeta(0,z)-log(2)

c(n*(mean(z)-(shape-1)*mean(w))/scale,-n*(1-mean(z^2)+(shape-1)*mean(z*w))/scale,n*(mean(u)+1/shape))

}

fitpn = function(x){

z=(x-mean(x))/sd(x) ; u=zeta(0,z)-log(2)

optim(c(mean(x),sd(x),-1/mean(u)),lpn,score,x=x,method="L-BFGS-B",

lower=c(-Inf,1e-22,1e-22),upper=rep(Inf,3),control=list(fnscale=-1))$par

}

Ejemplo:

data(ais,package="sn")

x = ais$BMI

param = fitpn(x)

hist(x,freq=F)

lines(d<-density(x))

lines(d$x,dpn(d$x,param[1],param[2],param[3]),col=4)

Referencias:

- Durrans SR (1992) Distributions of fractional order statistics in hydrology. Water Resour Res 28:1649–1655.

- Lehmann EL (1953) The power of rank tests. Ann Math Stat 24:23–43.

- Pewsey, A., Gómez, HW y Bolfarine, H. (2012) Likelihood-based inference for power distributions. Test 21:775-789.

Comparte esto: