Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
    
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. –  Pickett 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
add comment

2 Answers

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.

share|improve this answer
add comment

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.

share|improve this answer
    
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
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.