4

I am still pretty new to Mathematica, but have to solve a quite complex problem.

The function

kdis[x_, a_] := KernelMixtureDistribution[Inner[List, x, (1+x)^a, List]]

gives me a bivariate kernel density estimate of two gamma random variables given by

gamran[k_, t_, nr_] := RandomVariate[GammaDistribution[k, t], nr]

Given some observations

xy = {{x_1,y_1},{x_2,y_2},...,{x_n,y_n}}

, I seek to maximize the likelihood function

Likelihood[kdis[gamran[2, 0.5, 250000], gamran[k, t, 250000]], xy]

using either NMaximize or FindDistributionParameters.

NMaximize[{Likelihood[kdis[gamran[2, 0.5, 250000], gamran[k, t, 250000]], xy], 0 < k <= 10 && 0 < t <= 10}, {{k, 0.01, 10}, {t, 0.01,10}},StepMonitor :> Print[k]]

FindDistributionParameters[xy,kdis[gamran[2.7, 0.5, 250000], gamran[k,t,250000]],{{k,1},{t,1}}, ParameterEstimator -> "MaximumLikelihood"]

Both functions seem to run endlessly. In particular, NMaximize never even finishes the first step (StepMonitor never prints something). As the problem is constrained to positive values of $k$ and $t$, I define constraints and a positive initial region for NMaximize. However, I get the following error message:

GammaDistribution::posprm: "Parameter k at position 1 in GammaDistribution[k,t] is expected to be positive."

which makes no sense to me given the constraints and the intial region.

I know the problem is computationally very demanding. However, plotting the Likelihood with ContourPlot takes some time, but works and gives me relatively well defined maxima within my initial region. So I dont quite understand what is wrong with my script.

PERG
  • 51
  • 3
  • You define kdis but use bkdis. Also, for the second argument of Likelihood you probably meant {x,y}. It could help to use _?NumericQ as in this. – C. E. Nov 11 '13 at 10:07
  • Thank you... I edited the original post accordingly. However, inserting _?NumericQ seems to have no effect. I am still facing the error message... – PERG Nov 11 '13 at 10:20

2 Answers2

1

I finally solved the problem and the script is running now (though it still needs very long to get some results).

The problem was somehow related to the gamran function I defined. However, rewriting my kdis function and explicitely assigning _?NumericQ solved the problem:

kdis[x_, k_?NumericQ, t_?NumericQ, nr_?NumericQ] := KernelMixtureDistribution[Inner[List, x, (1 + x)^RandomVariate[GammaDistribution[k, t], nr], List]]

Most probably the composition of 3 different functions (Likelihood,kdis,gamran) was kind of awkward.

The error message is gone, as well... but you have to give the correct constraints k>0 and t>0.

PERG
  • 51
  • 3
0

My apologies if I have misunderstood your intention or your requirements: Consider test as product distribution of two gamma distributions,one with known parameters and the other unknown. The following test example may illustrate:

test = TransformedDistribution[{x, 
   y}, {x \[Distributed] GammaDistribution[2, 0.5], 
   y \[Distributed] GammaDistribution[0.2, 3]}]
data = RandomVariate[test, 10000];
EstimatedDistribution[data, 
 ProductDistribution[GammaDistribution[2, 0.5], 
  GammaDistribution[a, b]]]
FindDistributionParameters[data, 
 ProductDistribution[GammaDistribution[2, 0.5], 
  GammaDistribution[a, b]], ParameterEstimator -> "MaximumLikelihood"]

yields:

ProductDistribution[GammaDistribution[2, 0.5], 
 GammaDistribution[0.202338, 2.85498]]

{a -> 0.202338, b -> 2.85498}

Again apologies if I have misunderstood.

ubpdqn
  • 60,617
  • 3
  • 59
  • 148
  • Thank you for your suggestion... but unfortunately my problem is not related to a simple product distribution. The dependency structure is actually quite complicated. However, do you have a clue what leads to the error message? – PERG Nov 11 '13 at 13:17