Vectorizando dos muestras de prueba de Kolmogorov-Smirnov -- statistics campo con matlab campo con vectorization camp codereview Relacionados El problema

Vectorizing Two Sample Kolmogorov-Smirnov test


4
vote

problema

Español

Estoy implementando el prueba de dos muestras de Kolmogorov-Smirnov en Matlab. Debo admitir que conozco muy poco de estadísticas formales y simplemente estaba tratando de implementar la descripción en la página de Wikipedia.

En este momento, tomo como ingreso los dos vectores que se deben comparar x1 y x2 , y opcionalmente el valor p en el que se ejecutará la verificación. Produce 1 Si los vectores pasan la prueba y 9988777665544335 si no lo hacen.

Un par de preguntas que tengo:

¿Hay una mejor manera de calcular las funciones empíricas de distribución?

En este momento, evalué la función de distribución empírica en los puntos definidos en el vector t . En mi ingenua aproximación, lo hice ir desde el valor más pequeño de los vectores hasta el más grande, usando el doble de puntos que el vector con la mayoría de los puntos. No tengo idea, esta es la mejor manera de calcularlo, y sospechar que el vector no tiene que ser casi tan denso.

¿Puedo vectorizar la búsqueda a través del vector t ?

En este momento, tengo un bucle para el bucle que pasa por cada valor de t y calcula la función de distribución empírica de esa manera. Parece que puede haber una manera de vectorizar aún más esta operación y deshacerse de la bucle para el bucle. No he descubierto una buena manera de hacerlo sin embargo.

y, por supuesto, ¡cualquier otra sugerencia es bienvenida!

Mi código está abajo:

  function [OUT] = k_stest(x1,x2,p)  % Set default p-value     if nargin == 2          p = 0.05;     end   % Compute lenghts     N1 = length(x1);      N2 = length(x2);  % Set vector t over which the empirical distribution function will be computed.     t = linspace(min(min([x1 x2])),max(max([x1 x2])),2*max([N1 N2]));  % Initialize the statistic, negative values guarantees it will be overwritten by first call.     D = -1;   % Set c from tabulated values     if p == 0.10         c = 1.22;     elseif p == 0.05         c = 1.36;     elseif p == 0.025         c = 1.48;     elseif p == 0.01         c = 1.63;     elseif p == 0.005         c = 1.73;     elseif p == 0.001         c = 1.95;     else         disp('Invalid p-value. Only p = 0.10 0.05 0.025 0.01 0.005 0.001 are supported')         return     end  % Search though the vector t, computing the empirical distribution function for each t,  % and overwriting the statistic if a higher value is found.     for i = 1:length(t)         F1 = sum(x1<=t(i))/N1;         F2 = sum(x2<=t(i))/N2;         if abs(F2-F1) > D             D = abs(F2-F1);         end     end  % Compare the statistic to determine if the samples pass or fail.     if D == -1;        disp('Error, invalid input vectors')        return     elseif D > c*sqrt((N1+N2)/(N1*N2));         OUT = 0;     else         OUT = 1;     end  end   
Original en ingles

I'm implementing the Two-sample Kolmogorov-Smirnov test in MATLAB. I must admit I know very little of formal statistics and was simply trying to implement the description in wikipedia's page.

Right now I take as input the two vectors to be compared x1 and x2, and optionally the p-value that the check will run on. It outputs 1 if the vectors pass the test and 0 if they do not.

A couple of questions I have:

Is there a better way to compute the empirical distribution functions?

Right now I evaluate the empirical distribution function at the points defined in the vector t. In my naive approximation, I made it go from the smallest value of the vectors to the largest, using twice as many points as the vector with the most points. I have no idea is this is the best way to compute it, and suspect the vector doesn't have to be nearly as dense.

Can I vectorize the search through the t vector?

Right now I have a for loop that goes through every value of tand computes the empirical distribution function in that way. It seems that there might be a way to vectorize this operation further and get rid of the for loop. Haven't figured out a good way to do it however.

And of course, any other suggestions are welcome!

My code is below:

