Ajusta una Pareto a los siguientes datos

library(ggplot2)
library(actuar)
datos=c(10, 11, 15, 22, 28, 30, 32, 36, 38, 48,
        55, 56, 68, 68, 85, 87, 94, 103, 104, 105,
        109, 119, 121, 137, 178, 181, 226, 287, 310, 106,
        393, 438, 591, 1045, 1210, 1212, 2423, 321, 354, 51)

Reconocimiento de la familia.

Esto se realiza mediante estadística descriptiva: cálculo de medidas descriptivas y gráficos estadísticos

hist(datos,freq=F,breaks = 20,col="blue")#¿de qué le veo cara?
d<-density(datos)
lines(d)

a= boxplot(datos)

Tenemos 3 picos, dos datos atipicos 1210, 1212, 2423, Datos sesgados a la derecha

Medidas de tendencia centrales

mean(datos)
[1] 272.675
median(datos) 
[1] 104.5
var(datos)
[1] 212649
summary(datos)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  10.00   50.25  104.50  272.68  292.75 2423.00 
quantile(datos)
     0%     25%     50%     75%    100% 
  10.00   50.25  104.50  292.75 2423.00 

Tenemos indicios de asimetría

Medidas de dispersión

var(datos)
[1] 212649
sd(datos)
[1] 461.1389
range(datos)
[1]   10 2423
IQR(datos)
[1] 242.5
range(datos)
[1]   10 2423
diff(range(datos))
[1] 2413

Medidas de formas

Usamos libreria moments, para el sesgo y kurtosis

library(moments)
sesgo=skewness(datos)
sesgo
[1] 3.136381

sesgo > 0: La curva es asimétricamente positiva por lo que los valores se tienden a reunir más en la parte izquierda que en la derecha de la media

kurtosis_1=kurtosis(datos)
kurtosis_1
[1] 13.58214

kurtosis_1 > 3 la distribución es Leptocúrtica

Estimación de los parámetros de la densidad.

Podemos estimar los parámetros mediante MLE, Momentos, erc. Usamos la siguiente libreria

library(fitdistrplus)

Grafica el sesgo vs curtosis

descdist(datos)
summary statistics
------
min:  10   max:  2423 
median:  104.5 
mean:  272.675 
estimated sd:  461.1389 
estimated skewness:  3.259925 
estimated kurtosis:  15.20117 

Ajusta distribuciones univariadas

f1 <- fitdist(datos,"pareto")
$start.arg
$start.arg$shape
[1] 3.118228

$start.arg$scale
[1] 577.5879


$fix.arg
NULL
plot(f1)

f1
Fitting of the distribution ' pareto ' by maximum likelihood 
Parameters:

Aparentemente el ajuste es bueno para una Pareto (1.83, 247.60)

Verificación del ajuste

Realizado a los datos Mediante pruebas de hipótesis, qqplots, curvas de densidad, etc.

Pruebas de Hipotesis

Ho:Pertenecen a la misma distribución continua Pareto (p value > alfa=.05)) Ha:No Pertenecen a la misma distribución continua Pareto (p value < alfa=.05)

Kolmogorov-Smirnov

test de kolmogorof smirnov

ks.test(datos,ppareto,f1$estimate[1],f1$estimate[2])
Warning in ks.test(datos, ppareto, f1$estimate[1], f1$estimate[2]) :
  ties should not be present for the Kolmogorov-Smirnov test

    One-sample Kolmogorov-Smirnov test

data:  datos
D = 0.10741, p-value = 0.7454
alternative hypothesis: two-sided

como p-value = 0.7454 > alfa=.05 Aceptamos Ho

Anderson Darling

ADGofTest::ad.test(datos,ppareto,f1$estimate[1],f1$estimate[2])

    Anderson-Darling GoF Test

data:  datos  and  ppareto
AD = 0.46132, p-value = 0.7855
alternative hypothesis: NA

Como p-value = 0.7855 > alfa=.05 Aceptamos Ho

Ajustamos la curva de una Pareto

hist(datos,freq=F,breaks = 20,col="blue")
u<-seq(min(datos),max(datos),by=.01)
lines(u,dpareto(u,f1$estimate[1],f1$estimate[2]),col="orange",lwd=2)

