function [hFig, hAxes, hCurve] = accordPlotMaker(hAxes, fileToLoad,...
    passiveID, molID, customObsProp, customCurveProp, customFigProp, customAxesProp)
%
% The AcCoRD Simulator
% (Actor-based Communication via Reaction-Diffusion)
%
% Copyright 2016 Adam Noel. All rights reserved.
% 
% For license details, read LICENSE.txt in the root AcCoRD directory
% For user documentation, read README.txt in the root AcCoRD directory
%
% accordPlotMaker.m - draw simulation data as curves or surfaces. Options
%   include time-varying samples, probability mass functions (PMFs), and
%   cumulative
%   distribution functions (CDFs). Any or all realizations in the imported
%   simulation output file can be included. This function can be called
%   repeatedly to draw multiple curves on one axes. 2D plots are made using
%   the plot3 function (so that lines can be overlaid over 3D surfaces). 3D
%   plots are made use the surf function.
%
% INPUTS
% hAxes - axes to plot data. Should be "0" if figure does not yet exist
% fileToLoad - simulation file generated by accordImport. Ignored for
%   curves or surfaces that are some type "Expected"
% passiveID - ID of passive actor whose simulation data are being plotted.
%   Index is from recorded passive actor list.
% molID - ID of molecule observed by passive actor. Index is from list of
%   molecules recorded by passiveID
% customObsProp - structure of observer properties to change from AcCoRD
%   defaults. Can be passed as empty if no defaults are to be changed. See
%   accordBuildObserverStruct for structure fields and their default values.
% customCurveProp - structure of curve/surface properties to change from
%   AcCoRD defaults. Can be passed as empty if no defaults are to be
%   changed. Possible members depend on value of obsSpec.obsType. If plot
%   is 2D, then see accordBuildCurveStruct for structure fields and their
%   default values. If plot is 3D, then see accordBuildSurfStruct for
%   structure fields and their default values.
% customFigProp - structure of figure properties to change from AcCoRD
%   defaults. Can be passed as empty if no defaults are to be changed. See
%   accordBuildFigureStruct for structure fields and their default values.
% customAxesProp - structure of axes properties to change from AcCoRD
%   defaults. Can be passed as empty if no defaults are to be changed. See
%   accordBuildAxesStruct for structure fields and their default values.
%
% OUTPUTS
% hFig - handle(s) to plotted figure(s). Use for making changes.
% hAxes - handle(s) to axes in plotted figure(s). Use for making changes.
%
% Last revised for AcCoRD v1.2 (2018-05-30)
%
% Revision history:
%
% Revision v1.2 (2018-05-30)
% - corrected the range generated for 3D Empirical PMF (was off by 1)
%
% Revision v1.0 (2016-10-31)
% - changed mutual information calculations to use log2 (i.e., bits) instead of
% log (i.e., nats)
%
% Revision v0.7.0.1 (public beta, 2016-08-30)
% - added call to addpath to make JSON code visible if not already added
%
% Revision v0.7 (public beta, 2016-07-09)
% - Created file
%
% Created 2016-06-03


% Add subdirectory with JSONlab files to path (needed for loadjson)
addpath('JSONlab');

%% Load Default Display Properties and Apply Specified Changes
obsSpec = accordBuildObserverStruct(customObsProp);

