## 5259 Reputation

7 years, 110 days

## @NanaP Thanks.A quick reply to your...

@NanaP

ALL_Solutions file updated (see below)

Thanks.

---------------------------------------------------------------------
1
Let X some random variable anf f its PDF.
Let Y = a+b*X (b assumed <> 0) and g its PDF; the basic formula to keep in mind is

`g(t) = f((t-a)/b) / |b|`

---------------------------------------------------------------------
2
Concerning the CDF the things are a little bit more complex. See an example here answ_1.mw.
If the scaling factor b is assumed to be strictly positive and if F and G are the CDF of X and Y respectively, then the second basic formula is

`F(t) = F((t-a)/b)`

Demonstration of these formulas can be found in any text cours of Probability.

---------------------------------------------------------------------
3
Whatever the value of the scaling parameter Variance(Y) = b^2 * Variance(X).
As Variance is a central moment the shift parameter a doesn't matter.

---------------------------------------------------------------------
4
Skewness and Kurtosis are both defined based on central moments and the shift parameter a doesn't matter.
If you closely look to their formulas (see the Maple help pages for instance), you will see that they are dimensionless
(
Let MX(r) the central moment of order r of X and  MY(r) the central moment of order r of Y.
One has  MX(r) = br * MY(r).
Skewness(X) = MX(3) / MX(2)3/2 then Skewness(Y) = b3 *MX(3) / ( b2 *MX(2)) 3/2 = MX(3) / MX(2)3/2  = Skewness(X).
For the same reason he same Kurtosis(X) = Kurtosis(Y)
).

---------------------------------------------------------------------
​​​​​​​5
The Mean of a random variable (more precisely its Expectation), is linear operator over random variables.
Let E this operator.
Then E(Y) = E(a + b*X) = a + b*E(X).
Concerning the Mode: the transformation T : X --> Y being linear one easily deduces that Mode(Y) = a + b*Mode(X).
Concerning the Median, which is by defininition the value such that the PDF of a random variable takes the value 1/2, the demonstration that Median(Y) = a + b*Median(X) is a little bit more complex, see here answ_2.mw

FINALLY:
The following file contains all previous step-by-step demonstrations with Maple ​​​​​​​

 > restart

Probability Density Functions

X is a random variable with PDF f__X
Y = a + b . X
f__Y = PDF(Y)

 > with(IntegrationTools):
 > # X to Y relation # # For the sake of simplicity I consider only the case b > 0. Y = a+b*X
 (1)
 > # Defition Prob(X <= u) = Int(f__X(t), t=-infinity..u)
 (2)
 > # and Prob(Y <= v) = Int(f__Y(t), t=-infinity..v)
 (3)
 > # Use relation (1) eval((3), (1))
 (4)
 > # Rearrange the operand of Prob Prob(isolate(op(lhs((4))), X)) = rhs((4)) assuming b > 0
 (5)
 > # Which means, from definition (2) that: eval((2), u=rhs(op(lhs((5)))));
 (6)
 > # Thus equating rhs(6) and rhs(5) gives rhs((6)) = rhs((5))
 (7)
 > # Comparing the upper bounds of the two Integrals in (7) suggests # that the change of variable t --> (tau-a)/b in the rhs Integral # will change its upper bound (v-a) b Into v: Change(lhs((7)), t=(tau-a)/b, tau) = rhs((7)) assuming b > 0
 (8)
 > # Gather the two members within the same integral Combine((lhs-rhs)((8))) = 0
 (9)
 > # As this equality holds for any real value v one gets GetIntegrand((9)) = 0
 (10)
 > # And finally simplify(isolate((10), f__Y(tau)), size)
 (11)
 > # QED

Cumulative Density Functions

CDF(X) =F__X(t)
CDF(Y) = F__Y(t)

 > # Integrate both sides of (11) from -oo to some arbitrary value T map(z -> Int(z, tau=-infinity..T), (11))
 (12)
 > # By definition of F__Y: F__Y(T) = Change(rhs((12)), tau=a+b*t, t) assuming b > 0
 (13)
 > # By definition of F__X: lhs((13)) = F__X(op(2, GetRange(rhs((13)))))
 (14)
 > map(simplify, %);
 (15)
 > # QED

Means

mu__X = mean(X)
mu__Y = mean(Y)

 > # By definition: mu__X = Int(t*f__X(t), t=-infinity..+infinity);
 (16)
 > # and: mu__Y = Int(tau*f__Y(tau), tau=-infinity..+infinity);
 (17)
 > # Account for relation (11) eval((17), (11))
 (18)
 > # Expand this expression after having use this change of the integration # variable: tau = a + b*t Expand(Change((18), tau=a+b*t, t)) assuming b > 0
 (19)
 > # As f__Y is a PDF, one has Int(f__X(t), t = -infinity .. infinity) = 1
 (20)
 > # Simplify relation (19) using (20): eval((19), (20))
 (21)
 > # Use the definition (16) to simplify (21) eval((21), (rhs=lhs)((16)))
 (22)
 > # QED