LS0tDQp0aXRsZTogIlRlb3JpYSBkZWwgUmllc2dvIFRhcmVhIDIiDQphdXRob3I6ICJDcnV6IE1hdGVvIERhdmlkIg0KZGF0ZTogIiINCm91dHB1dDoNCiAgaHRtbF9ub3RlYm9vazoNCiAgICB0b2M6IHRydWUNCiAgICB0b2NfZGVwdGg6IDUNCiAgICB0b2NfZmxvYXQ6DQogICAgICBjb2xsYXBzZWQ6IGZhbHNlDQogICAgICBzbW9vb3RoX3Njcm9sbDogdHJ1ZQ0KICBodG1sX2RvY3VtZW50Og0KICAgIHRvYzogdHJ1ZQ0KICBwZGZfZG9jdW1lbnQ6DQogICAgdG9jOiB0cnVlDQpsYW5nOiBlcy1FUw0KLS0tDQpBanVzdGEgdW5hIFBhcmV0byBhIGxvcyBzaWd1aWVudGVzIGRhdG9zDQpgYGB7cn0NCmxpYnJhcnkoZ2dwbG90MikNCmxpYnJhcnkoYWN0dWFyKQ0KZGF0b3M9YygxMCwgMTEsIDE1LCAyMiwgMjgsIDMwLCAzMiwgMzYsIDM4LCA0OCwNCiAgICAgICAgNTUsIDU2LCA2OCwgNjgsIDg1LCA4NywgOTQsIDEwMywgMTA0LCAxMDUsDQogICAgICAgIDEwOSwgMTE5LCAxMjEsIDEzNywgMTc4LCAxODEsIDIyNiwgMjg3LCAzMTAsIDEwNiwNCiAgICAgICAgMzkzLCA0MzgsIDU5MSwgMTA0NSwgMTIxMCwgMTIxMiwgMjQyMywgMzIxLCAzNTQsIDUxKQ0KYGBgDQoNCiMgUmVjb25vY2ltaWVudG8gZGUgbGEgZmFtaWxpYS4NCg0KRXN0byBzZSByZWFsaXphIG1lZGlhbnRlIGVzdGFkw61zdGljYSBkZXNjcmlwdGl2YTogY8OhbGN1bG8gZGUgbWVkaWRhcyBkZXNjcmlwdGl2YXMgeSBncsOhZmljb3MgZXN0YWTDrXN0aWNvcw0KYGBge3J9DQpoaXN0KGRhdG9zLGZyZXE9RixicmVha3MgPSAyMCxjb2w9ImJsdWUiKQ0KZDwtZGVuc2l0eShkYXRvcykNCmxpbmVzKGQpDQphPSBib3hwbG90KGRhdG9zKQ0KDQpgYGANClRlbmVtb3MgMyBwaWNvcywgZG9zIGRhdG9zIGF0aXBpY29zIDEyMTAsIDEyMTIsIDI0MjMsDQpEYXRvcyBzZXNnYWRvcyBhIGxhIGRlcmVjaGENCg0KIyMgTWVkaWRhcyBkZSB0ZW5kZW5jaWEgY2VudHJhbGVzDQpgYGB7cn0NCm1lYW4oZGF0b3MpDQptZWRpYW4oZGF0b3MpIA0KdmFyKGRhdG9zKQ0Kc3VtbWFyeShkYXRvcykNCnF1YW50aWxlKGRhdG9zKQ0KYGBgDQpUZW5lbW9zIGluZGljaW9zIGRlIGFzaW1ldHLDrWENCg0KIyMgTWVkaWRhcyBkZSBkaXNwZXJzacOzbg0KYGBge3J9DQp2YXIoZGF0b3MpDQpzZChkYXRvcykNCnJhbmdlKGRhdG9zKQ0KSVFSKGRhdG9zKQ0KcmFuZ2UoZGF0b3MpDQpkaWZmKHJhbmdlKGRhdG9zKSkNCmBgYA0KIyMgTWVkaWRhcyBkZSBmb3JtYXMNClVzYW1vcyBsaWJyZXJpYSBtb21lbnRzLCBwYXJhIGVsIHNlc2dvIHkga3VydG9zaXMNCmBgYHtyfQ0KbGlicmFyeShtb21lbnRzKQ0Kc2VzZ289c2tld25lc3MoZGF0b3MpDQpzZXNnbw0KYGBgDQpzZXNnbyA+IDA6IExhIGN1cnZhIGVzIGFzaW3DqXRyaWNhbWVudGUgcG9zaXRpdmEgcG9yIGxvIHF1ZSBsb3MgdmFsb3JlcyBzZSB0aWVuZGVuIGEgcmV1bmlyIG3DoXMgZW4gbGEgcGFydGUgaXpxdWllcmRhIHF1ZSBlbiBsYSBkZXJlY2hhIGRlIGxhIG1lZGlhDQpgYGB7cn0NCmt1cnRvc2lzXzE9a3VydG9zaXMoZGF0b3MpDQprdXJ0b3Npc18xDQpgYGANCmt1cnRvc2lzXzEgPiAzIGxhIGRpc3RyaWJ1Y2nDs24gZXMgTGVwdG9jw7pydGljYQ0KDQojIEVzdGltYWNpw7NuIGRlIGxvcyBwYXLDoW1ldHJvcyBkZSBsYSBkZW5zaWRhZC4NClBvZGVtb3MgZXN0aW1hciBsb3MgcGFyw6FtZXRyb3MgbWVkaWFudGUgTUxFLCBNb21lbnRvcywgZXJjLg0KVXNhbW9zIGxhIHNpZ3VpZW50ZSBsaWJyZXJpYQ0KYGBge3J9DQpsaWJyYXJ5KGZpdGRpc3RycGx1cykNCmBgYA0KR3JhZmljYSBlbCBzZXNnbyB2cyBjdXJ0b3Npcw0KDQpgYGB7cn0NCmRlc2NkaXN0KGRhdG9zKQ0KYGBgDQpBanVzdGEgZGlzdHJpYnVjaW9uZXMgdW5pdmFyaWFkYXMNCmBgYHtyfQ0KZjEgPC0gZml0ZGlzdChkYXRvcywicGFyZXRvIikNCnBsb3QoZjEpDQpmMQ0KYGBgDQpBcGFyZW50ZW1lbnRlIGVsIGFqdXN0ZSBlcyBidWVubyBwYXJhIHVuYSBQYXJldG8gKDEuODMsIDI0Ny42MCkNCg0KIyBWZXJpZmljYWNpw7NuIGRlbCBhanVzdGUNClJlYWxpemFkbyBhIGxvcyBkYXRvcyBNZWRpYW50ZSBwcnVlYmFzIGRlIGhpcMOzdGVzaXMsIHFxcGxvdHMsIGN1cnZhcyBkZSBkZW5zaWRhZCwgZXRjLg0KDQojIyBQcnVlYmFzIGRlIEhpcG90ZXNpcw0KSG86UGVydGVuZWNlbiBhIGxhIG1pc21hIGRpc3RyaWJ1Y2nDs24gY29udGludWEgUGFyZXRvIChwIHZhbHVlID4gYWxmYT0uMDUpKQ0KSGE6Tm8gUGVydGVuZWNlbiBhIGxhIG1pc21hIGRpc3RyaWJ1Y2nDs24gY29udGludWEgUGFyZXRvIChwIHZhbHVlIDwgYWxmYT0uMDUpDQoNCiMjIyBLb2xtb2dvcm92LVNtaXJub3YNCnRlc3QgZGUga29sbW9nb3JvZiBzbWlybm92DQpgYGB7cn0NCmtzLnRlc3QoZGF0b3MscHBhcmV0byxmMSRlc3RpbWF0ZVsxXSxmMSRlc3RpbWF0ZVsyXSkNCmBgYA0KY29tbyBwLXZhbHVlID0gMC43NDU0ID4gYWxmYT0uMDUNCkFjZXB0YW1vcyBIbw0KDQojIyMgQW5kZXJzb24gRGFybGluZw0KYGBge3J9DQpBREdvZlRlc3Q6OmFkLnRlc3QoZGF0b3MscHBhcmV0byxmMSRlc3RpbWF0ZVsxXSxmMSRlc3RpbWF0ZVsyXSkNCmBgYA0KQ29tbyBwLXZhbHVlID0gMC43ODU1ID4gYWxmYT0uMDUNCkFjZXB0YW1vcyBIbw0KDQojIyBBanVzdGFtb3MgbGEgY3VydmEgZGUgdW5hIFBhcmV0bw0KDQpgYGB7cn0NCmhpc3QoZGF0b3MsZnJlcT1GLGJyZWFrcyA9IDIwLGNvbD0iYmx1ZSIpDQp1PC1zZXEobWluKGRhdG9zKSxtYXgoZGF0b3MpLGJ5PS4wMSkNCmxpbmVzKHUsZHBhcmV0byh1LGYxJGVzdGltYXRlWzFdLGYxJGVzdGltYXRlWzJdKSxjb2w9Im9yYW5nZSIsbHdkPTIpDQpgYGANCg0K