Calculating expectation for a custom distribution in Mathematica

Question

This question builds on the great answers I got on an earlier question:

Can one extend the functionality of PDF, CDF, FindDistributionParameters etc in Mathematica?

To start I have PDFs and CDFs for two custom distributions: nlDist and dplDist as you can see from the code dplDist builds upon nlDist.

    nlDist /: PDF[nlDist[alpha_, beta_, mu_, sigma_], 
   x_] := (1/(2*(alpha + beta)))*alpha* 
   beta*(E^(alpha*(mu + (alpha*sigma^2)/2 - x))* 
      Erfc[(mu + alpha*sigma^2 - x)/(Sqrt[2]*sigma)] + 
     E^(beta*(-mu + (beta*sigma^2)/2 + x))* 
      Erfc[(-mu + beta*sigma^2 + x)/(Sqrt[2]*sigma)]); 

nlDist /: 
  CDF[nlDist[alpha_, beta_, mu_, sigma_], 
   x_] := ((1/(2*(alpha + beta)))*((alpha + beta)*E^(alpha*x)* 
        Erfc[(mu - x)/(Sqrt[2]*sigma)] - 
       beta*E^(alpha*mu + (alpha^2*sigma^2)/2)*
        Erfc[(mu + alpha*sigma^2 - x)/(Sqrt[2]*sigma)] + 
       alpha*E^((-beta)*mu + (beta^2*sigma^2)/2 + alpha*x + beta*x)*
        Erfc[(-mu + beta*sigma^2 + x)/(Sqrt[2]*sigma)]))/ 
   E^(alpha*x);         

dplDist /: PDF[dplDist[alpha_, beta_, mu_, sigma_], x_] := 
  PDF[nlDist[alpha, beta, mu, sigma], Log[x]]/x;
dplDist /: CDF[dplDist[alpha_, beta_, mu_, sigma_], x_] := 
  CDF[nlDist[alpha, beta, mu, sigma], Log[x]];

Plot[PDF[dplDist[3.77, 1.34, -2.65, 0.40], x], {x, 0, .3}, 
 PlotRange -> All]
Plot[CDF[dplDist[3.77, 1.34, -2.65, 0.40], x], {x, 0, .3}, 
 PlotRange -> All]

In my earlier question, joebolte's and sasha's answers and recommendation to use TagSet helped me get this far. Now, my questions relate to the dplDist.

I now need to calculate expectation from some point on the x axis of the PDF. In survival analysis they refer to this as mean residual life. Something like the following:

Expectation[X \[Conditioned] X > 0.1, 
  X \[Distributed] dplDist[3.77, 1.34, -2.65, 0.40]] - 0.1

This doesn't work, essentially just returning the inputs as text.

I understand how I can use TagSet to define PDFs and CDFs for custom distributions, how do I do something similar for Expectation[]?

I'll post more about this follow up in a separate question, but I also need a strategy for calculating goodness of fit of the dplDist relative to some data to which I've fit the distribution.

Many thanks to all.

score 5 · Accepted Answer · answered Jun 09 '11 at 19:58

5

Although you have provided both PDF and CDF for your custom distribution to Mathematica, you have not given the domain, so it does not know boundaries of integration, and in fact whether to integrate or sum. Adding that makes things work:

In[8]:= nlDist /: 
 DistributionDomain[nlDist[alpha_, beta_, mu_, sigma_]] := 
 Interval[{-Infinity, Infinity}]

In[9]:= NExpectation[Log@X \[Conditioned] Log@X > 0.1, 
  X \[Distributed] nlDist[3.77, 1.34, -2.65, 0.40]] - 0.1

Out[9]= 0.199329

Compare this with ProbabilityDistribution that has the format ProbabilityDistribution[ pdf, {x, min, max}], where you explicitly indicate the domain.

In order for symbolic solvers like Probability, Expectation and their numeric counterparts work on these, it is also advised to set DistributionParameterQ and DistributionParameterAssumptions as well.

