Comando equavelent python per quantile in MATLAB

Question

Sto provando a replicare parte del mio codice in MATLAB a Python. Ho trovato che la funzione di quantili in MATLAB non ha uno "esattamente" corrispondente in Python. Quello che ho trovato più vicino sono i mQuantile di Python. per esempio.

per MATLAB:

quantile([ 8.60789925e-05, 1.98989354e-05 , 1.68308882e-04, 1.69379370e-04], 0.8) 

dà: 0.00016958

per Python:

scipy.stats.mstats.mquantiles([8.60789925e-05, 1.98989354e-05, 1.68308882e-04, 1.69379370e-04], 0.8) 

dà 0.00016912

Qualcuno sa come replicare esattamente quantile di MATLAB? molte grazie.

Answer 1

Il vettore di input ha solo 4 valori, che sono troppo pochi per ottenere una buona approssimazione dei quantili della distribuzione sottostante. La discrepanza è probabilmente il risultato di Matlab e SciPy che utilizzano diverse euristiche per calcolare quantili su distribuzioni campionate.

Answer 2

Il documentation for quantile (nella sezione More About => Algorithms) fornisce l'algoritmo esatto utilizzato. Ecco qualche codice Python che fa per un singolo quantile per una matrice piatta, utilizzando bottleneck fare ordinamento parziale:

import numpy as np 
import botteleneck as bn 

def quantile(a, prob): 
    """ 
    Estimates the prob'th quantile of the values in a data array. 

    Uses the algorithm of matlab's quantile(), namely: 
     - Remove any nan values 
     - Take the sorted data as the (.5/n), (1.5/n), ..., (1-.5/n) quantiles. 
     - Use linear interpolation for values between (.5/n) and (1 - .5/n). 
     - Use the minimum or maximum for quantiles outside that range. 

    See also: scipy.stats.mstats.mquantiles 
    """ 
    a = np.asanyarray(a) 
    a = a[np.logical_not(np.isnan(a))].ravel() 
    n = a.size 

    if prob >= 1 - .5/n: 
     return a.max() 
    elif prob <= .5/n: 
     return a.min() 

    # find the two bounds we're interpreting between: 
    # that is, find i such that (i+.5)/n <= prob <= (i+1.5)/n 
    t = n * prob - .5 
    i = np.floor(t) 

    # partial sort so that the ith element is at position i, with bigger ones 
    # to the right and smaller to the left 
    a = bn.partsort(a, i) 

    if i == t: # did we luck out and get an integer index? 
     return a[i] 
    else: 
     # we'll linearly interpolate between this and the next index 
     smaller = a[i] 
     larger = a[i+1:].min() 
     if np.isinf(smaller): 
      return smaller # avoid inf - inf 
     return smaller + (larger - smaller) * (t - i) 

ho fatto solo la singola-quantile, caso 1d perché è tutto il necessario. Se vuoi diversi quantili, probabilmente vale la pena fare l'ordinamento completo; per farlo per asse e sapeva che non avevi nessun nans, tutto ciò che dovresti fare è aggiungere un argomento sugli assi all'ordinamento e vettorizzare il bit di interpolazione lineare. Facendolo per asse con nans sarebbe un po 'più complicato.

Questo codice dà:

>>> quantile([ 8.60789925e-05, 1.98989354e-05 , 1.68308882e-04, 1.69379370e-04], 0.8) 
0.00016905822360000001 

e il codice MATLAB ha dato 0.00016905822359999999; la differenza è 3e-20. (Che è meno di precisione della macchina)

Answer 3

un po 'tardi, ma:

mquantiles è molto flessibile. Hai solo bisogno di fornire parametri alphap e betap. Qui, poiché MATLAB esegue un'interpolazione lineare, è necessario impostare i parametri su (0,5,0.5).

In [9]: scipy.stats.mstats.mquantiles([8.60789925e-05, 1.98989354e-05, 1.68308882e-04, 1.69379370e-04], 0.8, alphap=0.5, betap=0.5) 

EDIT: MATLAB dice che fa interpolazione lineare, tuttavia sembra che calcola il quantile attraverso piece-wise interpolazione lineare, che è equivalente al Tipo 5 quantile in R, e (0.5, 0.5) in scipy.