Variances

V__X = variance(X)
V__Y = variance(Y)

 > # By definition: V__X = Int((t-mu__X)^2*f__X(t), t=-infinity..+infinity);
 (23)
 > # and: V__Y = Int((tau-mu__Y)^2*f__Y(tau), tau=-infinity..+infinity);
 (24)
 > # Account for relation (11) eval((24), (11))
 (25)
 > # Use the change tau = a + b*t for the integration variable Change((25), tau=a+b*t, t) assuming b > 0
 (26)
 > # Account for relation (22) eval((26), (22))
 (27)
 > # Factorize the integrand lhs((27)) = Int(factor(GetIntegrand(rhs((27)))), t=GetRange(rhs((27))))
 (28)
 > # Transform relation (23) this way Expand((23)): isolate(%, select(has, [op(rhs(%))], t^2)[])
 (29)
 > # Expand relation (28) Expand((28))
 (30)
 > # Modify relation (30) while accounting for relation (29) simplify(eval((30), (29)))
 (31)
 > # QED

Medians

Medians

 > # By definition of the median: Def__X := Int(f__X(t), t=-infinity..Median__X) = 1/2;
 (32)
 > # and Def__Y := Int(f__Y(tau), tau=-infinity..Median__Y) = 1/2;
 (33)
 > # Account for relation (11) in (33) eval((33), (11))
 (34)
 > # Change the integration variable Change((34), tau=a+b*t, t) assuming b > 0
 (35)
 > # Use definition (32) eval((35), (rhs=lhs)((32)))
 (36)
 > # Equating the upper bounds of each member gives op(2, GetRange(lhs((36)))) = op(2, GetRange(rhs((36))))
 (37)
 > # Isolate Median__Y isolate((37), Median__Y)
 (38)

Central Moments anq Quantiles

 > # For any central moments or order r (r >= 2) the proof that CentralMonent(Y, r) = b^r * CentralMoment(X, r); # proceeds exactly as the one used for the variances # (central moment of order 2).
 (39)
 > # For any Quantile the proof that Quantile(Y, prob) = a + b * Quantile(X, prob); # proceeds exactly as the one used for the medians # (quantile associated to prob=1/2).
 (40)
 >

Here is the demonstration that your problem may have no solution in integer numbers.

I consider the elementary case of 2x2 random matrices of integer generated this way:

`RandomMatrix(2\$2, generator=-10..10);`

All these matrices have a determinant in the range -200..200.
Nevertheless one can proove these determinants can only take 333 values out of the 401 possible a priori.
This means, for instance, that there exist no such matrix whose determinant is 173.

In more general terms your problem doesn't always has a solution: this will depend on the value of the determinant you give.
(note that it has always an infinity of solutions in real numbers, see my answer).

 > restart:
 > randomize(168749831218663);
 (1)
 > with(LinearAlgebra): with(Statistics):

Preliminary investigation.

 > # It is obvious that the determinant of matrices like M below is in [-200, 200]. N := 2: M := RandomMatrix(N\$2, generator=-10..10);
 (2)
 > # Modify matrix M by randomly replacing one of it's entry by a name x. # For instance: r := rand(1..N^2): [indices(M)][r()]; M[op(%)] := x: d := Determinant(M): M, d;
 (3)
 > # Let T a target determinant and let's assume T is an integer. # # Then x must be such that SearchedValue := x = solve(d=T, x);
 (4)
 > # And if we still want x to be an integer, T/10+1/2 must be one too. # Given T is in [-100, 100], only values of T such that FeasibleT := {T=10*k+5, k >= -9, k <= 9}
 (5)
 > # Thus T must one of Dom := {seq(-95..95, 10)}
 (6)
 > # Let's take T=15: seq(eval(SearchedValue, T=guess), guess in Dom);
 (7)

A brute force method (before a smarter analysis).

 > r := rand(-10..10): NR := 10^6: dets := Vector(NR): for nr from 1 to NR do   dets[nr] := r()*r()-r()*r() end do:
 > Ta := Tally(dets, output=table);
 (8)
 > {indices(Ta, nolist)}; numelems(%);
 (9)
 > ho := Histogram(Select(z->is(z::odd), dets), discrete=true, color=red): he := Histogram(Select(z->is(z::even), dets), discrete=true, color=blue): plots:-display(ho, he, size=[1000, 400])
 > # Observe the pre-eminence of even determinants numelems( Select(t -> is(t::even), dets) ) / NR; evalf(%);
 (10)

