Crear un conjunto de datos de intervalo fijo a partir de un conjunto de datos de intervalo aleatorio con datos obsoletos

Question

Actualización: He proporcionado un breve análisis de las tres respuestas en la parte inferior del texto de la pregunta y he explicado mis elecciones.

Mi pregunta: ¿Cuál es el método más eficiente de construir un conjunto de datos de intervalo fijo a partir de un conjunto de datos de intervalo aleatorio que utiliza datos obsoletos?

Algunos antecedentes: Lo anterior es un problema común en las estadísticas. Con frecuencia, uno tiene una secuencia de observaciones que ocurren en momentos aleatorios. Llámelo Input. Pero uno quiere que ocurra una secuencia de observaciones, por ejemplo, cada 5 minutos. Llámalo Output. Uno de los métodos más comunes para construir este conjunto de datos es el uso de datos obsoletos, es decir, establecer cada observación en Output igual a la observación más reciente en Input.

lo tanto, aquí hay un código de ejemplo para construir conjuntos de datos:

TInput = 100; 
TOutput = 50; 

InputTimeStamp = 730486 + cumsum(0.001 * rand(TInput, 1)); 
Input = [InputTimeStamp, randn(TInput, 1)]; 

OutputTimeStamp = 730486.002 + (0:0.001:TOutput * 0.001 - 0.001)'; 
Output = [OutputTimeStamp, NaN(TOutput, 1)]; 

Ambos conjuntos de datos comienzan en cerca de la medianoche en el cambio de milenio. Sin embargo, las marcas de tiempo en Input se producen a intervalos aleatorios, mientras que las marcas de tiempo en Output se producen a intervalos fijos. Por simplicidad, he asegurado que la primera observación en Input siempre se produce antes de la primera observación en Output. Siéntase libre de hacer esta suposición en cualquier respuesta.

Actualmente, se soluciona el problema como este:

sMax = size(Output, 1); 
tMax = size(Input, 1); 
s = 1; 
t = 2; 
%#Loop over input data 
while t <= tMax 
    if Input(t, 1) > Output(s, 1) 
     %#If current obs in Input occurs after current obs in output then set current obs in output equal to previous obs in input 
     Output(s, 2:end) = Input(t-1, 2:end); 
     s = s + 1; 
     %#Check if we've filled out all observations in output 
     if s > sMax 
      break 
     end 
     %#This step is necessary in case we need to use the same input observation twice in a row 
     t = t - 1; 
    end 
    t = t + 1; 
    if t > tMax 
     %#If all remaining observations in output occur after last observation in input, then use last obs in input for all remaining obs in output 
     Output(s:end, 2:end) = Input(end, 2:end); 
     break 
    end 
end 

Seguramente hay una manera más eficiente, o por lo menos, forma más elegante de resolver este problema? Como mencioné, este es un problema común en las estadísticas. ¿Quizás Matlab tiene alguna función incorporada de la que no estoy al tanto? Cualquier ayuda sería muy apreciada ya que uso esta rutina MUCHO para algunos conjuntos de datos grandes.

LAS RESPUESTAS: Hola a todos, he analizado las tres respuestas, y tal como están, Angainor es la mejor.

La respuesta de ChthonicDaemon, si bien es claramente la más fácil de implementar, es realmente lenta. Esto es cierto incluso cuando la conversión a un objeto timeseries se realiza fuera de la prueba de velocidad. Supongo que la función resample tiene muchos sobrecarga en este momento. Estoy ejecutando 2011b, por lo que es posible que Mathworks lo haya mejorado en el tiempo transcurrido. Además, este método necesita una línea adicional para el caso en que Output extremos más de una observación después de Input.

La respuesta de Rody es solo un poco más lenta que la de Angainor (no es sorprendente dado que ambos emplean el enfoque histc), sin embargo, parece tener algunos problemas. Primero, el método de asignar la última observación en Output no es robusto a la última observación en Input que ocurre después de la última observación en Output. Esta es una solución fácil. Pero hay un segundo problema que creo proviene de tener InputTimeStamp como la primera entrada a histc en lugar del OutputTimeStamp adoptado por Angainor. El problema surge si cambia OutputTimeStamp = 730486.002 + (0:0.001:TOutput * 0.001 - 0.001)'; a OutputTimeStamp = 730486.002 + (0:0.0001:TOutput * 0.0001 - 0.0001)'; al configurar las entradas de ejemplo.

Angainor's parece robusto para todo lo que le lancé, además fue el más rápido.

Hice un montón de pruebas de velocidad para diferentes especificaciones de entrada - los siguientes números son bastante representativas:

Mi bucle ingenua: Elapsed time is 8.579535 seconds.

Angainor: Elapsed time is 0.661756 seconds.

Rody: Elapsed time is 0.913304 seconds.

ChthonicDaemon: Elapsed time is 22.916844 seconds.

Estoy + 1-ing la solución de Angainor y estoy marcando la cuestión resuelta.

Answer 1

Aquí está mi opinión sobre el problema. histc es el camino a seguir:

% find Output timestamps in Input bins 
N = histc(Output(:,1), Input(:,1)); 

% find counts in the non-empty bins 
counts = N(find(N)); 

% find Input signal value associated with every bin 
val = Input(find(N),2); 

% now, replicate every entry entry in val 
% as many times as specified in counts 
index = zeros(1,sum(counts)); 
index(cumsum([1 counts(1:end-1)'])) = 1; 
index = cumsum(index); 
val_rep = val(index) 

% finish the signal with last entry from Input, as needed 
val_rep(end+1:size(Output,1)) = Input(end,2); 

% done 
Output(:,2) = val_rep; 

he comprobado en contra de su procedimiento durante unos modelos de entrada diferentes (I cambió el número de marcas de tiempo de salida) y los resultados son los mismos. Sin embargo, todavía no estoy seguro de haber entendido tu problema, así que si algo está mal, házmelo saber.

Answer 2

Este enfoque de "datos obsoletos" se conoce como zero order hold en los campos de señal y horario. La búsqueda de esto rápidamente trae muchas soluciones. Si tiene Matlab 2012b, todo esto está integrado en la clase timeseries mediante el uso de la función resample, por lo que sólo tendría que hacer

TInput = 100; 
TOutput = 50; 

InputTimeStamp = 730486 + cumsum(0.001 * rand(TInput, 1)); 
InputData = randn(TInput, 1); 
InputTimeSeries = timeseries(InputData, InputTimeStamp); 

OutputTimeStamp = 730486.002 + (0:0.001:TOutput * 0.001 - 0.001); 
OutputTimeSeries = resample(InputTimeSeries, OutputTimeStamp, 'zoh'); % zoh stands for zero order hold