[transformedHistogram,newX]=transformHistogram(originalHistogram,originalX,formula)

function [transformedHistogram,newX]=transformHistogram(originalHistogram,originalX,formula)
% [transformedHistogram,newX]=transformHistogram(originalHistogram,originalX,formula)
% written by Peter Nagy, [e-mail]
%
% transformHistogram transforms the histogram stored in the 1D array
% 'orignalHistogram' according to the formula 'formula'.
% 'originalHistogram' is defined at bins stored in 'orignalX'.
% originalHistogram is assumed to be stored according to the following principle:
% f(X(i)) = frequency (X(i) < Y <= X(i+1))
% f(X(end) = frequency (Y = X(end))
% This is the principle how histc generates its output. 
% The best is if f(end)=0
% The output 'transformedHistogram' contains the transformed histogram at
% bins 'newX'. newX defines the same intervals as originalX, but the range
% of newX may be smaller, if the formula moves the histogram in the
% negative direction on the X scale. 
% The last bin of transformedHistogram is empty, i.e.
% length(transformedHistogram)=length(newX)-1 (the last frequency
% corresponding to Y = X(end) is missing).
% 'formula' has to be a string containing the transformation rules. The variable has to be x, e.g. 2x+3
% Example: transformHistogram(starved,xscale,'x.^2')
% Pay attention to the dot (.^), element-wise operation.

% check input
if ~isequal(size(originalHistogram),size(originalX))
    error('Sizes of originalHistogram and originalX have to be equal');
end

% generate new X axis
x=originalX;
transformingFormula=@(x) eval(formula);
transformedX=transformingFormula(x);

% calculate derivative of formula
syms x;
symbolicFormula=eval(formula);
derivative=diff(symbolicFormula);
derivative=@(x) eval(derivative);

% generate probability density function of original histogram
originalPDF=originalHistogram(1:end-1)./diff(originalX);

% calculate probability density function of transformed histogram
transformedPDF=originalPDF./derivative(originalX(1:end-1));

% if the range of the transformedX axis is smaller than that of the
% original, reduce it
if max(transformedX)<max(originalX)
    newX=originalX(1:find(originalX<=max(transformedX),1,'last'));
else
    newX=originalX;
end

% calculate the PDF of the transformed distribution using linear
% interpolation at values corresponding to the original X scale
transformedPDF2=interp1(transformedX(1:end-1),transformedPDF,newX);

% calculate the transformed histogram
transformedHistogram=transformedPDF2(1:end-1) .* diff(newX);