A smarter analysis.

 > # Consider all possible outcomes 's' of generator=-10..10 s := <[\$-10..10]>: # 's' contains 11 out of 21 even values. # Consider the product 's2 = s.s^+' (which represents for instance the product # of diagonal elements in the expression of the determinant). # The probability that this product is even is #       even * even       even * odd        odd  * even #     (11/21)*(11*21) + (11/21)*(10/21) + (10/21)*(11/21) Prob_Product_Even := (11/21)*(11/21) + (11/21)*(10/21) + (10/21)*(11/21); # Then the determinant is even if #    (1) either the products of the diagonal elements AND  of the off-diagonal #        elements are both even, #    (2) or if they are both odd. # # Thus Prob_Dterminant_Even := Prob_Product_Even^2 + (1-Prob_Product_Even)^2; evalf(%);
 (11)
 > # Visualization of the products of, let's say, diagonal elements. interface(rtablesize=21): s2 := s.s^+; Proportion_of_even := numelems( Select(t -> is(t::even), convert(s2, Vector)) ) / 21^2;
 (12)
 > # Even if the determinants are in [-200, 200], some values do not exist. # # First remark that the 21^2 = 441 products of integers in the range -10..10 are not 441 # different numbers sor some values are more frequent than others. # For instance, the value 0 appears 21+20 = 41 times. # # In fact the number of different product is much lower than 441: Number_of_different_products := numelems({entries(s2, nolist)})
 (13)
 > # There are potentially 401 different values the determinants can get (from -200 to +200). # # Let's take a simple example to prove that not allthese 401 values do exist. # It is clear that the product of two integers in the range -10..10 cannot be a number # in the open interval (-99, -91) union (91, 99). # Then, no determinant can have a value in the open range (-199, -191) union (191, 199). # # # It is quite complex to determine analytically the number (andthe values) of all feasible # determinants, but this becomes simple if we counts all of them. # First consider all the values of the product of two integers in the range -10..10 by # using 's2' above: All_values_of_products := {entries(s2, nolist)}; Navp := numelems(All_values_of_products);
 (14)
 > # Then form the (Navp, Navp) matrix whose each column are identical and contain the elements # of 'All_values_of_products'. Mavp := <[All_values_of_products[]]> . Vector[row](Navp, 1);
 (15)
 > # Form the matrix 'Mavp - Mavp^+'. # This matrix contains all the possible values of the determinants. Mdet := Mavp - Mavp^+;
 (16)
 > # Select all the different values in Mdet andcount them. # # Observe these outcomes correspond to those the brute force above gave us. All_Feasible_Determinants := { entries(%, nolist) }; Nafd := numelems(All_Feasible_Determinants);
 (17)
 > # What are the missing values? Missing_Determinants := {\$-200..200} minus All_Feasible_Determinants
 (18)
 > # This mean that there exist no 'RandomMatrix(2^2, generator=-10..10)' whose determinant # has a value in 'Missing_Determinants'.

## @Axel Vogt @one man You're rig...

@Axel Vogt @one man

You're right, it's pretty common, but one might expect all Mapleprimes users to be a step above.

## Not a bug....

The plot with scale x = 1850 .. 2100 is not a zoom of the plot with scale x = 1800 .. 2100 : the date are not the same given the way you buid them.

## @one man What's disturbing is that ...

@one man

What's disturbing is that you sometimes put a lot of effort into providing a complete answer, but after a few back-and-forths, the OP becomes silent, so you never know whether what you've done has been helpful or not.
Sometimes you feel like a sucker.

## @delvin Thanks.Were you able to go ...

Thanks.

Were you able to go forward from file ddd_2_mmcdara.mw ?

## @Carl Love Yiu're right Cerl, I...

@Carl Love

You're right Carl, I should have open a post.
Sorry for the inconvenience.

## @delvin Here is a correction up to ...

Here is a correction up to the evaluation of fin1.
Check it aand try to a step ahead.

ddd_2_mmcdara.mw

## I still haven't heard from you from abou...

Updated June 20, GMT 7pm, have you suddenly gone mute?

Can't you run my code after modifying the parameter values as described in the last few lines?
Do you still need help?

• Do you have any idea of alpha AP and aAlpha S upper bounds?

• What were the initial conditions for the experimental data you provide?

@delvin
I propose you a rewritting of your ddd.mw file.
I modified it up to the command

`fin1 := ...`

I have pointed out an inconsistency in the equations (a function depends on xi[n] and you take its derivative wrt t).
I corrected this in a way which seems correct to me, but I ask you to validate what i did.

If it is ok, one could go further following a step by step way.

ddd_1_mmcdara.mw

## Please, do not use the email,...

for I will not answer it.
Keep doing exchanges within this thread.

TIA

## @sursumCorda As the new parameters ...

@sursumCorda