function [OUT] = k_stest(x1,x2,p)  % Set default p-value     if nargin == 2          p = 0.05;     end   % Compute lenghts     N1 = length(x1);      N2 = length(x2);  % Set vector t over which the empirical distribution function will be computed.     t = linspace(min(min([x1 x2])),max(max([x1 x2])),2*max([N1 N2]));  % Initialize the statistic, negative values guarantees it will be overwritten by first call.     D = -1;   % Set c from tabulated values     if p == 0.10         c = 1.22;     elseif p == 0.05         c = 1.36;     elseif p == 0.025         c = 1.48;     elseif p == 0.01         c = 1.63;     elseif p == 0.005         c = 1.73;     elseif p == 0.001         c = 1.95;     else         disp('Invalid p-value. Only p = 0.10 0.05 0.025 0.01 0.005 0.001 are supported')         return     end  % Search though the vector t, computing the empirical distribution function for each t,  % and overwriting the statistic if a higher value is found.     for i = 1:length(t)         F1 = sum(x1<=t(i))/N1;         F2 = sum(x2<=t(i))/N2;         if abs(F2-F1) > D             D = abs(F2-F1);         end     end  % Compare the statistic to determine if the samples pass or fail.     if D == -1;        disp('Error, invalid input vectors')        return     elseif D > c*sqrt((N1+N2)/(N1*N2));         OUT = 0;     else         OUT = 1;     end  end 
        

Lista de respuestas

2
 
vote
vote
La mejor respuesta
 

Tomaré esto en orden cronológico, desde la parte superior hacia abajo. Creo que he cubierto todos los aspectos del código =)

  function [OUT] = k_stest(x1,x2,p)   

Esto se ve bien, pero sería mejor si tuvieras espacios entre los argumentos de entrada. Además, k_stest no es un nombre muy bueno, ya que es difícil entender lo que solo está mirando el nombre. Leí en otra revisión que "no necesita los corchetes alrededor de la variable de salida, pero siempre debe incluirlo". No necesariamente estoy de acuerdo, generalmente los omití, pero eso es solo una cuestión de preferencias personales.


  % Set default p-value     if nargin == 2          p = 0.05;     end    

Esto está bien. Dependiendo de cómo llame a la función, puede agregar algunos cheques adicionales. Es p un escalar? Son x1 y x2 vectores? Esto no es necesario, pero puede ser útil si desea usar esta función durante mucho tiempo y no recordar los detalles. assert puede ser muy útil aquí. < / p>


  % Compute lenghts     N1 = length(x1);      N2 = length(x2);   

numel8 (Número de elementos) es mejor que length . Es más robusto y mucho más rápido para los largos vectores.


  k_stest0  

Esto es peligroso, y un poco desordenado.

k_stest1

asume 99887776655443312 y k_stest3 son vectores de filas, I.E. Vectores horizontales. Muchas funciones en Matlab devuelven vectores como vectores de columnas, es decir, vectores verticales de forma predeterminada, incluida la k_stest4 Operador (:).

Si desea que llame a esta función con dos vectores verticales en lugar de dos vectores horizontales, esto fallaría (debido a 99887776655443315 ).

en su lugar, debe hacer:

  k_stest6  

De esta manera está convirtiendo k_stest7 y k_stest8 a vertical, y concatenarlos verticalmente, antes de la llamada 99887766655443319 . Ahora solo necesita una llamada a % Set default p-value if nargin == 2 p = 0.05; end 0 y % Set default p-value if nargin == 2 p = 0.05; end 1 .


La forma en que está inicializando % Set default p-value if nargin == 2 p = 0.05; end 2 está bien!


  % Set default p-value     if nargin == 2          p = 0.05;     end  3  

Comparación de valores de puntos flotantes usando % Set default p-value if nargin == 2 p = 0.05; end 4 es peligroso. Por ejemplo, % Set default p-value if nargin == 2 p = 0.05; end 5 devolverá % Set default p-value if nargin == 2 p = 0.05; end 6 . En su lugar, debe comparar los valores con el valor deseado, pero con cierta tolerancia (por ejemplo, % Set default p-value if nargin == 2 p = 0.05; end 7 ).

  % Set default p-value     if nargin == 2          p = 0.05;     end  8  