%% Create figure and apply properties if it does not exist
if hAxes == 0
    figureProp = accordBuildFigureStruct(customFigProp);
    axesProp = accordBuildAxesStruct(customAxesProp);
    
    %% Create figure and apply properties
    hFig = figure;
    if ~isempty(figureProp)
        figurePropFields = fieldnames(figureProp);
        numProp = numel(figurePropFields);
        for i = 1:numProp
            set(hFig, figurePropFields{i}, figureProp.(figurePropFields{i}));
        end
    end

    %% Create axes and apply properties
    hAxes = axes('Parent',hFig);
    if ~isempty(axesProp)
        axesPropFields = fieldnames(axesProp);
        numProp = numel(axesPropFields);
        for i = 1:numProp
            set(hAxes, axesPropFields{i}, axesProp.(axesPropFields{i}));
        end
    end
    hold(hAxes, 'on');
    
    % Set axes labels based on type of plot
    switch obsSpec.obsType
        case {'Sample','Expected'}
            xlabel(hAxes, 'Time [s]');
            ylabel(hAxes, 'Number of molecules');
        case '3D Expected'
            % No obvious labels to use in this case
        case {'Empirical CDF', 'Expected CDF'}
            xlabel(hAxes, 'Number of Molecules');
            ylabel(hAxes, 'Cumulative Distribution Function');
        case {'3D Empirical CDF', '3D Expected CDF'}
            xlabel(hAxes, 'Time [s]');
            ylabel(hAxes, 'Number of Molecules');
            zlabel(hAxes, 'Cumulative Distribution Function');
        case {'Empirical PMF', 'Expected PMF'}
            xlabel(hAxes, 'Number of Molecules');
            ylabel(hAxes, 'Frequency');
        case {'3D Empirical PMF', '3D Expected PMF'}
            xlabel(hAxes, 'Time [s]');
            ylabel(hAxes, 'Number of Molecules');
            zlabel(hAxes, 'Frequency');
        case 'Mutual Information'
            xlabel(hAxes, 'Time Offset [s]');
            ylabel(hAxes, 'Mutual Information');
        case '3D Mutual Information'
            xlabel(hAxes, 'Reference Time [s]');
            ylabel(hAxes, 'Time Offset [s]');
            zlabel(hAxes, 'Mutual Information');
        case 'Monte Carlo Mutual Information'
            xlabel(hAxes, 'Time Offset [s]');
            ylabel(hAxes, 'Mutual Information');
        case '3D Monte Carlo Mutual Information'
            xlabel(hAxes, 'Reference Time [s]');
            ylabel(hAxes, 'Time Offset [s]');
            zlabel(hAxes, 'Mutual Information');
    end
else
    hFig = get(hAxes, 'Parent');
end

%% Is plot 3D? If so then load surface properties instead of line properties
b3D = false;
switch obsSpec.obsType
    case {'3D Expected', ...
            '3D Empirical PMF', '3D Empirical CDF', ...
            '3D Expected PMF', '3D Expected CDF', ...
            '3D Mutual Information', '3D Monte Carlo Mutual Information'}
        b3D = true;
        curveSpec = accordBuildSurfStruct(customCurveProp);
    otherwise
        curveSpec = accordBuildCurveStruct(customCurveProp);
end


%% Load Data
switch obsSpec.obsType
    case {'Sample', 'Empirical CDF', '3D Empirical CDF', 'Empirical PMF', ...
            '3D Empirical PMF', 'Mutual Information', '3D Mutual Information'}
        % Plotting actual simulation data. Take samples from the data file
        
        load(fileToLoad, 'config', 'data');

        obsMatrix = data.passiveRecordCount{passiveID}(:,molID,:);
        numObs = size(obsMatrix,3);
        numRepeat = size(obsMatrix,1);

        % Determine all original simulation observation times
        actorID = 1;
        while config.actor{actorID}.passiveID ~= passiveID
            actorID = actorID + 1;
        end
        tArrayFull = config.actor{actorID}.startTime + ...
            (0:(numObs-1))*config.actor{actorID}.actionInterval;

        % Determine times that will be sampled
        if obsSpec.bMaxSample
            tPlotInd = obsSpec.firstSample:obsSpec.sampleInterval:numObs;
        else
            tPlotInd = obsSpec.firstSample:obsSpec.sampleInterval:obsSpec.lastSample;
        end
        
        % Determine what repeats to use (default is all)
        switch obsSpec.iterType
            case 'All'
                iter = 1:numRepeat;
            case 'Custom'
                iter = obsSpec.iterCustom;
            otherwise
        end
        
        % Extract the sampled observations
        obsMatrixSampled = obsMatrix(iter,:,tPlotInd);
        numPlotInd = length(tPlotInd);
    otherwise
end

