[R-es] Clasificacion de individuos
jluis.gilsanz en tasacionesh.com
jluis.gilsanz en tasacionesh.com
Mie Oct 1 14:23:39 CEST 2014
Hola de nuevo:
Tirando por el camino de la distribución Beta,(el de encontrar la
distribución subyacente porque de estimación Bayesiana como me sugería
Carlos estoy absolutamente pez), he alcanzado el siguiente punto.
He intentado ajustar una Beta de diversas formas:
1.-Usando el método de los momentos mediante la función rriskFitdist.cont
del paquete rriskDistributions (gracias por descubrirmelo Carlos)
Resultado: una B(0.1096258 , 1.5452894). Generada una simulación con esa
distribución KS rechaza el ajuste.
2.-Usando el método de los momentos mediante una estimación "manual" de
los parámetros.
Resultado: una B(0.1088774 , 1.53474). Generada una simulación con esa
distribución KS rechaza el ajuste.
3.-Usando el método de la máxima verosimilitud mediante la función
rriskFitdist.cont del paquete rriskDistributions.
No hay resultado puesto que no es posible obtener los parámetros por este
método. Informa de un error 100 que no tengo ni idea de a que se refiere.
Por otra parte he intentado ajustar, sin la premisa de que sea una Beta,
alguna otra distribución a mis datos así:
4.-Usando el método de los percentiles (deciles en mi caso) mediante la
función fit.perc del paquete rriskDistributions
Resultado: Aunque en la parte gráfica, de la GUI que la función genera,,
sugiere una Beta, en realidad no se "queda" con ninguna distribución.
5.-Usando la funcion fit.cont del paquete rriskDistributions.
En la GUI de la función sugiere una Normal aun cuando por KS se rechaza el
ajuste en esa y el resto de distribuciones ajustadas.
La cuestión es que he probado con una Cauchy porque tanto por el valor
del estadístico de KS (el menor de todos) como por los gráficos mostrados
parecía la menos mala.
Al generar una simulación con los parámetros obtenidos C(0.02681035
0.01475054) y ejecutar el test de KS tengo un resultado de que mi
variable si ajusta a esa distribución, contradiciendo el resultado previo
mostrado por la función y dejándome un poco "descolocado".
Visto lo visto parece ser que mi variable V3 finalmente no es una Beta si
no una Cauchy lo cual me deja aun mas confuso.
Aun me queda atacar el proyecto vía análisis cluster que tiene pinta de
que va a ser la opción implementada finalmente pero el tema del ajuste de
una distribucion subyacente con R me ha dejado bastante confuso.
La sintaxis de R, sus resultados y el fichero txt con los datos van a
continuación:
Mil gracias por llegar hasta aquí ;-)
> #####Cargamos datos y librerias#
> library(rriskDistributions)
> library(lattice)
> Datos <-read.delim("c:\\Datos_2.txt",header=TRUE,
dec=",",stringsAsFactors =FALSE)
> v3<-Datos$v3
>
>
> #####Histograma vbl original y estadisticos basicos#
> histogram(v3, type="percent", breaks=100, endpoints=c(0,1))
> mu <- mean(v3)
> sigma <-var(v3)
> num <- length(v3)
> pctiles <- quantile(v3, c(0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1))
> #####Funcion de optencion de alpha y beta según momentos "manual"#
> estBetaParams <- function(mu, var)
+ {
+ alpha <- ((1 - mu) / var - 1 / mu) * mu ^ 2
+ beta <- alpha * (1 / mu - 1)
+ return(params = list(alpha = alpha, beta = beta))
+ }
> #####Obtencion de ditribucion beta por metodo de momentos usando paquete
rriskDistributions#
> BetaAj_fit_mme <- rriskFitdist.cont(v3,"beta", method="mme")
> #####Generamos Beta simulada con parametros obtenidos#
> BetaAj_fit_mme_data <-
rbeta(num,BetaAj_fit_mme$estimate[1],BetaAj_fit_mme$estimate[2])
> histogram(BetaAj_fit_mme_data, type="percent", breaks=100,
endpoints=c(0,1))
> #####Comprobamos parametros por calculo "manual" de parametros por
momentos#
> BetaAj_manual_params <- estBetaParams(mu,sigma)
> BetaAj_manual_params
$alpha
[1] 0.1088774
$beta
[1] 1.53474
> BetaAj_fit_mme$estimate
shape1 shape2
0.1096258 1.5452894
> #####Los parametros estimados por rriskDistribution y de forma manual
son aprox identicos#
> #####alpha-> 0.1088774 vs0.1096258 beta -> 1.53474 vs 1.5452894 #
> #####Generamos Beta simulada con parametros obtenidos por momentos de
forma "manual"#
> BetaAj_manual_data <-
rbeta(num,BetaAj_manual_params$alpha,BetaAj_manual_params$beta )
> histogram(BetaAj_manual_data, type="percent", breaks=100,
endpoints=c(0,1))
> #####Comprobamos bondad de ajuste de Betas simuladas por metodo de
momentos generadas por rriskDistribution y de forma manual con nuestra v3#
> ks.Aj_fit_mme<- ks.test(v3,BetaAj_fit_mme_data)
Warning message:
In ks.test(v3, BetaAj_fit_mme_data) :
p-value will be approximate in the presence of ties
> ks.Aj_manual <- ks.test(v3,BetaAj_manual_data)
Warning message:
In ks.test(v3, BetaAj_manual_data) :
p-value will be approximate in the presence of ties
> ks.Aj_fit_mme
Two-sample Kolmogorov-Smirnov test
data: v3 and BetaAj_fit_mme_data
D = 0.5234, p-value < 2.2e-16
alternative hypothesis: two-sided
> ks.Aj_manual
Two-sample Kolmogorov-Smirnov test
data: v3 and BetaAj_manual_data
D = 0.4851, p-value < 2.2e-16
alternative hypothesis: two-sided
> #####KS rechazado en ambos casos. Nuestra v3 no ajusta con las dos
simulaciones generadas#
> #####Obtencion de ditribucion beta por maxima verosimilitud#
> BetaAj_fit_mle <- rriskFitdist.cont(v3,"beta", method="mle")
Error in rriskFitdist.cont(v3, "beta", method = "mle") :
the function mle failed to estimate the parameters,
with the error code 100
> BetaAj_fit_mle
Error: object 'BetaAj_fit_mle' not found
> #####No se puede obtener Beta por metodo maxima verosimilitud#
> #####Obtencion de ditribucion por percentiles usando rriskDistribution
> Beta_pctil <- fit.perc(
c(0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1),pctiles)
-----------------------------------------------------------------------
Begin fitting distributions...
Truncated normal distribution has been fitted successfully!
Warning: non-central chi-square distribution could not be fitted!
Warning: Pert distribution could not be fitted!
Warning: Triangular distribution could not be fitted!
Error: INVALID INPUT, percentiles are out of the domain (0,inf) =>
Gompertz distribution couldn't be fitted!
Warning: Gompertz distribution could not be fitted!
Normal distribution has been fitted successfully!
Beta distribution has been fitted successfully!
Cauchy distribution has been fitted successfully!
Warning: Chi-square distribution could not be fitted!
Logistic distribution has been fitted successfully!
Warning: Student's t distribution could not be fitted!
Exponential distribution has been fitted successfully!
Error: INVALID INPUT, percentiles are out of the domain [0,inf) => F
distribution couldn't be fitted!
Warning: F distribution could not be fitted!
Gamma distribution has been fitted successfully!
Weibull distribution has been fitted successfully!
Error: INVALID INPUT, percentiles are out of the domain (0,inf) =>
Lognormal distribution couldn't be fitted!
Warning: Lognormal distribution could not be fitted!
Error: INVALID INPUT, uniform distribution can be fitted only by TWO given
percentiles!
Warning: Uniform distribution could not be fitted!
End fitting distributions...
-----------------------------------------------------------------------
Chosen continuous distribution is: NA
Fitted parameters are:
[1] NA
--------------------------------------------------------------------
> #####Segun la parte grafica de la GUI elige una Beta como distribucion
ajustada#
> Beta_pctil
$`p/q`
p q
0% 0.0 0.0000
10% 0.1 0.0000
20% 0.2 0.0160
30% 0.3 0.0190
40% 0.4 0.0240
50% 0.5 0.0300
60% 0.6 0.0360
70% 0.7 0.0430
80% 0.8 0.0638
90% 0.9 0.1000
100% 1.0 1.0000
$chosenDistr
[1] NA
$fittedParams
[1] NA
> #####Realmente no elige ninguna distribucion
> #####Obtencion de ditribucion continua que mejor ajusta usando
rriskDistribution
> Ajust_Cont <- fit.cont(v3)
-------------------------------------------------------------------
Begin fitting distributions...
Normal distribution has been fitted successfully!
Cauchy distribution has been fitted successfully!
Logistic distribution has been fitted successfully!
Warning: Beta distribution couldn't be fitted!
Exponential distribution has been fitted successfully!
Warning: Chi-square distribution couldn't be fitted!
Uniform distribution has been fitted successfully!
Warning: gamma distribution couldn't be fitted!
Warning: Lognormal distribution couldn't be fitted!
Warning: Weibull distribution couldn't be fitted!
Warning: F distribution couldn't be fitted!
Student distribution has been fitted successfully!
Gompertz distribution has been fitted successfully!
Warning: Triangular distribution couldn't be fitted!
End fitting distributions...
-------------------------------------------------------------------
logL AIC BIC Chisq(value) Chisq(p) AD(value) H(AD)
KS(value) H(KS)
Normal 108.28 -212.55 -205.63 851.81 0 50.62
rejected 0.33 rejected
Cauchy 411.8 -819.61 -812.69 41.42 0 9.18
rejected 0.16 rejected
Logistic 240.66 -477.31 -470.39 243.54 0 24.42
rejected 0.27 rejected
Exponential 402.89 -803.78 -800.32 Inf 0 Inf
rejected 0.23 rejected
Uniform NULL NULL NULL 1206.99 0 Inf NULL
0.37 rejected
Student -219.33 440.65 444.11 6548.81 0 78.14 NULL
0.50 rejected
Gompertz 402.88 -801.76 -794.84 Inf 0 Inf NULL
0.23 rejected
-------------------------------------------------------------------
> #####KS rechazado en todos los casos. Nuestra variable original no
ajusta con ninguna de las distribuciones probadas #
> #####La funcion sugiere una Normal como la menos mala pero graficamente
y por el valor de KS la Cauchy parece la mas adecuada#
-------------------------------------------------------------------
Chosen continuous distribution is: Cauchy (cauchy)
Fitted parameters are:
location scale
0.02681035 0.01475054
--------------------------------------------------------------------
> #####Generamos una Cauchy simulada y probamos su ajuste #
> Ajust_cau <- rcauchy(235,0.0268103,0.01475054)
> histogram(Ajust_cau, type="percent", breaks=100, endpoints=c(0,1))
> ks.test(v3,Ajust_cau)
Two-sample Kolmogorov-Smirnov test
data: v3 and Ajust_cau
D = 0.1149, p-value = 0.0899
alternative hypothesis: two-sided
> #####Aparentemente KS acepta el ajuste. ¿Porque en la GUI de fit.cont
aparece como rechazado? #
{In Archive} Re: [R-es] Clasificacion de individuos
(Internet)
jorgeivanvelez
To:
ihidalgo
Cc:
Jluis GILSANZ, r-help-es
30/09/2014 15:41
Archive:
This message is being viewed in an archive.
Me parece muy buena la idea de Isidro. En
http://cran.r-project.org/web/views/Cluster.html hay infinidad de metodos
para realizar la clasificacion que quieres. Otra posibilidad es utilizar
un modelo de clases latentes o un cluster jerarquico usando el paquete
pvclust. En ninguno de los casos la normalidad es "necesaria".
Por otro lado, si tus variables son porcentajes, por que no trabajas con
la distribucion Beta? En R estan implementados por supuesto la densidad,
los percentiles y los numeros aleatorios (ver ?rbeta). Si quieres hacer
regresion, el paquete betareg es excelente y vale la pena darle una mirada
al articulo en JSS.
Si estamos todos perdidos, quizas un poco de contexto nos ayude a
ayudarte.
Saludos cordiales,
Jorge.-
2014-09-30 20:40 GMT+10:00 Isidro Hidalgo <ihidalgo en jccm.es>:
¿Has probado a hacer directamente una clasificación de los individuos con
3
clusters? Es muy sencillo y quizá te funcione. En un espacio
tridimensional un
cluster se hace más rápido que un disparo.
Un saludo.
Isidro
> -----Mensaje original-----
> De: r-help-es-bounces en r-project.org [mailto:r-help-es-bounces en r-
> project.org] En nombre de jluis.gilsanz en tasacionesh.com
> Enviado el: martes, 30 de septiembre de 2014 12:24
> Para: r-help-es en r-project.org
> Asunto: [R-es] Clasificacion de individuos
>
> Estimados apa eRos:
>
> La duda o propuesta que os voy a plantear es a la vez metodol gica y
> relacionada con R.
>
> Me encuentro trabajando con tres variables que son el resultado de un
> computo de porcentajes.
> Me explico, se toma una muestra de n casos (unos 6.500 aprox)
> pertenecientes a i individuos (unos 230 aprox) en la que se comprueba
> si un determinado evento ha ocurrido o no, anot ndose 1 en caso de
> dicha ocurrencia y 0 en caso de no ocurrencia.
> Algo as como:
> indiv ocurrencia
> -------- -----------------
> 1 0
> 1 0
> 1 1
> 2 0
> 2 1
> 3 0
> 3 0
> 3 0
> 4 1
> 4 1
> . .
> . .
> . .
> n
>
> Tras ello se computa el porcentaje de ocurrencias para cada individuo
> obteniendo las variables que ser n estudiadas, obteniendo algo as :
>
> indiv %
> ----- ------
> 1 0.333
> 2 0.5
> 3 0
> 4 1
> . .
> . .
> i
>
> Este mismo proceso se repite en cada una de las tres variables objeto
> de estudio.
>
> Lo que se pretende es clasificar los i individuos en tres grupos seg n
> sus resultados en cuanto a los porcentajes calculados (%):
> -Por debajo de la media: Individuos que forman parte de la cola con
> peores porcentajes. Aprox 10% -Por encima de la media: Individuos que
> forman parte de la cola con peores porcentajes. Aprox 10% -Acordes a
> la media: El 80% de individuos resultantes.
>
> Se trata b sicamente de "rega ar" a los del primer grupo y "felicitar"
> los del segundo grupo ;-)
>
> La cuesti n es que de las tres variables en estudio, las dos ultimas no
> son normales:
>
> >stem(v1)
>
> The decimal point is 1 digit(s) to the left of the |
>
> 0 | 0000000000466899
> 1 | 0133347777778999
> 2 | 0000011233344555667778889999
> 3 | 0001233333333334444567778888889999999999
> 4 | 000001122233333344444566788889999
> 5 | 000000000000011234444566667777889
> 6 | 00122233345555777777788899
> 7 | 00011222334455567779
> 8 | 1333336668
> 9 |
> 10 | 0000000000000
>
> >quantile(v1,c(0.1,0.9))
> 10% 90%
> 0.1670 0.7834
>
>
>
> >stem(v2)
>
> The decimal point is 1 digit(s) to the left of the |
>
> 0 |
> 00000000000000000000000000000000000000000000000000000000000000000000+75
> 1 | 00000000111122233333444445677788888999
> 2 | 0000122444557899
> 3 | 0001123378999
> 4 | 266
> 5 | 0000
> 6 | 57
> 7 |
> 8 | 3
> 9 |
> 10 | 000
>
> >quantile(v2,c(0.1,0.9))
> 10% 90%
> 0.000 0.304
>
>
>
> stem(V3)
>
> The decimal point is 1 digit(s) to the left of the |
>
> 0 |
> 00000000000000000000000000000000011111111122222222222222222222222222+12
> 8
> 1 | 000001133333477
> 2 | 000
> 3 | 3
> 4 |
> 5 | 000
> 6 |
> 7 |
> 8 |
> 9 |
> 10 | 00000
>
> >quantile(V3,c(0.1,0.9))
> 10% 90%
> 0.0 0.1
>
>
>
> La primera variable V1 aparece como Normal segun los test de Kolmogorv-
> Smirnov, Jarque-Bera (simetria), Agostino (simetria) y Anscombe
> (curtosis) pero como No Normal segun el test de Shapiro-Wilks.
> Las otras dos no aparecen como Normales en ninguno de los test,
> logicamente al tener una asimetria tan fuerte.
> He probado transformando mediante Box-Cox pero la ni la raiz
> cuadadrada, ni 1/Variable me solucionan el problema y al haber muchos
> casos con 0 la logaritmica tampoco me vale.
>
>
>
> Asumiendo Normalidad puedo emplear los intervalos de confianza para
> "rega ar" o "felicitar" a los individuos y Puedo utilizar los
> percentiles
> 0.1 y 0,9 para obtener las clasificaciones en cualquier caso (normal o
> no normal).
>
> Pero me surgen varias dudas:
> -La principal de todas es la metodologia (y paquetes/funciones de R
> asociados) para llevar a cabo tarea de clasificacion.Actualmente estoy
> valorando hacer la clasificacion mediante intervalos de confianza (en
> caso de poder transformar a "normales" las dos ultimas variables),
> percentiles o directamente aplicando la desigualdad Chebychev.
>
> - Que paquetes me ofrece R para obtener la distribucion de probabilidad
> subyacente en una determinada muestra?.
>
>
> - Que test, de los que se ofrecen en R, considerais como los mas
> potentes
> para comprobar si una muestra de estas caracteristicassigue una
> determinada distribucion?.Yo utilizo ks.test y shapiro.test pero igual
> hay
> otros mejores en algun paquete que conozcais.Por la construccion de las
> variables habia supuesto que las variables debian de seguir una
> distribucion U(0,1) pero el test ks solo concluye uniformidad en la
> primera variable.
>
>
> Muchas gracias a tod en s por leer hasta el final del ladrillo
>
>
> Saludos
>
>
>
>
>
>
>
>
>
>
> -- AVISO LEGAL --
>
> Los datos personales que en esta comunicaci n aparecen, as como los
> que nuestra
> empresa mantiene de Vd. y de su empresa, son tratados con la finalidad
> de mantener
> el contacto as como realizar las gestiones que en esta aparecen (Ley
> Org nica
> 15/1999, de 13 de diciembre, de Protecci n de Datos de Car cter
> Personal).
> Puede ejercer sus derechos de acceso, rectificaci n, cancelaci n y
> oposici n
> dirigi ndose a atencion.clientes en tasacionesh.com
> La utilizaci n de su direcci n de correo electr nico por parte de
> nuestra empresa
> queda sujeta a las disposiciones de la Ley 34/2002, de Servicios de la
> Sociedad de
> la Informaci n y el Comercio Electr nico. Si Vd. recibe comunicaci n
> comercial por
> nuestra parte y desea dejar de recibirla, rogamos nos lo comunique por
> v a electr nica
> a trav s de la direcci n atencion.clientes en tasacionesh.com
>
> [[alternative HTML version deleted]]
_______________________________________________
R-help-es mailing list
R-help-es en r-project.org
https://stat.ethz.ch/mailman/listinfo/r-help-es
------------ próxima parte ------------
Se ha borrado un adjunto en formato HTML...
URL: <https://stat.ethz.ch/pipermail/r-help-es/attachments/20141001/5edfb0a4/attachment-0001.html>
------------ próxima parte ------------
A non-text attachment was scrubbed...
Name: no disponible
Type: image/gif
Size: 4073 bytes
Desc: no disponible
URL: <https://stat.ethz.ch/pipermail/r-help-es/attachments/20141001/5edfb0a4/attachment-0007.gif>
------------ próxima parte ------------
A non-text attachment was scrubbed...
Name: no disponible
Type: image/gif
Size: 3655 bytes
Desc: no disponible
URL: <https://stat.ethz.ch/pipermail/r-help-es/attachments/20141001/5edfb0a4/attachment-0008.gif>
------------ próxima parte ------------
A non-text attachment was scrubbed...
Name: no disponible
Type: image/gif
Size: 3639 bytes
Desc: no disponible
URL: <https://stat.ethz.ch/pipermail/r-help-es/attachments/20141001/5edfb0a4/attachment-0009.gif>
------------ próxima parte ------------
A non-text attachment was scrubbed...
Name: no disponible
Type: image/gif
Size: 40525 bytes
Desc: no disponible
URL: <https://stat.ethz.ch/pipermail/r-help-es/attachments/20141001/5edfb0a4/attachment-0010.gif>
------------ próxima parte ------------
A non-text attachment was scrubbed...
Name: no disponible
Type: image/gif
Size: 32878 bytes
Desc: no disponible
URL: <https://stat.ethz.ch/pipermail/r-help-es/attachments/20141001/5edfb0a4/attachment-0011.gif>
------------ próxima parte ------------
A non-text attachment was scrubbed...
Name: no disponible
Type: image/gif
Size: 32675 bytes
Desc: no disponible
URL: <https://stat.ethz.ch/pipermail/r-help-es/attachments/20141001/5edfb0a4/attachment-0012.gif>
------------ próxima parte ------------
A non-text attachment was scrubbed...
Name: no disponible
Type: image/gif
Size: 3349 bytes
Desc: no disponible
URL: <https://stat.ethz.ch/pipermail/r-help-es/attachments/20141001/5edfb0a4/attachment-0013.gif>
------------ próxima parte ------------
An embedded and charset-unspecified text was scrubbed...
Name: Datos_2.txt
URL: <https://stat.ethz.ch/pipermail/r-help-es/attachments/20141001/5edfb0a4/attachment-0001.txt>
------------ próxima parte ------------
-- AVISO LEGAL --
Los datos personales que en esta comunicación aparecen, así como los que nuestra
empresa mantiene de Vd. y de su empresa, son tratados con la finalidad de mantener
el contacto así como realizar las gestiones que en esta aparecen (Ley Orgánica
15/1999, de 13 de diciembre, de Protección de Datos de Carácter Personal).
Puede ejercer sus derechos de acceso, rectificación, cancelación y oposición
dirigiéndose a atencion.clientes en tasacionesh.com
La utilización de su dirección de correo electrónico por parte de nuestra empresa
queda sujeta a las disposiciones de la Ley 34/2002, de Servicios de la Sociedad de
la Información y el Comercio Electrónico. Si Vd. recibe comunicación comercial por
nuestra parte y desea dejar de recibirla, rogamos nos lo comunique por vía electrónica
a través de la dirección atencion.clientes en tasacionesh.com
Más información sobre la lista de distribución R-help-es