A finales del siglo XIX el físico norteamericano Newbold descubrió que la proporción de datos que empiezan por una cifra \(d\), \(p(d)\), en listas de datos correspondientes a muchos fenómenos naturales y demográficos es aproximadamente: \[ p(d) = \log_{10}\left(\frac{d+1}{d}\right), \ \ d=1,2,\ldots,9. \] Por ejemplo, \(p(1)=\log_{10} 2 \approx 0.301030\) es la frecuencia relativa de datos que empiezan por 1. A raíz de un artículo publicado en 1938 por Benford, la fórmula anterior se conoce como ley de Benford. Con el siguente código R dibujamos la función de masa de esta distribución:

d = 1:9 # Numeros del 1 al 9
PBenford = log10((d+1)/d) # Fn de masa Ley de Benford
plot(d,PBenford,pch=19,cex=2,
     ylab="Función de masa Ley de Benford",
     cex.axis=1.5,cex.lab=1.5, xaxt="n")
axis(1,at=d,labels=d,cex.axis=1.5)

El fichero poblacion.RData incluye un fichero llamado poblaciones con la población total de los municipios españoles, así como su población de hombres y de mujeres.

a) Contrasta a nivel \(\alpha=0.05\) la hipótesis nula de que la población total se ajusta a la ley de Benford.

b) Repite el ejercicio pero considerando sólo los municipios de más de 1000 habitantes.

c) Considera las poblaciones totales (de los municipios con 10 o más habitantes) y contrasta a nivel \(\alpha=0.05\) la hipótesis nula de que el primer dígito es independiente del segundo.

(Indicación: Puedes utilizar, si te sirven de ayuda, las funciones del fichero benford.R).

Solución:

a) Sea \(D\) el primer dígito de la población total en un municipio español. Planteamos el contraste de bondad de ajuste:

\(H_{0}:\) \(D\) satisface la ley de Benford

source("benford.R")

load("poblacion.RData")
# Población total en los municipios españoles:
PobT = poblaciones$pobtotal
# Frecuencias observadas de primeros dígitos:
OT = benford(PobT)
OT

## [1] 2554 1441 1023  741  637  531  462  379  344

Podemos hacer los cálculos “a mano”:

# Tamaño muestral:
nT = sum(OT)
nT

## [1] 8112

# Frecuencias esperadas bajo la ley de Benford ($H_{0}$):
eT = nT*PBenford
eT

## [1] 2441.9553 1428.4523 1013.5030  786.1340  642.3183  543.0724  470.4307
## [8]  414.9493  371.1848

# Estadístico del contraste ji-cuadrado:
TestEst = sum((OT-eT)^2/eT)
TestEst

## [1] 13.50036

Como \(\chi^2_{8;0.05}=15.51\), no rechazamos \(H_{0}\) a nivel \(\alpha=0.05\).

También podemos hacer el contraste de bondad de ajuste con R:

chisq.test(OT,p=PBenford)

## 
##  Chi-squared test for given probabilities
## 
## data:  OT
## X-squared = 13.5, df = 8, p-value = 0.09575

barplot(OT/nT,names.arg=d,main="Poblaciones totales",space=0)
lines(d-1/2,PBenford,type="b")

b)

PobT1000 = PobT[PobT>1000]
OT1000 = benford(PobT1000)
chisq.test(OT1000,p=PBenford)

## 
##  Chi-squared test for given probabilities
## 
## data:  OT1000
## X-squared = 298.91, df = 8, p-value < 2.2e-16

c)

# Poblaciones de más de 10 habitantes:
Masde10 = PobT>=10
PobT10 = PobT[Masde10]
# Frecuencias observadas de primer y segundo dígito
OT10 = benford3(PobT10) # Tabla de contingencia
OT10

##     0   1   2   3   4   5   6   7   8   9
## 1 377 355 283 254 266 207 216 206 201 189
## 2 172 211 159 157 136 144 107 144 116  95
## 3 126 138 113  75 113  95  96  85  96  86
## 4  88  82  81  84  66  61  67  75  73  64
## 5  78  69  71  70  67  51  61  52  61  57
## 6  65  53  61  53  73  45  54  43  43  40
## 7  59  48  45  49  44  49  56  42  38  32
## 8  42  48  39  40  30  38  40  39  31  31
## 9  29  39  35  26  28  24  30  48  45  39

chisq.test(OT10)

## 
##  Pearson's Chi-squared test
## 
## data:  OT10
## X-squared = 120.52, df = 72, p-value = 0.0002974