%% Determine data to plot
switch obsSpec.obsType
    case 'Sample'
        % Plotting (average) time-varying observations
        xData = tArrayFull(tPlotInd);
        yData = reshape(mean(obsMatrixSampled(:,:),1),1,[]);
    case 'Expected'
        % Plotting specified data
        xData = obsSpec.data1;
        yData = obsSpec.data2;
    case '3D Expected'
        % Plotting specified data
        xData = obsSpec.data1;
        yData = obsSpec.data2;
        zData = obsSpec.data3;
    case 'Empirical CDF'
        % Plotting empirical cumulative distribution function over all
        % sampled observations
        [yData, xData] = ecdf(obsMatrixSampled(:));
    case '3D Empirical CDF'
        % Plotting time-varying empirical cumulative distribution function.
        % Every sampled time has its own CDF
        numYMax = max(obsMatrixSampled(:))+1;
        
        xData = zeros(numPlotInd,numYMax);
        yData = zeros(numPlotInd,numYMax);
        zData = zeros(numPlotInd,numYMax);
        
        for i = 1:numPlotInd
            xData(i,:) = tArrayFull(tPlotInd(i));
            [zDataCur, yDataCur] = ecdf(obsMatrixSampled(:,:,i));
            numY = length(yDataCur);
            if yDataCur > 0
                % Force first observation to be 0, even if it isn't
                % observed
                firstPoint = 2;
                lastPoint = 1+numY;
            else
                firstPoint = 1;
                lastPoint = numY;
            end
            yData(i,firstPoint:lastPoint) = yDataCur;
            zData(i,firstPoint:lastPoint) = zDataCur;
            % "Extend" display of CDF to global maximum observation (so
            % that CDF appears square from "above")
            yData(i,(lastPoint+1):end) = numYMax;
            zData(i,(lastPoint+1):end) = zDataCur(end);
        end
    case 'Expected CDF'
        % Plotting expected cumulative distribution function based on
        % chosen distribution. If data1 is non-empty, then it is used to
        % define the bin values. Otherwise, we center the bins around the
        % expected mean
        meanExp = obsSpec.trialProbability*obsSpec.numTrials;
        varExp = obsSpec.trialProbability*obsSpec.numTrials*(1-obsSpec.trialProbability);
        if isempty(obsSpec.data1)
            xData = linspace(meanExp-4*sqrt(varExp),meanExp+4*sqrt(varExp),obsSpec.numMeasurements);
        else
            xData = obsSpec.data1;
        end
        xData = unique(floor(xData(xData>=0))); % Ignore negative and repeated elements
        
        yData = zeros(1,length(xData));
        switch obsSpec.distribution
            case 'Binomial'
                yData = betainc(1-obsSpec.trialProbability,obsSpec.numTrials-xData, xData+1);
            case 'Poisson'
                yData = gammainc(meanExp, floor(xData+1), 'upper');
            case 'Gaussian'
                yData = 0.5 * (1 + erf((xData+0.5-meanExp)./sqrt(2*varExp)));
        end
    case '3D Expected CDF'
        % Plotting expected cumulative distribution function based on
        % chosen distribution. If data1 is non-empty, then it is used to
        % define the bin values. Otherwise, we center the bins around the
        % maximum expected mean
        % data2 vector is used to define the time-values, so its length
        % should equal that of trialProbability
        numTimeInd = length(obsSpec.trialProbability);
        if length(obsSpec.data2) ~= numTimeInd
            error('Length of observer data2 vector does not equal number of trial probabilities');
        end
        meanExp = obsSpec.trialProbability*obsSpec.numTrials;
        varExp = obsSpec.trialProbability*obsSpec.numTrials.*(1-obsSpec.trialProbability);
        
        if isempty(obsSpec.data1)
            binVals = linspace(min(meanExp)-4*sqrt(min(varExp)),max(meanExp)+4*sqrt(max(varExp)),obsSpec.numMeasurements);
        else
            binVals = obsSpec.data1;
        end
        binVals = unique(floor(binVals(binVals>=0))); % Ignore negative and repeated elements
        
        numBin = length(binVals);
        
        xData = repmat(obsSpec.data2',1,numBin);
        yData = repmat(binVals,numTimeInd,1);
        zData = zeros(numTimeInd, numBin);
        
        for j = 1:numTimeInd
            switch obsSpec.distribution
                case 'Binomial'
                    zData(j,:) = betainc(1-obsSpec.trialProbability(j),obsSpec.numTrials-binVals, binVals+1);
                case 'Poisson'
                    zData(j,:) = gammainc(meanExp(j), floor(binVals+1), 'upper');
                case 'Gaussian'
                    zData(j,:) = 0.5 * (1 + erf((binVals+0.5-meanExp(j))./sqrt(2*varExp(j))));
            end
        end
    case 'Empirical PMF'
        % Plotting Empirical PMF (i.e., probability mass function) over
        % all sampled observations
        [yData, xData] = histcounts(obsMatrixSampled(:), 'Normalization', 'probability', 'BinMethod','integers');
        numX = length(xData);
        xData = (xData(1:(numX-1))+xData(2:numX))/2; % Convert bin location coordinates to midpoints
    case '3D Empirical PMF'
        % Plotting time-varying PMF. Every sampled time has its own
        % PMF (i.e., probability mass function)
        numYMax = max(obsMatrixSampled(:))+2;
        
        xData = zeros(numPlotInd,numYMax);
        yData = zeros(numPlotInd,numYMax);
        zData = zeros(numPlotInd,numYMax);
        
        for i = 1:numPlotInd
            xData(i,:) = tArrayFull(tPlotInd(i));
            [zDataCur, yDataCur] = histcounts(obsMatrixSampled(:,:,i), 'Normalization', 'probability', 'BinMethod','integers');
            numY = length(yDataCur);
            yDataCur = (yDataCur(1:(numY-1))+yDataCur(2:numY))/2; % Convert bin location coordinates to midpoints
            if yDataCur(1) < 1
                firstPoint = 1;
                lastPoint = numY-1;
            else
                % Force first observation to be 0, even if it isn't
                % observed (will have probability 0 in that case)
                firstPoint = 2;
                lastPoint = numY;
            end
            yData(i,firstPoint:lastPoint) = yDataCur;
            zData(i,firstPoint:lastPoint) = zDataCur;
            % "Extend" display of PMF to global maximum observation (so
            % that time-varying PMF appears square from "above")
            yData(i,(lastPoint+1)) = yDataCur(end)+1;
            yData(i,(lastPoint+2):end) = numYMax;
        end
    case 'Expected PMF'
        % Plotting expected probability mass function based on
        % chosen distribution. If data1 is non-empty, then it is used to
        % define the bin values. Otherwise, we create bin range based on
        % expected mean and variance
        meanExp = obsSpec.trialProbability*obsSpec.numTrials;
        varExp = obsSpec.trialProbability*obsSpec.numTrials*(1-obsSpec.trialProbability);
        if isempty(obsSpec.data1)
            xData = (meanExp-4*sqrt(varExp)):(meanExp+4*sqrt(varExp));
        else
            xData = obsSpec.data1;
        end
        xData = unique(floor(xData(xData>=0))); % Ignore negative and repeated elements
        yData = zeros(1,length(xData));
        switch obsSpec.distribution
            case 'Binomial'
                yData = betainc(1-obsSpec.trialProbability,obsSpec.numTrials-xData, xData+1) - ...
                    betainc(1-obsSpec.trialProbability,obsSpec.numTrials-(xData-1), xData);
            case 'Poisson'
                yData = gammainc(meanExp, floor(xData+1), 'upper') - ...
                    gammainc(meanExp, floor(xData), 'upper');
            case 'Gaussian'
                yData = 0.5 * erf((xData+0.5-meanExp)./sqrt(2*varExp)) - ...
                    0.5 * erf((xData-0.5-meanExp)./sqrt(2*varExp));
        end
    case '3D Expected PMF'
        % Plotting expected probability mass function based on
        % chosen distribution. If data1 is non-empty, then it is used to
        % define the bin values. Otherwise, we center the bins around the
        % maximum expected mean
        % data2 vector is used to define the time-values, so its length
        % should equal that of trialProbability
        numTimeInd = length(obsSpec.trialProbability);
        if length(obsSpec.data2) ~= numTimeInd
            error('Length of observer data2 vector does not equal number of trial probabilities');
        end
        meanExp = obsSpec.trialProbability*obsSpec.numTrials;
        varExp = obsSpec.trialProbability*obsSpec.numTrials.*(1-obsSpec.trialProbability);
        
        if isempty(obsSpec.data1)
            binVals = (min(meanExp)-4*sqrt(min(varExp))):(max(meanExp)+4*sqrt(max(varExp)));
        else
            binVals = obsSpec.data1;
        end
        binVals = unique(floor(binVals(binVals>=0))); % Ignore negative and repeated elements
        % Determine range for each bin
        numBin = length(binVals);
        
        xData = repmat(obsSpec.data2',1,numBin);
        yData = repmat(binVals,numTimeInd,1);
        zData = zeros(numTimeInd, numBin);
        
        for j = 1:numTimeInd
            switch obsSpec.distribution
                case 'Binomial'
                    zData(j,:) = betainc(1-obsSpec.trialProbability(j),obsSpec.numTrials-binVals, binVals+1) - ...
                        betainc(1-obsSpec.trialProbability(j),obsSpec.numTrials-(binVals-1), binVals);
                case 'Poisson'
                    zData(j,:) = gammainc(meanExp(j), floor(binVals+1), 'upper') - ...
                        gammainc(meanExp(j), floor(binVals), 'upper');
                case 'Gaussian'
                    zData(j,:) = 0.5 * erf((binVals+0.5-meanExp(j))./sqrt(2*varExp(j))) - ...
                        0.5 * erf((binVals-0.5-meanExp(j))./sqrt(2*varExp(j)));
            end
        end
    case 'Mutual Information'
        % Plotting sample mutual information based on a reference sample
        xData = tArrayFull(tPlotInd)-tArrayFull(tPlotInd(1));
        yData = zeros(1,numPlotInd);
        
        curMin = min(obsMatrixSampled(:,1,1));
        curMax = max(obsMatrixSampled(:,1,1));
        numCur = curMax-curMin+1;
        for i = 1:numPlotInd
            % 'current' samples are obsMatrixSampled(:,1,1)
            % 'future' samples are obsMatrixSampled(:,1,i)
            futureMin = min(obsMatrixSampled(:,1,i));
            futureMax = max(obsMatrixSampled(:,1,i));
            numFuture = futureMax-futureMin+1;
            
            % Calculate joint PMF between current and future sample times
            jointPMF = zeros(numCur,numFuture);
            for r = 1:numRepeat
                curInd = obsMatrixSampled(r,1,1)-curMin+1;
                futureInd = obsMatrixSampled(r,1,i)-futureMin+1;
                jointPMF(curInd,futureInd) = jointPMF(curInd,futureInd) + 1;
            end
            jointPMF = jointPMF ./ numRepeat;
            
            curPMF = sum(jointPMF,2); % Current marginal PMF
            futurePMF = sum(jointPMF,1); % Future marginal PMF
            
            % Calculate mutual information
            for curVal = curMin:curMax
                for futureVal = futureMin:futureMax
                    curInd = curVal - curMin + 1;
                    futureInd = futureVal - futureMin + 1;
                    if jointPMF(curInd,futureInd) > 0
                        yData(i) = yData(i) + jointPMF(curInd,futureInd) * ...
                            log2(jointPMF(curInd,futureInd)/curPMF(curInd)/futurePMF(futureInd));
                    end
                end
            end
        end
    case '3D Mutual Information'
        % Plotting sample mutual information based on reference samples
        numOffset = length(obsSpec.data1);
        
        xData = zeros(numPlotInd,numOffset); % Reference times
        yData = zeros(numPlotInd,numOffset); % Offset times
        zData = zeros(numPlotInd,numOffset)*NaN; % Mutual information
        
        for t = 1:numPlotInd
            xData(t,:) = tArrayFull(tPlotInd(t));
            yData(t,:) = tArrayFull(tPlotInd(t)+obsSpec.data1)-tArrayFull(tPlotInd(t));
            curMin = min(obsMatrixSampled(:,1,t));
            curMax = max(obsMatrixSampled(:,1,t));
            numCur = curMax-curMin+1;
            for i = 1:numOffset
                % 'current' samples are obsMatrixSampled(:,1,1)
                % 'future' samples are obsMatrixSampled(:,1,i)
                futureMin = min(obsMatrix(iter,1,tPlotInd(t)+obsSpec.data1(i)));
                futureMax = max(obsMatrix(iter,1,tPlotInd(t)+obsSpec.data1(i)));
                numFuture = futureMax-futureMin+1;
            
                % Calculate joint PMF between current and future sample times
                jointPMF = zeros(numCur,numFuture);
                for r = 1:numRepeat
                    curInd = obsMatrixSampled(iter(r),1,t)-curMin+1;
                    futureInd = obsMatrix(iter(r),1,tPlotInd(t)+obsSpec.data1(i))-futureMin+1;
                    jointPMF(curInd,futureInd) = jointPMF(curInd,futureInd) + 1;
                end
                jointPMF = jointPMF ./ numRepeat;
                
                curPMF = sum(jointPMF,2); % Current marginal PMF
                futurePMF = sum(jointPMF,1); % Future marginal PMF
                
                % Calculate mutual information
                zData(t,i) = 0;
                for curVal = curMin:curMax
                    for futureVal = futureMin:futureMax
                        curInd = curVal - curMin + 1;
                        futureInd = futureVal - futureMin + 1;
                        if jointPMF(curInd,futureInd) > 0
                            zData(t,i) = zData(t,i) + jointPMF(curInd,futureInd) * ...
                                log2(jointPMF(curInd,futureInd)/curPMF(curInd)/futurePMF(futureInd));
                        end
                    end
                end
            end
        end
    case 'Monte Carlo Mutual Information'
        % Plotting Monte Carlo mutual information of independent RVs
        numOffset = length(obsSpec.trialProbability);
        numInd = length(obsSpec.data1);
        if ~isempty(obsSpec.data1) && (numOffset == numInd || numOffset <=2)
            xData = obsSpec.data1;
        else
            xData = obsSpec.trialProbability;                
        end
        yData = zeros(1,length(xData));
        numObs = obsSpec.data2(1);
        numRepeat = obsSpec.data2(2);
        
        % Generate random samples
        for r = 1:numRepeat
            curData = zeros(1,numObs);
            for s = 1:numObs
                curData(s) = sum(rand(1,obsSpec.numTrials) < obsSpec.trialProbability(1));
            end
            
            for i = 1:numOffset
                if numOffset <= 2
                    yInd = 1:length(xData);
                    i = numOffset;  % Force a single offset to be checked
                                    % If numOffset == 1, then observations
                                    % will be compared to those with same
                                    % trial probability.
                                    % If numOffset == 2, then comparison
                                    % with observations with same trial
                                    % probability is avoided
                else
                    yInd = i;
                end
                
                if i == 1
                    futureData = curData;
                else
                    futureData = zeros(1,numObs);
                    for s = 1:numObs
                        futureData(s) = sum(rand(1,obsSpec.numTrials) < obsSpec.trialProbability(i));
                    end
                end

                curMin = min(curData);
                curMax = max(curData);
                numCur = curMax-curMin+1;
                futureMin = min(futureData);
                futureMax = max(futureData);
                numFuture = futureMax-futureMin+1;
                
                % Calculate joint PMF between current and future sample
                % "times"
                jointPMF = zeros(numCur,numFuture);
                for s = 1:numObs
                    curInd = curData(s)-curMin+1;
                    futureInd = futureData(s)-futureMin+1;
                    jointPMF(curInd,futureInd) = jointPMF(curInd,futureInd) + 1;
                end
                jointPMF = jointPMF ./ numObs;
                
                % "Marginal" distributions (these are actually independent)
                curPMF = sum(jointPMF,2);
                futurePMF = sum(jointPMF,1);
                
                % Determine empirical mutual information
                for curVal = curMin:curMax
                    for futureVal = futureMin:futureMax
                        curInd = curVal - curMin + 1;
                        futureInd = futureVal - futureMin + 1;
                        if jointPMF(curInd,futureInd) > 0
                            yData(yInd) = yData(yInd) + jointPMF(curInd,futureInd) * ...
                                log2(jointPMF(curInd,futureInd)/curPMF(curInd)/futurePMF(futureInd));
                        end
                    end
                end
            end
        end
        yData = yData./numRepeat;
    case '3D Monte Carlo Mutual Information'
        % Plotting Monte Carlo mutual information of independent RVs
        [numAnchor, numOffset] = size(obsSpec.trialProbability);
        numIndX = length(obsSpec.data1);
        numIndY = length(obsSpec.data2);
        if numAnchor ~= numIndX || numOffset ~= numIndY
            error('Incorrect size of data1 and data2 vectors. Must equal rows and columns of trialProbability, respectively.');
        end
        xData = repmat(obsSpec.data1',1,numOffset);
        yData = repmat(obsSpec.data2,numAnchor,1);
        zData = zeros(numIndX, numIndY);
        numObs = obsSpec.data3(1);
        numRepeat = obsSpec.data3(2);
        
        % Generate random samples
        for r = 1:numRepeat
            for a = 1:numAnchor
                curData = zeros(1,numObs);
                for s = 1:numObs
                    curData(s) = sum(rand(1,obsSpec.numTrials) < obsSpec.trialProbability(a,1));
                end

                for i = 1:numOffset
                    if i == 1
                        futureData = curData;
                    else
                        futureData = zeros(1,numObs);
                        for s = 1:numObs
                            futureData(s) = sum(rand(1,obsSpec.numTrials) < obsSpec.trialProbability(a,i));
                        end
                    end

                    curMin = min(curData);
                    curMax = max(curData);
                    numCur = curMax-curMin+1;
                    futureMin = min(futureData);
                    futureMax = max(futureData);
                    numFuture = futureMax-futureMin+1;
                
                    % Calculate joint PMF between current and future sample
                    % "times"
                    jointPMF = zeros(numCur,numFuture);
                    for s = 1:numObs
                        curInd = curData(s)-curMin+1;
                        futureInd = futureData(s)-futureMin+1;
                        jointPMF(curInd,futureInd) = jointPMF(curInd,futureInd) + 1;
                    end
                    jointPMF = jointPMF ./ numObs;
                    
                    % "Marginal" distributions (these are actually independent)
                    curPMF = sum(jointPMF,2);
                    futurePMF = sum(jointPMF,1);
                
                    % Determine empirical mutual information
                    for curVal = curMin:curMax
                        for futureVal = futureMin:futureMax
                            curInd = curVal - curMin + 1;
                            futureInd = futureVal - futureMin + 1;
                            if jointPMF(curInd,futureInd) > 0
                                zData(a,i) = zData(a,i) + jointPMF(curInd,futureInd) * ...
                                    log2(jointPMF(curInd,futureInd)/curPMF(curInd)/futurePMF(futureInd));
                            end
                        end
                    end
                end
            end
        end
        zData = zData./numRepeat;
    otherwise
        error('Observation type "%s" invalid', obsType)
end

%% Normalize data if specified

% Normalize the x-data?
if obsSpec.bNormalizeX
    % Normalize with respect to what?
    switch obsSpec.normalizeTypeX
        case 'Max'
            maxValX = max(yData(:));
        case 'Custom'
            maxValX = obsSpec.normalizeCustomX;
    end
else
    maxValX = 1;
end

% Normalize the y-data?
if obsSpec.bNormalizeY
    % Normalize with respect to what?
    switch obsSpec.normalizeTypeY
        case 'Max'
            maxValY = max(yData(:));
        case 'Custom'
            maxValY = obsSpec.normalizeCustomY;
    end
else
    maxValY = 1;
end

% Normalize the z-data?
if obsSpec.bNormalizeZ
    % Normalize with respect to what?
    switch obsSpec.normalizeTypeZ
        case 'Max'
            maxValZ = max(zData(:));
        case 'Custom'
            maxValZ = obsSpec.normalizeCustomZ;
    end
else
    maxValZ = 1;
end

%% Plot Curve and Format it
if b3D
    hCurve = surf(hAxes, xData./maxValX, yData./maxValY, zData./maxValZ);
else
    zData = zeros(1,length(xData));
    zData(:) = obsSpec.data3;
    hCurve = plot3(hAxes, xData./maxValX, yData./maxValY, zData./maxValZ);
end
if ~isempty(curveSpec)
    curvePropFields = fieldnames(curveSpec);
    numProp = numel(curvePropFields);
    for i = 1:numProp
        set(hCurve, curvePropFields{i}, curveSpec.(curvePropFields{i}));
    end
end

end