Si dos valores están más cerca entre sí que % Set default p-value if nargin == 2 p = 0.05; end 9 entonces son esencialmente iguales.


  p0  

¡Tienes razón, la vectorización es la forma de ir! Use p1 Para crear una matriz lógica donde cada elemento que satisfaga la condición p2 son verdaderas y las otras son falsas. Luego suma toda la máscara, y toma la diferencia máxima.

  p3  

  p4  

La verificación de vectores de entrada no válidos debe estar en el inicio de la función. No puedo ver una razón por la que lo quieres al final en su lugar. Y por lo que puedo decir, nunca fallará (si la entrada no es válida, obtendrá un error más arriba).

en lugar del p5 y p66 , puede simplemente hacer:

  p7  

Para sumarlo todo, su función se puede escribir así. He omitido el comentario para la simplicidad. Debe incluirlos en su código.

  p8  
 

I'll take this in chronological order, from the top down. I think I have covered all aspects of the code =)

function [OUT] = k_stest(x1,x2,p) 

This looks nice, but it would be better if you had spaces between the input arguments. Also, k_stest is not a very good name, as it's hard to understand what it does just be looking at the name. I read in another review that "You don't need the brackets around the output variable, but you should always include it". I don't necessarily agree, I usually omit them, but that's just a matter of personal preferences.


% Set default p-value     if nargin == 2          p = 0.05;     end  

This is OK. Depending on how you call the function you might add some additional checks. Is p a scalar? Are x1 and x2 vectors? This is not necessary, but can be useful if you want to use this function a long time from now and don't remember the specifics. assert can be very useful here.


% Compute lenghts     N1 = length(x1);      N2 = length(x2); 

numel (number of elements) is better than length. It's more robust and a lot faster for long vectors.


% Set vector t over which the empirical distribution function will be computed.     t = linspace(min(min([x1 x2])),max(max([x1 x2])),2*max([N1 N2])); 

This is dangerous, and a bit messy.

[x1 x2] assumes both x1 and x2 are row-vectors, i.e. horizontal vectors. A lot of functions in Matlab returns vectors as column-vectors i.e. vertical vectors by default, including the colon operator (:).

If you happen to call this function with two vertical vectors instead of two horizontal vectors, this would fail (due to [x1 x2]).

Instead, you should do:

x12 = [x1(:); x2(:)];  t = linspace(min(x12), max(x12), 2*max([N1, N2])); 

This way you're converting x1 and x2 to vertical vectors, and concatenate them vertically, prior to the linspace call. Now you only need one call to max and min.


The way you're initializing D is OK!


% Set c from tabulated values     if p == 0.10         c = 1.22;     elseif p == 0.05         c = 1.36;     ... 

Comparing floating point values using == is dangerous. For instance 0.1 + 0.2 == 0.3 will return false. Instead, you should compare the values with the desired value, but with some tolerance (for instance eps).

    if abs(p - 0.10) < eps         c = 1.22;      elseif abs(p-0.05) < eps         c = 1.36; 

If two values are closer to each other than +/- eps then they are essentially equal.


% Search though the vector t, computing the empirical distribution function for each t,  % and overwriting the statistic if a higher value is found.     for i = 1:length(t)         F1 = sum(x1<=t(i))/N1;         F2 = sum(x2<=t(i))/N2;         if abs(F2-F1) > D             D = abs(F2-F1);         end     end 

You are right, vectorization is the way to go! Use bsxfun to create a logical matrix where each element that satisfy the condition x1 <= t(ii) are true and the others are false. Then sum the entire mask, and take the maximum difference.

mask1 = bsxfun(@le, x1(:), t); mask2 = bsxfun(@le, x2(:), t); s1 = sum(mask1) / N1; s2 = sum(mask2) / N2; D = max(abs(s1 - s2)); 

% Compare the statistic to determine if the samples pass or fail.     if D == -1;        disp('Error, invalid input vectors')        return     elseif D > c*sqrt((N1+N2)/(N1*N2));         OUT = 0;     else         OUT = 1;     end     end 

The check for invalid input vectors should be in the start of the function. I can't see a reason why you want it in the end instead. And as far as I can tell it will never fail (if the input is invalid then you will get an error further up).

