2
$\begingroup$

I want to calculate a distribution from a function that depends on several variables, where there is an interdependency between the variables. This works:

TransformedDistribution[(a - b)* c, 
 {a \[Distributed] NormalDistribution[20, 5], 
  b \[Distributed] NormalDistribution[30, 2], 
  c \[Distributed] NormalDistribution[100, 8]}]

RandomVariate[%]
(* -442.531 *)

Since there are no interdependencies. This, however, does not:

TransformedDistribution[(a - b)* c,
 {a \[Distributed] NormalDistribution[20, 5], 
  b \[Distributed] NormalDistribution[a, 2], 
  c \[Distributed] NormalDistribution[100, 8]}]

RandomVariate[%]

During evaluation of In[175]:= NormalDistribution::realprm: Parameter [FormalX]1 at position 1 in NormalDistribution[[FormalX]1,2] is expected to be real.

I found that ParameterMixtureDistribution works for this:

ParameterMixtureDistribution[
 TransformedDistribution[(a - b)*c,
  {b \[Distributed] NormalDistribution[0.75*a, (0.5 a - 1)/6], 
   c \[Distributed] NormalDistribution[100, 8]}], 
   a \[Distributed] NormalDistribution[20, 5]]

RandomVariate[%] 
(* 741.52 *)

But, not for this:

ParameterMixtureDistribution[
 TransformedDistribution[(a - b)*Log[c],
  {b \[Distributed] NormalDistribution[0.75*a, (0.5 a - 1)/6], 
   c \[Distributed] NormalDistribution[100, 8]}], 
   a \[Distributed] NormalDistribution[20, 5]]

RandomVariate[%] 

I don't get an error...but I'm not getting useful output: it just returns the input.

$\endgroup$
1
  • $\begingroup$ I think my question may not have been well formulated. What's the etiquette? Should I delete the question altogether and ask a new questions? Severely edit it? $\endgroup$
    – Brian
    Commented Jun 30, 2023 at 12:12

2 Answers 2

4
$\begingroup$

This is an extended comment echoing the good advice from @heropup. While ParameterMixtureDistribution function is very handy, your second example doesn't need it. Writing things out in terms of the random variables we have

$$A\sim N(20,5)$$ $$C\sim N(100,8)$$ $$Z\sim N(0,2)$$ $$B=A+Z$$

So

$$(A-B)C=(A-(A+Z))C=Z C$$

And you can get the desired result with

TransformedDistribution[z * c,
 {z \[Distributed] NormalDistribution[0, 2], 
  c \[Distributed] NormalDistribution[100, 8]}]

For your second example suppose you fix the issue that @heropup pointed out with making the standard deviation always non-negative with Abs[a/2 - 1]/6 and changing $\log C$ to $\log{|C|}$ (as I don't thing you want to deal with complex numbers - at least I assume you don't):

Then

$$B = 3A/4 + Z |A/2 - 1|/6$$

The Mathematica code becomes

dist = TransformedDistribution[(aa - 3 aa/4 - zz Abs[aa/2 - 1]/6) Log[Abs[cc]],
  {aa \[Distributed] NormalDistribution[20, 5], 
   cc \[Distributed] NormalDistribution[100, 8], 
   zz \[Distributed] NormalDistribution[0, 1]}]
