Model parameter distributions from Gamma distributed data

[mikeclimber1]

I have a set of data points {x_i,y_i} where each y_i is a Gamma distributed variable where both the shape k and scale \theta depend on i.

I then fit the data points with a power law model y=a(x)^b.

I would like to know the probability distributions for the fit parameters a and b.

Is there an analytical approach for this problem? The only method I can think of is to simulate a bunch of data sets and manually build the distributions of a and b, at which point I could fit the distributions.

 

[pmsrw3]

What objective function did you use? How many points, and are they independent? There's a general theorem that says that maximum likelihood estimates based on independent variates are asymptotically normal.

Other than that, I would be surprised if there's anything substantially better than the Monte Carlo approach you mention.

 

[mikeclimber1]

By taking the log of each side ln(y)=b ln(x)+ln(a), than an ordinary least square fit can be used I think.

There are four data points and they are independent. To give you an idea of what I'm dealing with, the last point {x₄,y₄}, is distributed with a Gamma function with a shape parameter k=1, meaning that it's exponential, whereas for {x₁,y₁}, the Gamma distribution is closer to normal. I also know that shape and scale parameters vary with i in such a way that the mean scales linearly with i. This means that on average, the four data points are linear and b\approx1. However, I believe that the distribution of b is non-normal, given the strong asymmetry in how some of the data points are distributed. Hopefully this makes sense.

I've used the Monte Carlo approach and found that distribution of b has a mean value of \approx1, but it is very asymmetric. I was hoping that there was some analytical technique to determine the functional form of the distribution of b.

 

[pmsrw3]

Oh, only four points. Yeah, the asymptotic properties aren't going to be much help.

Well, I'm going to tentatively suggest, instead of a least squares fit, you look at a http://en.wikipedia.org/wiki/Generalized_linear_model" . Common stats packages will do this, there's some analytical backing (but don't ask me about it -- I'm ignorant, you'll have to read up), and it will allow you to deal explicitly with the fact that you have gamma-distributed data.

Mucho gusto.

 

[mikeclimber1]

Thanks for the suggestion. I looked a little at generalized linear models because it allows for the dependent variables to be generated from any distribution of the exponential family. But I'm also ignorant, so I haven't figured out how to use it yet.

Thanks again.

 

[pmsrw3]

I actually have done GLM fits. Mathematica, R, and SPSS all have GLM fit functions -- I just used them (Mathematica to be precise). For my application, it was simple and worked very well -- much better than a least squares fit.

 

[bpet]

If using linear regression the explicit formula would be

a = exp(B-AC/D)
b = C/D

where A=E[log(X)], B=E[log(Y)], C=E[log(X)log(Y)]-AB, D=E[log(X)^2]-A^2 where the E[...] is the sample mean, so I'd be surprised if there is an explicit formula for the distribution.

Edit: if the Xj are fixed then effectively log(a) and b are linear combinations of the log(Yj), so perhaps you could try inverting the characteristic function, which can be written in terms of

E[\exp(itY_j)]=\theta_j^{it}\Gamma(it+k_j)/\Gamma(k_j)