@janhardo 


(there is a little piece of C(t)=0 in the bottom left corner)

So I foccused only on the code AFTER your SECOND restart.

In the loop you substitute function f₁(t) by Cforms[i](t)... but u_expt depends on the arbitrary function _F1(t) not on f₁(t).
I corrected this the same way you did in yout procedure PlotPDEsurface.

I understand you want to add legends to each surface in plots_list. am I right?
So I propose you to proceed as described in this sand15.mw file.
Please refer to the comments in worksheet BarGraph_3_Categories.mw  in My Post to see how to adjust manually the rendering of the plot.

@acer 

Thank you on behalf of the author of the message, and I apologize for  having suggested that he should start a new thread.
I suggested this because I wasn't sure whether anyone other than the both of us would be looking at this thread. 

@JoyDivisionMan 

Point 1: "I did not understand the section below, but I think it is a numerical estimate"
Yes: there is no closed form of the CDF so we need to use numerical integration (evalf/Int) instead of int.
The procedure F__U uses two different evalf/Int depending on the interval ([0, 1] ot [1, 2]) the upper bound of the integration belons to.

 

Point 2: "But when I ask for the PDF of Z, I get a bunch of garbage.  But I do get the expected result when I entered plot(f__U(u), u = 0 .. 2)."
PDF(U, u) has support (-2, 2) and 2*PDF(U, u) doesn't change this.
The PDF of Z has support (0, 2) and has the same shape than PDF(U, u) restricted to the interval (0, 2)... but a "vertical scale" has to be done in order that int(PDF(U, u), u=0..2) = 1. Thus the meaning of my writting f__U(u) := 2*PDF(U, u) assuming u > 0.

Maybe this will be more clear to you if you give a look to this more carefully rewitten worksheet where I paid attention to write f__Z instead of f__U everywhere it was needed:
PDF_Z_sand15.mw
(see also the ADD-ON topic below)

Point 3: About your first file:
I use Maple 2014 (see the head of the attached file) and I get no errors: PDFDerivation_MAPLE_2015.mw
This means I can't help you to fix the error you get. Suggestion: open a new thread by asking a question like "Can you help me make yjis code run properly with Maple 2023 ?"


Point 4: About your convolution file:
(Still with Maple 2015), here is a correct way to compute the correlation Convolution_sand15.mw



ADD-ON
Here is a detailed explanation about the the way to get the PDF of Z = |X-Y| for arbitrary random variables X and Y, even if they are not necessary symmetric. This worksheet should also Point 2 clarify above.

 Detailed_explanation_about_convolution.mw

@dharr 

Eriksson et al. write page 3

But from the online site numerical.recipes where the Abramowitz and Stegun's Handbook of Mathematical Functions book can be accessed for free fom its index. But browsing this index I did not find any explicit mention of a modified Bessel functions of the third kind.
So I suspect that some author af one of the earliest paper on NIG distribution used this designation for some reason and that it has persisted over time without anyone really checking its foundations (which is not that unusual in scientific publications).

A little bit of archeology.
_________________________________________________________________________________________
Early appearance of the modified Bessel functions of the third kind:
1979: Halgreen 


__________________________________________________________________________________________

A few years after: a "simple" Bessel function:
1982: Barndorff-Nielsen et al. 
___________________________________________________________________________________________

Fifteen years after :
In 1997 Jaschke p. 5 K_𝝀 is simply "a Bessel function":
___________________________________________________________________________________________

From Bessel to modified Bessel:
In his 1999 PhD about "Generalized Hyperbolic Models" (distributions which encompass NIG distribution for parameter 𝝀=-1/2) Prause writes


Note that papers about NIG distributions use K₁ and instead of K_-1 but assuming that  K₁(z) corresponds to Maple BesselK(𝝀, z) one has  BesselK(-𝝀, z) = BesselK(𝝀, z).
________________________________________________________________________________________

Close encounter of the (modified Bessel of the) third kind
In 2010, Hammerstein's PhD on a quite similar subject writes in Annex 3:

So in a decade second order has become third order.
_________________________________________________________________________________________

The truth is elsewhere (X files)
Finally let us talk Jim Killingsworth, a computer programmer that has been likely as puzzled as I was:
From its online post  Killingsworth  one reads

___________________________________________________________________________________________
So let's say that depending on sense from which the wind blows, or maybe if they want to be felt as a more prominent matheaticians, some will use complicated things to designate simple ones.

@dharr 

Thanks for all these details.

Presently I'm working on those NIG (Normal Inverse Gaussian) distributions and all the papers about them I know use the notation K₁(...). but, despite this notations suggests they refer to the index 1 Bessel function of second kind, all the authors  tell about the incomplete Bessel function of the third kind of index 1.

I know statisticians often use their own notations or definition (and so are their definitions of the erf function and of the Hermite polynomials), but here it is quite confusing.

@dharr 

Thank you for your reply. Far be it from me to overwhelm you with dozens of pages of articles ( ;-) ), but please, just take a quick look at the references in these two.

Garcia Definitions 2.5 and 2.6 (p 26) plus lines which follow in the case the order (aka the index) is an integer ('1' in my case).

