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: