function [spurt_data] = ...
                get_spurt_data_bootstrap( t0, y0, max_min_data, fdobj0, degree, ...
                                          boot_samples, lambda_fct, lambda_stderr, ...
                                          error_level, pvalue )

    % This function returns mean and std error information for the maxima in 
    % max_min_data using a numerical (bootstrap) method.
    %
    % See Dandurand & Shultz paper on AMD for details
    %
    % Input parameters
    % ----------------
    % t0 and y0 - original raw data 
    % expected_spurts_count - how many spurts this curve contains
    % max_min_data - structure containing info on locations of x values where to find maxima and mimima 
    %                3 fields:
    %                * location: where extrenum is located (t value)
    %                * type: what kind of extrenum (see global_constants for details)
    %                * significant: if the extrenum is statistically
    %                significant
    %                We will use newton optimization method starting from those points to get
    %                the new ones: location, x1 and x2
    % fdobj0 - the FDA fit of the raw data
    % degree - in which derivative we are looking for max and mins
    % boot_samples - how many iterations to do (the more, the better the
    %               estimate of error, but the more computation it takes)
    % lambda_fct    smoothing parameter for the function curve
    % lambda_stderr smoothing parameter for the std error curve
    % error_level   required precision for zero crossing estimates 
    %               that is: abs(F(estimate_xvalue)) < error_level
    % 
    % Output value
    % ------------
    % spurt_data: a structure with the following fields:
    % * begin_location - location of the beginning of the spurt
    % * center_location - Spurt center location (x value)
    % * stderr_location - Standard error on the location estimate
    % * amplitude - Spurt amplitude (or intensity)
    % * stderr_amplitude - Standard error on the amplitude estimate 
    % * duration - Duration of the spurt, in x axis units (e.g., time)
    % * stderr_duration - Standard error on the duration estimate
    % (One entry per maximum (spurt))
    
    % TO DO: output the proportion of samples that didn't have the same
    % number of spurts (skipped samples). This may be indicative of the 
    % sensitivity of the analysis.

    global_constants;
    
    MAX_RETRY= 5*boot_samples;

    est0 = eval_fd(t0, fdobj0, 0);  % FDA estimates of the original curve
    r0   = y0 - est0;                  % residuals of the original fit
    
    last_boot_index = 1;
    
    orig_spurt_data  = ...
        get_spurt_data_means_only(max_min_data, fdobj0, degree, false);
    
    expected_spurts_count = size(orig_spurt_data,2);
    
    begin_data = zeros(expected_spurts_count, boot_samples);
    center_data = zeros(expected_spurts_count, boot_samples);
    amplitude_data = zeros(expected_spurts_count, boot_samples);
    duration_data = zeros(expected_spurts_count, boot_samples);
    
    retry = 0;
    
    while(last_boot_index <= boot_samples)
        % generate a new vector of observations based on est0, but with random
        % noise estimated using the residuals
        %
        % randn = generates random values from normal distributions
        % From Matlab help: Generate a random distribution with a specific mean 
        % and variance (sigma2). To do this, multiply the output of randn by 
        % the standard deviation (sigma), and then add the desired mean
        new_t = t0;

        % Because e0 should be very small (close to zero), we don't
        % include it in the computation of new_y.
        resampleIndex = randi(size(r0,1), size(r0,1), 1);
        new_y = est0 + r0(resampleIndex);

        % generate a new FDA curve to model these samples
        [ fdobj stderr y2cMap ] = ...
            create_fd_objects( new_t, new_y, lambda_fct, lambda_stderr, degree );

        % find the locations (i.e., x values) where maxima and minima in the
        % degree'th and degree+1'th derivatives of the fitted function (e.g., if
        % degree = 0, find max and mins in the function and the velocity)
        max_min_mat_loc = find_min_max_locations(new_t, fdobj, degree, error_level);
        x_values        = [max_min_mat_loc.location]';   % extract locations of the maxs and mins

        % compute the standard error for the degree'th and degree+1'th derivatives
        stderr_degree = get_stderr_estimates(x_values, fdobj, stderr, y2cMap, degree);

        % determine which spurts are significant
        significant_spurt_entry = determine_spurts_significance(max_min_mat_loc, pvalue, ...
            fdobj, stderr_degree, degree);

        % we do not require these spurts to be significant in the resampled version, so use "true"
        spurt_data_subject  = ...
            get_spurt_data_means_only(significant_spurt_entry, fdobj, degree, true);
        
        spurts_count = size(spurt_data_subject,2);
        
        if (spurts_count == expected_spurts_count)
            % process spurts for this curve. There can be made spurts for a single
            for spurt_id=1:spurts_count
                begin_data(spurt_id, last_boot_index) = spurt_data_subject(spurt_id).begin_location;
                center_data(spurt_id, last_boot_index) = spurt_data_subject(spurt_id).center_location;
                amplitude_data(spurt_id, last_boot_index) = spurt_data_subject(spurt_id).amplitude;
                duration_data(spurt_id, last_boot_index) = spurt_data_subject(spurt_id).duration;
            end
            last_boot_index = last_boot_index + 1;
        else
            % skip
            retry = retry +1;
        end
        
        if (retry > MAX_RETRY)
            disp('Unable to compute variance')
            break
        end
    end

    spurt_data = orig_spurt_data;
    
    for spurt_id = 1:expected_spurts_count
        begin_var         = var(begin_data(spurt_id, :));
        begin_stderr      = sqrt(begin_var / boot_samples);
        spurt_data(spurt_id).stderr_begin     = begin_stderr;

        center_var         = var(center_data(spurt_id, :));
        center_stderr      = sqrt(center_var / boot_samples);
        spurt_data(spurt_id).stderr_center    = center_stderr;

        amp_var         = var(amplitude_data(spurt_id, :));
        amp_stderr      = sqrt(amp_var / boot_samples);
        spurt_data(spurt_id).stderr_amplitude = amp_stderr;

        dur_var         = var(duration_data(spurt_id, :));
        dur_stderr      = sqrt(dur_var / boot_samples);
        spurt_data(spurt_id).stderr_duration  = dur_stderr;
    end
return;