SeedRandom[12345];
RandomVariate[dist]
(* 25.4829 *)
$\endgroup$
3
  • $\begingroup$ Thanks for the answer. As I mentioned in my comment to @heropup, I think my question has been too watered down. In your answer, you took advantage of the fact that a Normal distribution simply has its mean shifted when adding a value to it. Clever! My intent, however, was to understand how to generate a distribution that depends on the output of another distribution. What if, for example, the mean and standard deviation of b depended on a? Or if b were a Transformed distribution, with a transformation that used the output of a? $\endgroup$
    – Brian
    Commented Jun 30, 2023 at 12:10
  • $\begingroup$ "It depends." It depends on the specifics of the example. Some scenarios have nice results for the moments but not for the resulting pdf. Also, there are "You can't get there from here" issues. An example is that ParameterMixtureDistribution gives you the marginal distribution for $B$. But you need the joint distribution of $B$ and $A$ in your examples. In short, you'll need to give relatively specific examples. If you edit your question, I'm happy to delete my "extended comment." I don't know the recommended etiquette. $\endgroup$
    – JimB
    Commented Jun 30, 2023 at 13:19
  • $\begingroup$ "My intent, however, was to understand how to generate a distribution that depends on the output of another distribution." Understood. There are statistical theory approaches and what Mathematica has currently implemented in nice compact functions and other more basic Mathematica functions that require more thought. $\endgroup$
    – JimB
    Commented Jun 30, 2023 at 13:25
4
$\begingroup$

Well, you have some fundamental problems with the model specification even when you are able to obtain a result: your standard deviation is negative if $a < 2$, which occurs with approximate probability $0.000159109$. That may seem negligible but it is once in approximately every $6285$ simulations.

In order to address this you would need to adjust your definitions; e.g., left-truncate $a$. A similar but less egregious problem happens if you attempt to compute $\log c$: this will not be real-valued with probability $3.73256 \times 10^{-36}$ which is computationally negligible. In any case, even if $c < 0$, Mathematica would simply compute a complex-valued logarithm for the realization. But it will refuse to generate a random variate from a normal distribution with negative or zero standard deviation.

As a result, for your specific inputs, I would not trust the output of RandomVariate even when it does give you one. If your intent is to generate random variates from such a model, you are better off doing it manually; i.e., generate a variate for $a$, then use this to generate $b$ (where truncation could be implemented by discarding any realizations where $a \le 2$); then generate $c$, and finally compute whatever function of $a, b, c$ is desired. It is also likely to be faster.


Running the following code on a 3-year old laptop takes about 7 seconds and generates nearly $10^7$ variates, a full order of magnitude more than the number you attempted to generate.

Select[ParallelTable[If[# > 2, (# - RandomVariate[NormalDistribution[.75 #, (.5 # - 1)/6]]) RandomVariate[NormalDistribution[100, 8]]] &[RandomVariate[NormalDistribution[20, 5]]], 10^7], NumberQ]

Meanwhile, my attempt to use your code with ParameterMixtureDistribution is so slow, it takes 16 seconds to just get $10^4$ variates. I don't claim that my code is optimized; in fact, I think there are many other people here who could do it better and faster than this. Since I don't know what code you used for your manual attempt I can't say why your implementation was slow.

If you are needing to combine this with other distributions, you should have included this in the scope of your original question, rather than asking about one part that might not fully capture the issue you are facing. But as I have already stated, you have more fundamental problems with your model specification.

$\endgroup$
3
  • $\begingroup$ Thanks for the response. I hope you'll forgive me that I tried to reduce my problem to the essential problem I was having to keep things simpler, but I think I didn't do that well. I have been using TruncatedDistribution[] and other means to prevent issues as you described, but left them out to simplify the code here. $\endgroup$
    – Brian
    Commented Jun 30, 2023 at 11:57
  • $\begingroup$ I also tried doing it "manually," but found this was actually very, very slow. I created a Module[], and generated RandomVariates for a, b, and c, in sequence, then calculated (a-b)*Log(c). Finally, I wrapped that in Table[] to generate a list of results. 1,000,000 iterations takes > 260 seconds. Since I want to convert this to a distribution (with EmpiricalDistribution[]?) to combine with still more distributions, the time is prohibitive. Is there a faster way? $\endgroup$
    – Brian
    Commented Jun 30, 2023 at 12:00
  • $\begingroup$ @Brian See my edits above. $\endgroup$
    – heropup
    Commented Jun 30, 2023 at 15:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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