How can I model a distribution for a data set if I know how each sample behaves

I will attempt to describe what the data is although it isn't central to the question I'm asking.

My data set is the spectral attenuation of a particular type of brain tissue as recorded on a dual energy CT machine. Essentially each same records the Hounsfield unit from a particular ROI at 7 different x-ray energy levels.

My data can be split into 2 subsets. For one set the sample can be modeled by the function y=a*x^b, and the other can be modeled by y=a*x^-b.

Equivalently I can say that the logarithm of the data is linear, according to log(y)=a*log(x)±b.

Each particualar sample in the set is characterized by its coefficients (a,b). I would like to find a linear regression model that will describe the entire data set, thus allowing me to judge for a particular pair of (a,b) where it falls relative to the distribution, or general line.

I hope the question is clear. Any help would be appreciated.

Answer 1

[John D'Errico]

So you have two datasets, thus sets of (x,y) pairs. Almost always it will be appropriate to log the model. Data that follows models like this often has proportional error structures anyway, thus noise that follows a lognormal distribution, or something close. This is especially important if your data spans more than a power of 10 or so in y. Otherwise, some data points will have too large of an impact on the fit.

So definitely log your data. If the models are:

y1 = a*x1^b
y2 = a*x2^(-b)

Then the logged models will be

log(y1) = log(a) + b*log(x1)
log(y2) = log(a) - b*log(x2)

Note that you wrote the log of the model incorrectly.

I'll assume that x1, y1, x2, y2 are column vectors containing the two data sets.

n1 = length(x1);
n2 = length(x2);

coefs = [ones(n1+n2,1),[log(x1);-log(x2)]]\[log(y1);log(y2];
a = exp(coefs(1));
b = coefs(2);