Komenda equavelent python dla kwantyla w matlab

Question

Próbuję replikować niektóre z mojego kodu w Matlab do Pythona. Znalazłem funkcję kwantyle w Matlab nie ma "dokładnie" odpowiadające w python. To, co znalazłem najbliżej, to mquantiles Pythona. na przykład

dla Matlab:

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

daje: 0.00016958

dla Pythona:

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

daje 0.00016912

Czy ktoś wie, jak dokładnie replikować kwantyl Matlaba? wielkie dzięki.

Answer 1

Twój wektor wejściowy ma tylko 4 wartości, co stanowi o wiele za mało, aby uzyskać dobre przybliżenie kwantyli leżących u podstaw dystrybucji. Rozbieżność jest prawdopodobnie wynikiem Matlab i SciPy przy użyciu różnych heurystyk do obliczania kwantyli w próbkowanych dystrybucjach.

Answer 2

W documentation for quantile (w sekcji Więcej o = = Algorytmy) podano dokładny algorytm. Oto niektóre kodu Pythona, który robi to za pomocą pojedynczego kwantyl na płaskiej tablicy, używając bottleneck zrobić częściową Sortowanie:

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) 

ja tylko zrobiłem jednego kwantyli, 1d sprawy, ponieważ to wszystko, co potrzebne. Jeśli chcesz mieć kilka kwantyli, prawdopodobnie warto po prostu dokonać pełnego sortowania; zrobić to na oś i wiedziałem, że nie masz żadnych nans, wszystko co musisz zrobić, to dodać argument osi do sortowania i wektoryzować bit interpolacji liniowej. Wykonanie tego na osi z nans byłoby nieco trudniejsze.

Ten kod daje:

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

i kod MATLAB dał 0.00016905822359999999; różnica wynosi 3e-20. (Która jest mniejsza niż precyzyjnej)

Answer 3

Nieco późno, ale:

mquantiles jest bardzo elastyczny. Musisz tylko podać parametry alphap i betap. Tutaj, ponieważ MATLAB wykonuje interpolację liniową, musisz ustawić parametry na (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) 

Edycja: MATLAB mówi, że wykonuje interpolację liniową, ale wydaje się, że oblicza kwantylu przez kawałek mądry interpolacji liniowej, który jest równoważny typ 5 kwantylem w R i (0,5, 0,5) w scipy.