Instead of the if and elseif, you can simply do:

OUT = D <= c*sqrt((N1+N2)/(N1*N2)); 

To sum it all up, your function can be written like this. I've omitted the comment for simplicity. You should include them in your code.

function OUT = k_stest(x1, x2, p)  if nargin == 2     p = 0.05; end  % Add check for invalid input vectors!  N1 = numel(x1); N2 = numel(x2);  x12 = [x1(:); x2(:)];  t = linspace(min(x12), max(x12), 2*max([N1, N2]));  if abs(p - 0.10) < eps     c = 1.22;  elseif abs(p-0.05) < eps     c = 1.36; % Continue with the rest of the values % end  mask1 = bsxfun(@le, x1(:), t); mask2 = bsxfun(@le, x2(:), t); s1 = sum(mask1) / N1; s2 = sum(mask2) / N2; D = max(abs(s1 - s2));  OUT = D <= c*sqrt((N1+N2)/(N1*N2));  end 
 
 
   
   

Relacionados problema

6  ¿Cómo hacer que mis operaciones de Groupby y Transponen eficientes?  ( How to make my groupby and transpose operations efficient ) 
Tengo un DataFrame 3,745,802 filas y 30 columnas. Me gustaría realizar cierto 9988776665544331 y 9988776655544332 terminará finalmente con más columnas ...

3  Generar vectores de la unidad aleatoria alrededor del círculo  ( Generate random unit vectors around circle ) 
Estoy tratando de generar un montón de vectores de unidad de distribución uniformemente alrededor del círculo de la unidad. Aquí está mi código, que está func...

4  Vectorizando un bucle de conversión de coordenadas polares en un proceso de coincidencia de huellas dactilares  ( Vectorizing a polar coordinate conversion loop in a fingerprint matching process ) 
Estoy trabajando en una técnica de coincidencia de huellas dactilares y se necesita mucho tiempo de computación para obtener el resultado. Estoy tratando de i...

6  Cálculo de estacionamientos y estadísticas de carga para un servicio compartido de automóviles  ( Calculating parking and charging statistics for a car sharing service ) 
Estoy buscando ayuda para diseñar mi código de manera más eficiente y también en mantener la pendiente de mi curva de aprendizaje personal empinada: Para un...

2  Encuentra el primer umbral de cruce en la matriz NUTPY  ( Find first threshold crossing in numpy array ) 
Dado los siguientes datos generados artificialmente: t_steps = 30 data = np.array([ np.arange(t_steps) * .05, np.arange(t_steps) * .1, np.aran...

1  Vectorizamos gradiente 2D con contenedores espacialmente variables  ( Vectorize 2d gradient with spatially varying bins ) 
El siguiente código lleva a algunos valores ssh y resuelve las ecuaciones para el movimiento geostrófico (diapositiva 8 aquí ). La parte principal del Có...

2  Escalado e incremento de elementos no cero de una matriz adorable  ( Scaling and incrementing non zero elements of a numpy matrix ) 
Tengo una matriz adorable C y desea crear una copia cPrime , que tiene algún funcionamiento de la matriz original a todos los valores no cero. En el código...

4  Crosstabulación vectorizada en Python para dos matrices con dos categorías cada una  ( Vectorized crosstabulation in python for two arrays with two categories each ) 
Tengo dos listas de Python, label y presence1 . Quiero hacer la tabulación cruzada y obtener el recuento para cada bloque de cuatro, como A, B, C y D en ...

5  Vectorización para simulación de temperatura  ( Vectorization for temperature simulation ) 
Soy nuevo en Matlab y me gustaría saber consejos para reducir el tiempo de procesamiento de dicho programa. El esquema del código real es muy similar al códig...

5  Evaluando una serie de polinomios legendra  ( Evaluating a series of legendre polynomials ) 
La siguiente función representa el potencial electrostático, en coordenadas esféricas, debido a un anillo de carga $ Q = 1 $ y RADIUS $ R = 1 $, colocado ...




© 2022 respuesta.top Reservados todos los derechos. Centro de preguntas y respuestas reservados todos los derechos