DistributionParameterQ should give False is parameters explicitly violate assumptions, and DistributionParameterAssumptions should return the boolean expression representing assumptions on parameters of your distribution.

answered Jun 09 '11 at 19:58

Sasha

5,935
1
25
33

@Sasha: This is not related to the question or the answer, but are you at liberty to comment if v9 is expected this year? I'm asking because I've been thinking about getting v8, but if v9 is around the corner, I'd rather wait for it... – abcd Jun 09 '11 at 20:18
@Sasha Many thanks for the answer and the additional insight and heads up on the broader issues. J – Jagra Jun 09 '11 at 20:22
@yoda I am not in position to make such statements, however I would like to point out that purchase of v8 will give you [Premier Service](http://www.wolfram.com/mathematica/how-to-buy/premier-service.html) for 1 year that includes free upgrade to the most recent version released. Historic data on spacings between releases of prior version of [Mathematica](http://en.wikipedia.org/wiki/Mathematica) might help you decide. – Sasha Jun 09 '11 at 21:33
@Sasha: No problem, I understand. Thanks for the additional info re: premier service! – abcd Jun 10 '11 at 00:32
@Sasha: Something doesn't seem to work right. I copied your code and use of DistributionDomain along with my original code into a new notebook and restarted the computer just to have a clean environment, but when I run NExpectation[Log@X \[Conditioned] Log@X > 0.1, X \[Distributed] nlDist[3.77, 1.34, -2.65, 0.40]] - 0.1 I get the following output: -0.1 + NExpectation[Log[X] \[Conditioned] Log[X] > 0.1, X \[Distributed] nlDist[3.77, 1.34, -2.65, 0.4]] instead of what you got: Out[9]= 0.199329 I should have caught this sooner. I'd appreciate any thoughts. Thank. J – Jagra Jun 10 '11 at 16:51
@Jarga Yes, you're right. It turns out both `DistributionParameterQ` as well `DistributionParameterAssumptions` are essential. – Sasha Jun 10 '11 at 21:23

score 1 · Answer 2 · answered Jun 09 '11 at 17:17

1

I'm not sure I really understand your question... The expected value or the mean, is the first moment of the distribution and can be calculated as

expectation := Integrate[x #, {x,-Infinity,Infinity}]&;

and use it as expectation[f[x]], where f[x] is your pdf.

Your last code snippet doesn't work for me. I don't know if it is v8 code or if it is custom defined or if you're trying to say that is what you'd like your function to be like...

You can also try looking into Mathematica's ExpectedValue function.

ExpectedValue[x, NormalDistribution[m, s], x]
Out[1] = m

answered Jun 09 '11 at 17:17

abcd

41,765
7
81
98

The function Expectation[] in version 8 replaces ExpectedValue[] I hoped for a way to use TagSet in something like the following: – Jagra Jun 09 '11 at 17:53
I meant to comment: @yoda The function Expectation[] in version 8 replaces ExpectedValue[] I hoped for a way to use TagSet or something similar in something like the following: dplDist /: Expectation[X, X \[Distributed] dplDist[alpha_, beta_, mu_, sigma_]] := ?; – Jagra Jun 09 '11 at 17:59

score 0 · Answer 3 · edited Apr 13 '17 at 12:55

0

The following page contains some tips on enabling custom distributions (i.e. written from scratch without TransformedDisribution or ProbabilityDistribution) for use in CopulaDistribution, RandomVariate, etc: https://mathematica.stackexchange.com/questions/20067/efficient-generation-of-random-variates-from-a-copula-distribution/26169#26169

edited Apr 13 '17 at 12:55

Community

1
1

answered May 30 '13 at 21:42

doug fromoz

1
2

Please summarise the link in your answer; that way, if the link goes stale the answer won't be completely useless. – michaelb958--GoFundMonica May 30 '13 at 21:58

Calculating expectation for a custom distribution in Mathematica

3 Answers3

Linked