# Maximum Likelihood Estimation

Posted 10 years ago
21475 Views
|
9 Replies
|
6 Total Likes
|
 How come there is no discussion about MLE in the community?Is it so straightforward to implement it in Mathematica?I have been trying for a while to replicate the results from Barlow (2002)  - "A diffusion model for electricity prices", but I didn't succeed.Does anyone have any similar code, or any suggestions?Letting DD be the vector of prices, here is what I have done:And then, assuming some initial values, I use:FindMaximum[{obslogL && \[Sigma] > 0 && \[Lambda] != 0 && \[Alpha] != 0}, {{\[Lambda], 10}, {a, 4}, {\[Sigma], 2}, {\[Alpha], 1.3}}]However, I get the following error:"FindMaximum::nrnum: The function value -False is not a real number at {\,a,\,\} = {9.0541,1.43219,-0.40076,0.265587}. >>"and the hopeless result: {-737.874, {\ -> 9.99596, a -> 4.00037, \ -> 2.67958, \ -> -4.17363}}
9 Replies
Sort By:
Posted 9 years ago
 I have a somewhat related question, which could also provide an answer to another unanswered question of mine: TimeSeriesModelFit // (G)ARCH estimating methods Hopefully this thread will be revised by someone after so long time. How do I maximize the MLE which contains many variables in a reasonable time? Below I have uploaded a notebook where I try to find 103 parameters.Excel's Solver can do the maximization with less parameters relatively easy (I didn't try with so many), but of course there is a precision concern. MaxBFGS from OxMetrics, on the other hand, does it reasonably fast.FindMaximum in Mathematica, however, takes years... Attachments:
Posted 10 years ago
 There may not be "one theta to rule them all"  Compare LogLikelihood values of the results: In[1]:= fun[   x_] := (1 - c)*(1/(E^((-a + x - c*x)^2/(2*b^2))*(b*Sqrt[2*Pi])))  In[2]:= dist =    ProbabilityDistribution[fun[x], {x, -Infinity, Infinity}];  In[3]:= data =    RandomVariate[dist /. {a -> 0, b -> 1, c -> 1/3}, 10^3,     WorkingPrecision -> 20];In[4]:= res = Table[Quiet@   FindDistributionParameters[data,     dist, {{a, .2}, {b, .8}, {c, start}},     ParameterEstimator -> "MaximumLikelihood",     WorkingPrecision -> 20], {start, .1, .9, .1}]Out[4]= {{a -> -0.019932674343740414493, b -> 1.1096360920952880069,   c -> 0.26954206522212784697}, {a -> -0.018766451751580295519,   b -> 1.0447134039438880896,   c -> 0.31227976365858204554}, {a -> -0.016674487536403968415,   b -> 0.92825542430454599255,   c -> 0.38894242441693089468}, {a -> -0.015633567435080450467,   b -> 0.87030823228157333632,   c -> 0.42708825124255829112}, {a -> -0.014913446774718998486,   b -> 0.83021968892034345328,   c -> 0.45347797919258945110}, {a -> -0.014144317149532428768,   b -> 0.78740285539643202694,   c -> 0.48166370243199674812}, {a -> -0.013811227971511007229,   b -> 0.76886004700210827827,   c -> 0.49387017410287984506}, {a -> -0.015406220411774267890,   b -> 0.85765200449162658321,   c -> 0.43541966915346278263}, {a -> -0.015791913582145423886,   b -> 0.87912323840309263789, c -> 0.42128545588647653513}}In[5]:= LogLikelihood[dist /. #, data] & /@ resOut[5]= {-1837.0542834440804521, -1837.0542834440804521, \-1837.0542834440804521, -1837.0542834440804521, \-1837.0542834440804521, -1837.0542834440804521, \-1837.0542834440804521, -1837.0542834440804521, -1837.0542834440804521}(*compare to the "original" distribution*)In[6]:= LogLikelihood[dist /. {a -> 0, b -> 1, c -> 1/3}, data]Out[6]= -1837.3811539068860566(*the original distribution may not be the best fit to the data, even \though it was used to generate said data*)
Posted 9 years ago
 Hi mam what about simulation of likelihood eatimates and their variances?
Posted 10 years ago
 Thank you very much, Peter!It does help.. somehow, but I still have problems with it.Namely, I get different results if I start with a different "theta" for example.
Posted 10 years ago
 Hi,a simple made up example to illustrate use of ProbabilityDistribution: In[49]:= fun[   x_] := (1 - c)*(1/(E^((-a + x - c*x)^2/(2*b^2))*(b*Sqrt[2*Pi])))  In[50]:= Integrate[fun[x], {x, -Infinity, Infinity},   Assumptions -> b > 0 && 0 < c < 1] Out[50]= 1  In[51]:= dist =    ProbabilityDistribution[fun[x], {x, -Infinity, Infinity}];In[52]:= data =   RandomVariate[dist /. {a -> 0, b -> 1, c -> 1/3}, 10^3];In[53]:= FindDistributionParameters[data, dist, {{a, .2}, {b, .8}, {c, .3}},  ParameterEstimator -> "MaximumLikelihood"]Out[53]= {a -> 0.0312187, b -> 0.900318, c -> 0.368235}As for FindMaximum I would try to avoid != and use inequalities instead, maybe even split the region.
Posted 10 years ago
 Sandu,ProbabilityDistribution allows you to define your own Distribution from any given pdf or cdf, symbolic or numeric, while EmpiricalDistribution (and cousins) creates distributions from given data.http://reference.wolfram.com/mathematica/ref/ProbabilityDistribution.htmlhttp://reference.wolfram.com/mathematica/ref/EmpiricalDistribution.htmlHope this helps.
Posted 10 years ago
 Dear Sandu, we recommend avoiding using Subscripts in community post - they increase code volume and make it unreadable. We suggest to use, for example, instead of thisSubscript[g, ?]simplyga- plain English names for variables.
Posted 10 years ago
 EstimateDistribution requires a specified distribution.Mine, has the following PDF:p(Si,Si-1) = |g'[Si| * 1/Sqrt[2 Pi Theta] * E^(-(g[Si] - rho*g[Si-1] - b)/(2 Theta)))^2Where, theta, rho and b are known fucntions of my parameters.It looks like a normal one, but it isn't.If could define a distribution from this PDF, then I guess I'm done.
Posted 10 years ago
 Sandu,I'm not sure I completely understand what you are trying to do, but would EstimatedDistribution with ParameterEstimator -> "MaximumLikelihood" be of any help?see http://reference.wolfram.com/mathematica/ref/EstimatedDistribution.html section Details and OptionsPeter