Testes de modalidade e estimativas de densidade do kernel

Ideas Lab Blocked Desbloquear Seguir Seguindo 30 de dezembro de 2018

Ao processar um grande número de conjuntos de dados que podem ter diferentes distribuições de dados, somos confrontados com as seguintes considerações:

  • A distribuição de dados é unimodal e, se for o caso, qual modelo melhor se aproxima dela (distribuição uniforme, distribuição T, distribuição qui-quadrado, distribuição cauchy, etc)?
  • Se a distribuição de dados for multimodal, podemos identificar automaticamente o número de modos e fornecer estatísticas descritivas mais detalhadas?
  • Como podemos estimar a função de densidade de probabilidade de um novo conjunto de dados?

Este caderno aborda os seguintes assuntos:

  • Histogramas vs aproximação da função densidade de probabilidade
  • Estimativas de densidade de kernel
  • Escolha da largura de banda ideal: Silverman / Scott / Grid Search Cross Validation
  • Testes estatísticos para distribuições unimodais
  • Teste DIP para unimodality
  • Identificação do número de modos de distribuição de dados com base na estimativa da densidade do núcleo

Histogramas e PDFs

Como explicado neste post do blog https://mglerner.github.io/posts/histograms-and-kernel-density-estimation-kde-2.html histogramas têm a desvantagem de esconder em caixas de tamanho inadequado alguns dos detalhes do distribuição de dados real.

 def plotHistogramAndPdf (dados, x, pdf): 
ax = plt.gca ()
plt.hist (dados, caixas = 4, alfa = 0,4, label = 'histograma de valores de entrada');
plt.ylabel ('Frequency')
plt.xlabel ('x values')
ax2 = ax.twinx ()
plt.plot (x, pdf, c = 'vermelho', label = 'função de densidade de probabilidade');
plt.ylabel ('PDF')
[tl.set_color ('r') para tl em ax2.get_yticklabels ()]
ax.legend (bbox_to_anchor = (0,4, 1,15))
ax2.legend (bbox_to_anchor = (1.15,1.15))
plt.savefig ('figures / hist.jpg', bbox_inches = 'tight')

plotHistogramAndPdf (dados, x, true_pdf)

Estimativas de densidade de kernel

As estimativas de densidade do kernel dependem de uma largura de banda arbitrária que rege a suavidade da aproximação retornada. O exemplo abaixo ilustra o efeito de vários valores de largura de banda:

 def getKernelDensityEstimation (valores, x, largura de banda = 0.2, kernel = 'gaussian'): 
model = KernelDensity (kernel = kernel, largura de banda = largura de banda)
model.fit (valores [:, np.newaxis])
log_density = model.score_samples (x [:, np.newaxis])
retornar np.exp (log_density)
 para largura de banda em np.linspace (0.2, 3, 3): 
kde = getKernelDensityEstimation (dados, x, largura de banda = largura de banda)
plt.plot (x, kde, alfa = 0.8, label = f'bandwidth = {arredondado (largura de banda, 2)} ')
plt.plot (x, true_pdf, label = 'True PDF')
plt.legend ()
plt.title ('Efeito de vários valores de largura de banda n Quanto maior a largura de banda, mais suave a aproximação se torna');
plt.savefig ('figures / bw.jpg', bbox_inches = 'tight')

Métodos para selecionar a largura de ban ótima para a estimativa da densidade do kernel

Para identificar a largura de banda ideal, existem alguns tipos de abordagens:

  • Regra de ouro de Silverman: assume distribuição gaussiana para a densidade desconhecida. Não é o seletor de largura de banda mais ideal, mas é usado como um estimador muito rápido razoavelmente bom ou como um primeiro estimador em seletores de largura de banda de múltiplos estágios. Regras de plug-in solve-the-equation mais precisas usam estimativa de derivada de densidade quadrada integrada funcional para estimar a largura de banda ideal. Eles exigem altos cálculos para resolver uma equação não linear usando métodos iterativos. Eles usam o ROT como primeira estimativa
  • Regra geral de Scott: é ideal para amostras aleatórias de dados normalmente distribuídos, no sentido de que minimiza o erro quadrático médio integrado da estimativa de densidade.

Esses dois métodos têm a vantagem de serem computacionalmente rápidos, mas geralmente fornecem muito poucos escaninhos e é provável que ele ofereça a distribuição de dados subjacente. Ambos os métodos foram implementados no pacote statsmodels e ilustrados abaixo.

  • Métodos baseados em validação cruzada: statsmodels vem com um parâmetro cv bandwidth. Como alternativa, podemos implementar uma validação cruzada de pesquisa de grade. Ao contrário dos dois primeiros métodos, a realização de uma pesquisa de grade provavelmente será mais computacionalmente cara, especialmente para conjuntos de dados maiores
 de statsmodels.nonparametric.bandwidths import bw_silverman, bw_scott, select_bandwidth 
 silverman_bandwidth = bw_silverman (dados) 
 largura de banda # select permite definir um kernel diferente 
silverman_bandwidth_gauss = select_bandwidth (dados, bw = 'silverman', kernel = 'gauss')
 scott_bandwidth = bw_scott (dados) 
 def bestBandwidth (dados, minBandwidth = 0.1, maxBandwidth = 2, nb_bandwidths = 30, cv = 30): 
