Valores |
Frecuencias |
0 |
100 |
1 |
250 |
2 |
300 |
3 |
500 |
4 |
450 |
5 |
2000 |
#Importación de librerías para resolver los problemas.
%pylab inline
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
from __future__ import division
Populating the interactive namespace from numpy and matplotlib
xi = np.array([0,1,2,3,4,5])
fi = np.array([100,250,300,500,450,2000])
f_acum = np.cumsum(fi)
Para generar el conjunto de datos utilizamos un bucle for
creando en primer lugar un array vacío (mejora de velocidad)
data = np.empty([fi.sum(),1])
for i in range(len(xi)):
if i<1:
data[0:f_acum[i]] = xi[i]
else:
data[f_acum[i-1]:f_acum[i]] = xi[i]
data
array([[ 0.], [ 0.], [ 0.], ..., [ 5.], [ 5.], [ 5.]])
# Para trabajar con pandas creamos el dataframe
df = pd.DataFrame(data)
df
<class 'pandas.core.frame.DataFrame'> Int64Index: 3600 entries, 0 to 3599 Data columns (total 1 columns): 0 3600 non-null values dtypes: float64(1)
df.mean() #media de la columna 1
0 3.930556 dtype: float64
df.describe()
0 | |
---|---|
count | 3600.000000 |
mean | 3.930556 |
std | 1.446714 |
min | 0.000000 |
25% | 3.000000 |
50% | 5.000000 |
75% | 5.000000 |
max | 5.000000 |
Pandas nos da bastantes datos estadísticos pero si necesitamos conocer más, usamos el modulo de SciPy (stats)
m_geometrica = stats.gmean(data) print("Media geométrica = "+ str(m_geometrica))
La media armónica no se puede determinar en las distribuciones con algunos valores iguales a cero, veamos que ocurre si intentamos obtenerla:
m_armonica = stats.hmean(data) print('Media armónica = '+ str(m_armonica))
Efectivamente, obtenemos un ValureError indicando lo anterior.
Para obtener la moda, hacemos uso de las herramientas de estadística de SciPy, previamente se calcuó con la función describe de pandas
(el valor de la mediana corresponde con el que indica el 50%):
moda = stats.mode(data)
print('Moda = '+ str(moda[0]))
Moda = [[ 5.]]
Como se observa, existen diferencias singificativas entre los valores de las medias obtenidos, para tener una mejor visualización de los datos se representa un diagrama de barras.
plt.show()
plt.bar(xi, fi)
plt.title('Diagrama de barras')
plt.ylim(0,2500)
plt.xlabel(u'días de mejora '+'$x_i$')
plt.ylabel('Frecuencias '+'$f_i$')
plt.xticks(xi+0.4,xi)
([<matplotlib.axis.XTick at 0x1bd2b5f8>, <matplotlib.axis.XTick at 0x1bffab70>, <matplotlib.axis.XTick at 0x1c20fc18>, <matplotlib.axis.XTick at 0x1c211128>, <matplotlib.axis.XTick at 0x1c211780>, <matplotlib.axis.XTick at 0x1c211dd8>], <a list of 6 Text xticklabel objects>)
Como se observa en el gráfico, el valor predominante corresponde al de 5 días de mejora. No obstante, el resto de resultados desplaza el valor de la media a uno menor. Sólo con la media y sin conocer los datos, pensaríamos que la mayoría de días de mejora fue de 4 días cuando en realidad era de 5. Por tanto, el valor más representativo es el de la mediana ya que es una medida estadísitca más robusta.
porc_pacientes_exito = fi[5]/fi.sum()*100
print('Resultado '+str(round(porc_pacientes_exito,2))+'%')
Resultado 55.56%
El coeficiente de variación se calcula cuando se desea comparar las dispersiones entre dos o más distribuciones. En este caso sólo tenemos una por lo que carece de sentido hallarlo.
El decil 3 (Percentil 30%) es uno de los 9 valores (Deciles) que divide los datos en 10 partes iguales. Así:
D_3 = np.sum(fi)*0.3+1/2
print(data[D_3])
[ 3.]
SciPy también incluye una función para obtener el percentil dado:
stats.scoreatpercentile(data, 30)
3.0
Tiempo de reacción |
Nº de pacientes |
0-10 |
300 |
10-20 |
500 |
20-30 |
400 |
30-40 |
500 |
40-60 |
300 |
t_inf = np.array([0, 10, 20, 30, 40])
t_sup = np.array([10, 20, 30, 40, 60])
n_pacientes = np.array([300, 500, 400, 500, 300])
freq_pacientes = n_pacientes
Calculamos las marca de clase:
x_i = (t_sup+t_inf)/2
x_i
array([ 5., 15., 25., 35., 50.])
Y la longitud del intervalo de clase
intervalo = t_sup-t_inf
intervalo
array([10, 10, 10, 10, 20])
Para obtener el $n$
n = np.sum(freq_pacientes)
n
2000
Calculamos la frecuencia acumulada de pacientes
f_pacientes_acum = np.cumsum(n_pacientes)
f_pacientes_acum
array([ 300, 800, 1200, 1700, 2000])
Y su frecuencia relativa junto con la frecuencia relativa acumulada
freq_rel = freq_pacientes/np.sum(freq_pacientes)
freq_rel
# np.sum(freq_rel) ## Comprobación
array([ 0.15, 0.25, 0.2 , 0.25, 0.15])
freq_acum_rel = np.cumsum(freq_rel)
freq_acum_rel
array([ 0.15, 0.4 , 0.6 , 0.85, 1. ])
El cálculo de las medias se hará mediante el uso de las fórmulas.
Media aritmética:
media_aritmetica = np.sum(freq_pacientes*x_i)/n
media_aritmetica
25.75
Media geométrica
media_geometrica = (np.prod(x_i**freq_pacientes))**(1/n)
media_geometrica
inf
Media cuadrática
media_cuadratica = np.sqrt((np.sum(freq_pacientes*x_i))/n)
media_cuadratica
5.0744457825461096
Media armónica
media_armonica = np.sum(freq_pacientes)/np.sum(freq_pacientes/x_i)
media_armonica
15.429831006612782
Estos mismos resultados se pueden obtener directamente si creamos un data_frame
data2 = np.zeros([freq_pacientes.sum(),1])
for j in range(len(x_i)):
if j<1:
data2[0:f_pacientes_acum[j]] = x_i[j]
else:
data2[f_pacientes_acum[j-1]:f_pacientes_acum[j]] = x_i[j]
data2
array([[ 5.], [ 5.], [ 5.], ..., [ 50.], [ 50.], [ 50.]])
df2 = pd.DataFrame(data2)
df2.describe()
0 | |
---|---|
count | 2000.00000 |
mean | 25.75000 |
std | 14.25795 |
min | 5.00000 |
25% | 15.00000 |
50% | 25.00000 |
75% | 35.00000 |
max | 50.00000 |
stats.hmean(data2) #media armonica de SciPy
array([ 15.42983101])
stats.gmean(data2) #media geometrica de SciPy
array([ 20.86021734])
plt.hist(data2, bins=(0,10,20,30,40,60))
plt.title('Histograma')
plt.ylim(0,600)
plt.xlabel(u'Tiempo de reacción '+'$x_i$')
plt.ylabel(u'nº de pacientes '+'$f_i$')
<matplotlib.text.Text at 0x222aa7f0>
Dada la distribución de los datos, el uso de la media aritmética en este caso no supone mucho error y de hecho coincide con la mediana. La media armónica se desplaza a la izquierda debido a su sensibilidad con valores pequeños al conjunto. La media geométrica en este caso puede ser también útil. Pese a todo ello, hay que descatar que las medidas de tendencia central previamente nombradas, proporcionan sólo información parcial y se requieren medidas de dispersión para determinar su calidad. Por ejemplo, la desviación estándar es elevada $\sigma = 14.26$ respecto al valor de la media, lo que indica que la dispersión de la muestra es elevada. La construcción geométrica para hallar las medias aritmética (A), geométrica (G), armónica (H) y cuadrática (Q) de dos números $a$ y $b$. Este diagrama ayuda a comprender la relación entre las mismas.
from IPython.display import Image
Image(url='http://upload.wikimedia.org/wikipedia/commons/thumb/f/f7/MathematicalMeans.svg/500px-MathematicalMeans.svg.png')
Dado que los valores en tiempo de reacción se dan en intervalos y que hasta ahora se han considerado marcas en el valor medio, se coniderará que la mitad del Nº de pacientes del intervalo 30-40' está por encima de 35'. Así, el número:
n_pacientes_medicamento2 = n_pacientes[3]/2+n_pacientes[4]
n_pacientes_medicamento2
550.0
Como ya se comentó, la distribución de tiempo (ver histograma) es plana por lo que una medida central no será representativa. El valor elevado de la desviación estándar lo indica a su vez.
Para comprender mejor la representatividad de los valores centrales, se puede observar el siguiente gráfico que representa dos distribuciones. La centrada tiene valores próxmias de media, mediana y moda así como un valor bajo de desv. estándar. Sin embargo, la distribución más alargada (línea discontínua) representa claras diferencias entre sus valores centrales y como consecuencia una desv. estándar elevada.
Image(url='http://upload.wikimedia.org/wikipedia/commons/thumb/d/de/Comparison_mean_median_mode.svg/500px-Comparison_mean_median_mode.svg.png')
Para poder comparar dispersiones entre las distribuciones de los tres medicamentos se debe de calcular el coeficiente de variación de Pearson.
# Medicamento 1: Sindolorcabezon
media_med1 = np.mean(data2)
std_med1 = np.std(data2)
CV_med1 = std_med1/media_med1*100
print('SINDOLORCABEZON: media: '+str(media_med1)+', CV: '+str(round(CV_med1,2)))
# Medicamento 2: Palacabeza
media_med2 = 25
std_med2 = np.sqrt(200)
CV_med2 = std_med2/media_med2*100
print('PALACABEZA: media: '+str(media_med2)+', CV: '+str(round(CV_med2,2)))
# Medicamento 3: Sinjaquecahoy
media_med3 = 30
std_med3 = np.sqrt(300)
CV_med3 = std_med3/media_med3*100
print('SINJAQUECAHOY: media: '+str(media_med3)+', CV: '+str(round(CV_med3,2)))
SINDOLORCABEZON: media: 25.75, CV: 55.36 PALACABEZA: media: 25, CV: 56.57 SINJAQUECAHOY: media: 30, CV: 57.74
De los resultados se observa que SINDOLORCABEZON y PALACABEZA obtienen resultados similares. Pese a que a primera vista PALACABEZA tenga un mejor tiempo las dispersión es elevada. Por lo que estadisticamente y con los datos facilitados sería díficil hacer una selección apropiada. En cualquier caso y dado los tiempos que corren, la recomendación sería el que suponga menor coste de los dos a la Seguridad Social.