As the new parameters a and b are of very different orders of magnitude, a second transformation b -> 1/beta helps finding values closer to three (error about 1e-6 to 3e-6).
With this transformation  a and b then have relatively closed values and a rand()() inititialpoint seems to become more "robust".

This transformation comes from this observation: if one chooses an initial point such that a=r, b=1/r, then the solution is extremely close to a=r, b=1/r with a value of the objective function closed to 3.

 > restart:
 > randomize()
 (1)
 > fun := ((2*(x+y+z))*(sqrt(y*z*(z+x)*(x+y))/(z+2*x+y)+sqrt(z*x*(x+y)*(y+z))/(x+2*y+z)+sqrt(x*y*(y+z)*(z+x))/(y+2*z+x))-9*x*y*z/(x+y+z)-2*(x*y+x*z+y*z))/(sqrt(x*y*z/(x+y+z))*(x+y+z-sqrt(27*x*y*z/(x+y+z))))
 (2)
 > J := eval(fun, [y=a*x, z=b*x]): J := normal(J) assuming x > 0: indets(%, name); K := simplify(J, size): simplify(eval(K, b=1/beta), size): opt := Optimization['Maximize'](%, assume = nonnegative, initialpoint = ({seq})(w = (evalf(rand()())), w in [a, beta])); printf("%1.3e", 3-opt[1])
 1.04e-06
 > J := eval(fun, [x=a*y, z=b*y]): J := normal(J) assuming y > 0: indets(%, name); K := simplify(J, size): simplify(eval(K, b=1/beta), size): opt := Optimization['Maximize'](%, assume = nonnegative, initialpoint = ({seq})(w = (evalf(rand()())), w in [a, beta])); printf("%1.3e", 3-opt[1])
 1.774e-06
 > J := eval(fun, [x=a*z, y=b*z]): J := normal(J) assuming z > 0: indets(%, name); K := simplify(J, size): simplify(eval(K, b=1/beta), size): opt := Optimization['Maximize'](%, assume = nonnegative, initialpoint = ({seq})(w = (evalf(rand()())), w in [a, beta])); printf("%1.3e", 3-opt[1])
 1.035e-06
 >

What initial value of a can we use as a good guess?
The answer is given by the following result

```J := eval(fun, [y=a*x, z=x/a]):
K := simplify(normal(J), size) assuming x::real, x > 0:
asympt(K, a, 3)
(1/2)
/1\         /1\
3 - |-|      + O|-|
\a/         \a/
```

which suggests that the higher the value of a (together with beta=1/b=a) the smaller the error.

A few numerical investigations lead to infer this result:

Led D the value of Digits, (D "large") and set the initial value of a to 10D/4, then the error is 10D/8 :

```Digits:=400:

simplify(eval(K, b=1/beta), size):
simplify(eval(%, beta=a), size) assuming a > 0:
opt := Optimization:-Maximize(%, initialpoint={a=10^100}, iterationlimit=10^3):
printf("%1.3e", 3-opt[1]);
1.000e-50

```

## For information...

```restart
interface(version)
Standard Worksheet Interface, Maple 2015.2, Mac OS X, ...
solve({a > 0, ln(a*(1 + a)) >= 0}, a);

allvalues(%);
/      /  2                    \     \
{ RootOf\_Z  + _Z - 1, index = 1/ <= a }
\                                    /
/1  (1/2)   1     \
{ - 5      - - <= a }
\2          2     /

```

## I agree with @Rouben Rostamian : s...

I agree with @Rouben Rostamian : some clarification is needed.

A few remarks:

• Do not load useless packages, this is the best way to generate conflicts.
• Do not try to use odeplot if you didn't use successfuly dsolve/numeric before.
• If you use dsolve/numeric, you must provide initial and/or boundary conditions.
• Before using a package, CurveFitting for instance, look at the corresponding help page to see what procedures it contains instead of imagining them (NonlinearFit is part of the Statistics package).

You will find in the attached file some elements which lead to a solution of the differential system.
I took some liberties in choosing arbitrary the initial conditions and the values of the parameters (which I all rewrote in, IMO, a more convenient form).

I discuss your "fit to data problem", but this is not conclusive at all: you need to be more precise in your request.

@Rouben Rostamian : I hope I didn't step on your toes.

## No error of your own in fsh.mw...

@delvin

but it"s likely there are som conflicts:

• between 2D input mode and 1D input more,
• the Maple versions (I use Maple 2015).
When I open your file I get a warning message, something like " this code has been written in a newer version and could not run correctly".
If I agree to "activate" the worksheet, then when I save it under another name, the m file is a Maple 2015 mw file.
When you will open it you won't get any warning on your side. But the fact is hat there is also a version mixing between 2023 and 2015.

Here is fsh.mw corrected, at least in the sense that the error is removed.
fsh_eq02.mw

 2 3 4 5 6 7 8 Last Page 4 of 113
﻿