"" "
Executar uma pesquisa de grade de validação cruzada para identificar a largura de banda ideal para a densidade do kernel
estimativa.
"" "
de sklearn.model_selection import GridSearchCV
model = GridSearchCV (KernelDensity (),
{'largura de banda': np.linspace (minBandwidth, maxBandwidth, nb_bandwidths)}, cv = cv)
model.fit (data [:, None])
return model.best_params _ ['largura de banda']
 cv_bandwidth = bestBandwidth (dados) 
 print (f "Largura de banda de Silverman = {silverman_bandwidth}") 
print (f "largura de banda Scott = {scott_bandwidth}")
print (f "largura de banda CV = {cv_bandwidth}")

Como esperado, o primeiro Silverman e Scott retornam valores de largura de banda maiores que levam a caixas maiores e, portanto, a perda de informações sobre a distribuição de dados.

O Statsmodels permite uma busca automática da largura de banda ideal baseada na validação cruzada e um operador de máxima verossimilhança:

 de statsmodels.nonparametric.kernel_density import KDEMultivariate 
stats_models_cv = KDEMultivariate (dados, 'c', bw = 'cv_ml'). pdf (x)

Trace as diferentes aproximações

 figura.figura (figsize = (14, 6)) 
plt.plot (x, true_pdf, label = 'True PDF')
 kde = getKernelDensityEstimation (dados, x, largura de banda = silverman_bandwidth) 
plt.plot (x, kde, alpha = 0.8, label = largura de banda f'Silverman ')
 kde = getKernelDensityEstimation (dados, x, largura de banda = scott_bandwidth) 
plt.plot (x, kde, alfa = 0.8, label = largura de banda f'Scott ')
 kde = getKernelDensityEstimation (dados, x, largura de banda = cv_bandwidth) 
plt.plot (x, kde, alpha = 0.8, label = largura de banda f'CV ')
 plt.plot (x, stats_models_cv, alpha = 0,8, label = f'Matilidade máxima de CV dos estatísticos ') 
 plt.legend () 
plt.title ('Comparativo de várias estimativas de largura de banda para o KDE');
plt.savefig ('figures / comp_bw.jpg', bbox_inches = 'tight')

Testes estatísticos para distribuições unimodais

Existem vários testes estatísticos que abordam o problema da modalidade de dados:

  • Teste DIP
  • teste de excesso de massa
  • Teste MAP
  • teste de existência de modo
  • teste runt
  • teste de amplitude
  • teste de sela

Infelizmente, muitos não foram implementados em bibliotecas de código aberto em Python.

Teste DIP

O seguinte pacote python https://github.com/BenjaminDoran/unidip fornece uma implementação do teste de imersão e também uma funcionalidade para extrair de forma ecuosa os picos de densidade nos dados que utilizam o teste Hartigan Dip de Unimodality.

 de unidip import UniDip 
import unidip.dip como dip
 data = np.msort (data) 
imprimir (dip.diptst (data))
intervalos = UniDip (data) .run ()
imprimir (intervalos)

Identificar e plotar valores máximos locais do KDE

Uma vez que tenhamos uma estimativa da funç˜ao de densidade do núcleo, podemos determinar se a distribuiç˜ao é multimodal e identificar os valores máximos ou picos correspondentes aos modos.
Isso pode ser feito identificando os pontos em que a primeira derivada altera o sinal. O método getInflexion points pode retornar por padrão todos os pontos de inflexão (minima + maxima) ou apenas uma seleção (typeOfInflexion = 'max' / 'min').
O gráfico abaixo descreve esses valores máximos que podem corresponder a vários modos da distribuição de dados. A análise pode ser continuada definindo um limite baseado na altura do pico para filtrar alguns valores menos significativos.

 def getInflexionPoints (dados, typeOfInflexion = Nenhum, maxPoints = Nenhum): 
"" "
Este método retorna os índices onde há uma mudança na tendência da série de entrada.
typeOfInflexion = None retorna todos os pontos de inflexão, max apenas valores máximos e min
só min
"" "
a = np.diff (data)
asign = np.sign (a)
signchange = ((np.roll (asign, 1) - asign)! = 0) .asty (int)
idx = np.where (troca de sinal == 1) [0]
 if typeOfInflexion == 'max' e dados [idx [0]] <data [idx [1]]: 
idx = idx [1:] [:: 2]

elif typeOfInflexion == 'min' e dados [idx [0]]> data [idx [1]]:
idx = idx [1:] [:: 2]
elif typeOfInflexion não é nenhum:
idx = idx [:: 2]

# id de classificação por valor mínimo
se 0 em idx:
idx = np.delete (idx, 0)
if (len (data) -1) em idx:
idx = np.delete (idx, len (data) -1)
idx = idx [np.argsort (data [idx])]
# Se tivermos maxpoints, queremos ter certeza de que as timeseries têm um ponto de corte
# em cada segmento, nem todos em um pequeno intervalo
se maxPoints não for nenhum:
idx = idx [: maxPoints]
se len (idx) <maxPoints:
return (np.arange (maxPoints) + 1) * (len (dados) // (maxPoints + 1))

return idx

Percebemos que os valores obtidos correspondem às âncoras iniciais da distribuição gerada.

Recursos