In Kucharska & Pielaszkiewicz, § 2.2 The Normal Inverse Gaussian, is given the definition of the Normal Inverse Gaussian distribution which involves the modified Bessel function of the third kind with index 1 (K₁).

The formula corresponds to the probability function and Figure 2.2 on the same page gives an example of I would like to obtain.

In Garcia it is written 

and above

Meanwhile, starting from these two excerpts I did some stuff and finally got results which compare favorably to those of Kucharska & Pielaszkiewicz, Figure 2.2.

Nevertheless they raise questions about the definition of these modified Bessel functions of the third kind which, for integer index/order seem to be nothing but Bessel functions of the second kind with same index/order.

Could it be that the former are extensions to any real index of the latter?

Download Investigation_and_perplexity.mw

What does that mean???

Do you mean that you want to determine its mean, its variance, any other moemnt, its PDF, its CDF? Be more explicit

Note that if "Random variable Y2 is the same as Y1 above" then abs(Y1-Y2) = 0. 
I guess you mean "Y2 has the same distribution than Y1", correct?

"Via simulation, I know that |Y1 - Y2| takes on a logrithmic form".
Seems very strange because the support of Z is obviously (0, 2) and the logarithm function is defined over (0, +oo). Where this knowledge come from? I'm interested because it's wrong. 

Indeed, here is the comparison between a 10⁶-size sample of Z and its analytic PDF(Z):

Note that you can get this result thanks to a simple convolution product.

You say that you "need to get a mathematical solution" (expression of the PDF, CDF, ...?).
Here again be more explicit. Whatever, the expression of the PDF is that complicated that I don' really understand why you need it. More of this the CDF doesn't have a closed form, but cn be assessed numerically:

Feel free to tell me if you are interested to know much about how to get these results.

@Ronan 

Seems simple and effective.

@dharr 

The time scale is indeed important.
What I was referring to was the OP's uncommented answers within, say, a working week (but it may be on vacation).

About your second point "... if there is not an answer in a day or two, the person doesn't check back (and didn't subscribe)."

This is indeed a serious problem, and I must admit that I myself probably neglected to check whether I hadn't received any late replies.  Could it be that some region on the right margin of the "Question" tab signals the arrival of a new contributiuon for, let's say, the ten last questions?
If I'm not mistaken, there is (was?) also the option to be notified directly when a new contribution is published, but this is configured (was configured?) in your profile and therefore applies to all the questions you ask and It's easy to get overwhelmed by the messages that arrive in your mailbox (I used this option years ago but rapidly switched it off for that reason).
One solution could be to access your account via the “Users” tab, where a tab entitled “New contributions” would collect late contributions.

@Andiguys 

I updated my last comment meanwhile.
Sorry for the inconvenience

Reply to Your comment
You write "If I understand correctly, the maximum occurs on the boundary. In that case, what are the analytical first-order/optimality equations for Pn, Pr​, and w?"
The First Order Conditions (FOC) remain the FOC wether the maximum is on the boundary or not. FOC just say that the gradient of the objective function is null where the (unconstrained) objective function reaches an extremum (note the reciprocal is false as the objective function f : x --> x^3 proves, indeed the second derivative must also be non null).

Let's take the simple example of the function f : --> 1-x^2. FOC write 2*x = 0, which implies x*=0 is a mmmiser of f.
But if you search the maximum of f in the interval [1, 2], then the maximiser is x*=1: this does not change the fact that the FOC is still 2*x=0 ... simply the value x=0 cannot be reached given the constraints.
I believe you make a confusion: FOC is a condition for an extrema, but this condition does not have to be fullfilled for the maximum of a constrained problem to exist.


You write "So, does this imply that in both cases there is no maximum solution?"
What does  "there is no maximum solution" mean to you?
If you mean "the objective function does not have a maximum value" you 're wron: any non constant function has one.
If you mean "there is no point M in the interior domain such that, f(m) < f(M) for m in a small neighborhood V(M) of M" you're right": simplythe point M is located on the boundary of the domain defined by the constraints.
So, unless the objective function f is constant, there is always (at least) one point M where f(M) has the largest value.
 

@Andiguys 

This can be easily proved by examining the Hessian of Pi.
It has a null eigenvalue and the condition for tit to be (semi)negative definiteis never realized with your data.
So there is no maximum inside the domain D = (0, 1)x(0, 1)x(0, 1) of variation of (Pn, Pr, w) and Pi reaches its maximum value on some point of the boundary of D.

More of this one can prove that a critical point in D is always a saddle point whatever the data you use.
So, to be clear: you will never get a macimum of Pi within D whatever your data (unless I did something wrong in my demonstration)

Proposal_follow_up_sand15.mw

@JoyDivisionMan

I see you didn't even take time to say if my answer satisfied you.
Is that a simple omission, or rudeness?

@Andiguys 

And whart are the variables relarive to which you want to minimize the objective function?
What about the others: can you provide numeric values or them?

@Andiguys 

Could you simply post a worksheet wihich contains only the objective function and the constraints instead of mixing your new problem with what I did previously?

By the way: it seems quite common for you to ask a question and, once you have received an answer, to ask a completely different question without even saying whether the solution we gave you suits you or not.
I find this particularly unpleasant.