Advanced Search

mmcdara

7891 Reputation

22 Badges

9 years, 59 days
Contact mmcdara

MaplePrimes Activity


These are replies submitted by mmcdara

@acer Hanks for the tip.Have a good...

mmcdara 7891
August 06 2019
0 0

@acer 

Thanks for the tip.

In Maple 2015 
v := Distribution(PDF = (t -> piecewise(0 <= t and t < 1, 1, 0)), ':-Mean'=()->1/Pi):
returns an error message : Error, `->` unexpected


We'll see this later, have a good week

@acer  "Your problematic code...

mmcdara 7891
August 06 2019
0 0

@acer 

"Your problematic code does not ...."
Sorry acer, I've just sent a reply to Carl Love where I waid I've done some typos the code I presented.

In  this same reply there is a "typo-free example" of what I really want to do.
It will help you to understand what is my real poblem.



 

@Carl Love Concerning our 5 advices...

mmcdara 7891
August 06 2019
0 0

@Carl Love 

Concerning our 5 advices:

  1. Thanks, I will remember.
  2. You're right, it was a typo.
  3. It was a typo too, but using ParentName doesn't change the point: the same error is generated.
  4. I've read them but I I will read them again more carefully.
  5. "avoid compound consitions ...": another good advice, thank you



"...That is your goal, correct?"
One hundred percent, you put your finger on it!

I'm presently re-writting my own goodness-of-fit tests, in particular a ChiSquare one, and I would like to have the sample S and the target distribution T as argument.
In ChiSquareSuitableModelTest these arguments are S and something like Normal(0, 1).

Of course, beyond some loss of abstraction, there is no difference between 
T := RandomVariable(Normal(0, 1)) : GoFTest(S, T, ...) ; 
and 
ChiSquareSuitableModelTest(S, Normal(0, 1), ...).

(I will show you later another situation where this difference matters)

To be able to realize the calculations within the procedure GoFTest, I must extract two informations from T: its ParentName, and its Parameters.
More of this I want to be capable to handle situations where T can be Normal(0, 1), Normal(0, b), Normal(a, 1), Normal(a, b): if I can extract the attribute 'Parameters' from T I will be able to detect which of them are formal (and need to be fitted from S) and which of them are numerics and correspond to prior requirements (a thing ChiSquareSuitableModelTest does not do and replaces by requesting the number of parameters to fit...which then goes wrong [one of my previous unanswerd question]).

I know how to do this for any type of distributions Maple knows.
Now a slightly different problem.
Let us suppose that S has ben drawn from a population that have  a generalized Beta distribution T. For instance:
T := shift + scale * RandomVariable(BetaDistribution(a, b));
There is no problem to sample T in Maple. 
The problem is that T, defined this way, is not a random variable (just type attributes(T)) and thus I can't used it as a parameter of GofTest.

The code below will give you an hint of what I'd want to do.

> 

restart:
> 

with(Statistics):
> 

randomize():

 

> 

N := 100:


SamplingDistribution := 4 + 2 * RandomVariable(BetaDistribution(3, 3));

S := Sample(SamplingDistribution, N):
Mean(S);

4+2*_R

  

HFloat(4.993588162631437)

 (1)
> 

Chi2GoF := proc(S, Law)

local N, NC:
local Law_p, Law_Parent, Law_Pars, Law_ParsNum, Law_ParsToFit, Law_NbParsToFit, Prior_Pars:
local MomentEquations, ik, ir, Law_Fit, TargetLaw:
local S_min, S_max, S_gap, S_bins, omega, OBS, TargetExc, TargetStat, Target_p_value:
local n, LocalChiSquare, mag, format:

#DEBUG():

N  := numelems(S):
NC := floor(N/5):
NC := min(NC, floor(evalf(N^(3/5))));

Law_p           := [attributes(Law)][3];
Law_Parent      := Law_p:-ParentName:
Law_Pars        := Law_p:-Parameters;
Law_ParsNum     := numelems(Law_p:-Parameters);
Law_ParsToFit   := map(n -> if type(n, {symbol, indexed}) then n end if, Law_Pars);
Law_NbParsToFit := numelems(Law_ParsToFit);


Prior_Pars := copy(Law_Pars):

if Law_NbParsToFit > 0 then
  MomentEquations := NULL:
  ik := 1:
  for ir from 1 to Law_NbParsToFit do
     if indets(Moment(Law, ik), name) = {} then
        ik := ik+1;
        ir := ir-1:
     else
        MomentEquations := MomentEquations, (Moment(Law, ik) = Moment(S, ik)):
        ik := ik+1:
     end if:
  end do:
  MomentEquations := [MomentEquations];

  Law_Fit   := solve( MomentEquations, Law_ParsToFit )[1];
  TargetLaw := Specialize(Law, Law_Fit):

elif Law_NbParsToFit = 0 then
  TargetLaw := 1*Law
end if:


S_min  := min(S):
S_max  := max(S):
S_gap  := (S_max - S_min)/NC:
S_bins := [S_min-S_gap/2, seq(Quantile(TargetLaw, n/NC, numeric), n=1..NC-1), S_min+S_gap/2]:
omega  := zip((u,v) -> u..v, %[1..-2], %[2..-1]);

OBS := rhs~( TallyInto(S, omega, bins=1) ):

TargetExc := [(N/NC)$NC]:
OBS       := rhs~( TallyInto(S, omega, bins=1) );

# < Vector[column](2, [Oberved, Posterior]) | Matrix(2, NC, [OBS, evalf[3](TargetExc)]) >;

TargetStat := add((OBS -~ TargetExc)^~2 /~ TargetExc):

Target_p_value := Probability(RandomVariable(ChiSquare(NC-1-Law_NbParsToFit)) > TargetStat, numeric);

Law_p      := [attributes(TargetLaw)][3];
Law_Parent := Law_p:-ParentName:
Law_Pars   := Law_p:-Parameters;
# Law_Pars   := seq(cat(Law_Pars[k], ","), k=1..numelems(Law_Pars)-1), Law_Pars[numelems(Law_Pars)];


if Law_NbParsToFit > 0 then
   printf("Prior shape distribution               = %a\n", Law_Parent(op(Prior_Pars)));
   printf("Target distribution fitted from sample = %a\n\n", Law_Parent(op(Law_Pars)));
else
   printf("Prior distribution = %a\n", Law_Parent(op(Prior_Pars)));
end if:
printf("%s\n", "____________________________________________________________"):
printf("%s\n", "          bins            | Expected | Observed |   CHI^2  |"):
printf("%s\n", "--------------------------|----------|----------|----------|"):
for n from 1 to numelems(omega) do
   LocalChiSquare := (OBS[n]-evalf(TargetExc[n]))^2 / evalf(TargetExc[n]):
   mag            := floor(log[10](LocalChiSquare))+1:
   format         := cat(
                          "%1.4e .. %1.4e  |  %6.2f  |  %6d  |  ",
                           piecewise(
                                      LocalChiSquare > 1,
                                      cat(seq(" ", k=1..4-mag), "%", mag, ".2f |\n"),
                                      "   %0.2f |\n"
                                    )
                        );
   printf(
           format,
           op(omega[n]),
           evalf(TargetExc[n]),
           OBS[n],
           LocalChiSquare
         );
end do:
printf("%s\n", "--------------------------|----------|----------|----------|"):
mag    := floor(log[10](TargetStat))+1:
format := cat(
               "%1.4e .. %1.4e  | %6d   |  %6d  |  ",
                piecewise(
                           TargetStat  > 1,
                           cat(seq(" ", k=1..4-mag), "%", mag, ".2f |\n"),
                           "   %0.2f |\n"
                         )
             );
printf(
        format,
        op(1,omega[1]), op(2,omega[-1]),
        N,
        N,
        TargetStat
      );
printf("%s\n", "--------------------------|----------|----------|----------|"):
format := cat(seq(" ", k=1..37), "|  p.value |   %0.4f |\n"):
printf(format, Target_p_value);
format := cat(seq(" ", k=1..37), "-----------------------\n"):
format := cat(seq(" ", k=1..37), "|__________|__________|\n"):
printf(format):

return [omega, evalf(TargetExc), OBS, cat(convert(Law_Parent, string), "(", Law_Pars, ")"), Target_p_value]

end proc:
 
> 

Target_1 := RandomVariable(Normal(a,b)):
res := Chi2GoF(S, Target_1):

Prior shape distribution               = NormalDistribution(a,b)
Target distribution fitted from sample = NormalDistribution(4.993588163,.3724480077)

____________________________________________________________
          bins            | Expected | Observed |   CHI^2  |
--------------------------|----------|----------|----------|
4.2164e+00 .. 4.4345e+00  |    6.67  |       6  |     0.07 |
4.4345e+00 .. 4.5799e+00  |    6.67  |      12  |     4.27 |
4.5799e+00 .. 4.6801e+00  |    6.67  |       5  |     0.42 |
4.6801e+00 .. 4.7616e+00  |    6.67  |       5  |     0.42 |
4.7616e+00 .. 4.8332e+00  |    6.67  |       5  |     0.42 |
4.8332e+00 .. 4.8992e+00  |    6.67  |      10  |     1.67 |
4.8992e+00 .. 4.9624e+00  |    6.67  |       6  |     0.07 |
4.9624e+00 .. 5.0247e+00  |    6.67  |       5  |     0.42 |
5.0247e+00 .. 5.0879e+00  |    6.67  |       7  |     0.02 |
5.0879e+00 .. 5.1540e+00  |    6.67  |       6  |     0.07 |
5.1540e+00 .. 5.2256e+00  |    6.67  |       4  |     1.07 |
5.2256e+00 .. 5.3070e+00  |    6.67  |       5  |     0.42 |
5.3070e+00 .. 5.4073e+00  |    6.67  |      11  |     2.82 |
5.4073e+00 .. 5.5527e+00  |    6.67  |       7  |     0.02 |
5.5527e+00 .. 4.3259e+00  |    6.67  |       0  |     6.67 |
--------------------------|----------|----------|----------|
4.2164e+00 .. 4.3259e+00  |    100   |     100  |    18.80 |
--------------------------|----------|----------|----------|
                                     |  p.value |   0.0935 |
                                     |__________|__________|

  
> 

Target_2 := RandomVariable(Normal(5, 0.5)):
res := Chi2GoF(S, Target_2):

Prior distribution = NormalDistribution(5,.5)
____________________________________________________________
          bins            | Expected | Observed |   CHI^2  |
--------------------------|----------|----------|----------|
4.2164e+00 .. 4.2495e+00  |    6.67  |       0  |     6.67 |
4.2495e+00 .. 4.4446e+00  |    6.67  |       6  |     0.07 |
4.4446e+00 .. 4.5792e+00  |    6.67  |      11  |     2.82 |
4.5792e+00 .. 4.6885e+00  |    6.67  |       7  |     0.02 |
4.6885e+00 .. 4.7846e+00  |    6.67  |       5  |     0.42 |
4.7846e+00 .. 4.8733e+00  |    6.67  |       9  |     0.82 |
4.8733e+00 .. 4.9582e+00  |    6.67  |      10  |     1.67 |
4.9582e+00 .. 5.0418e+00  |    6.67  |       8  |     0.27 |
5.0418e+00 .. 5.1267e+00  |    6.67  |       6  |     0.07 |
5.1267e+00 .. 5.2154e+00  |    6.67  |       9  |     0.82 |
5.2154e+00 .. 5.3115e+00  |    6.67  |       5  |     0.42 |
5.3115e+00 .. 5.4208e+00  |    6.67  |      11  |     2.82 |
5.4208e+00 .. 5.5554e+00  |    6.67  |       7  |     0.02 |
5.5554e+00 .. 5.7505e+00  |    6.67  |       3  |     2.02 |
5.7505e+00 .. 4.3259e+00  |    6.67  |       0  |     6.67 |
--------------------------|----------|----------|----------|
4.2164e+00 .. 4.3259e+00  |    100   |     100  |    25.55 |
--------------------------|----------|----------|----------|
                                     |  p.value |   0.0295 |
                                     |__________|__________|

  
> 

Target_3 := RandomVariable(Normal(5, b)):
res := Chi2GoF(S, Target_3):

Prior shape distribution               = NormalDistribution(5,b)
Target distribution fitted from sample = NormalDistribution(5,.2732036904)

____________________________________________________________
          bins            | Expected | Observed |   CHI^2  |
--------------------------|----------|----------|----------|
4.2164e+00 .. 4.5899e+00  |    6.67  |      18  |    19.27 |
4.5899e+00 .. 4.6965e+00  |    6.67  |       6  |     0.07 |
4.6965e+00 .. 4.7701e+00  |    6.67  |       4  |     1.07 |
4.7701e+00 .. 4.8298e+00  |    6.67  |       4  |     1.07 |
4.8298e+00 .. 4.8823e+00  |    6.67  |       8  |     0.27 |
4.8823e+00 .. 4.9308e+00  |    6.67  |       6  |     0.07 |
4.9308e+00 .. 4.9771e+00  |    6.67  |       5  |     0.42 |
4.9771e+00 .. 5.0229e+00  |    6.67  |       3  |     2.02 |
5.0229e+00 .. 5.0692e+00  |    6.67  |       5  |     0.42 |
5.0692e+00 .. 5.1177e+00  |    6.67  |       3  |     2.02 |
5.1177e+00 .. 5.1702e+00  |    6.67  |       5  |     0.42 |
5.1702e+00 .. 5.2299e+00  |    6.67  |       5  |     0.42 |
5.2299e+00 .. 5.3035e+00  |    6.67  |       4  |     1.07 |
5.3035e+00 .. 5.4101e+00  |    6.67  |      11  |     2.82 |
5.4101e+00 .. 4.3259e+00  |    6.67  |       0  |     6.67 |
--------------------------|----------|----------|----------|
4.2164e+00 .. 4.3259e+00  |    100   |     100  |    38.05 |
--------------------------|----------|----------|----------|
                                     |  p.value |   0.0003 |
                                     |__________|__________|

  
> 

Target_4 := RandomVariable(Normal(a, 0.5)):
res := Chi2GoF(S, Target_4):

Prior shape distribution               = NormalDistribution(a,.5)
Target distribution fitted from sample = NormalDistribution(4.993588163,.5)

____________________________________________________________
          bins            | Expected | Observed |   CHI^2  |
--------------------------|----------|----------|----------|
4.2164e+00 .. 4.2430e+00  |    6.67  |       0  |     6.67 |
4.2430e+00 .. 4.4382e+00  |    6.67  |       6  |     0.07 |
4.4382e+00 .. 4.5728e+00  |    6.67  |      10  |     1.67 |
4.5728e+00 .. 4.6821e+00  |    6.67  |       8  |     0.27 |
4.6821e+00 .. 4.7782e+00  |    6.67  |       4  |     1.07 |
4.7782e+00 .. 4.8669e+00  |    6.67  |       9  |     0.82 |
4.8669e+00 .. 4.9518e+00  |    6.67  |      11  |     2.82 |
4.9518e+00 .. 5.0354e+00  |    6.67  |       7  |     0.02 |
5.0354e+00 .. 5.1203e+00  |    6.67  |       7  |     0.02 |
5.1203e+00 .. 5.2090e+00  |    6.67  |       9  |     0.82 |
5.2090e+00 .. 5.3051e+00  |    6.67  |       5  |     0.42 |
5.3051e+00 .. 5.4144e+00  |    6.67  |      11  |     2.82 |
5.4144e+00 .. 5.5490e+00  |    6.67  |       7  |     0.02 |
5.5490e+00 .. 5.7441e+00  |    6.67  |       3  |     2.02 |
5.7441e+00 .. 4.3259e+00  |    6.67  |       0  |     6.67 |
--------------------------|----------|----------|----------|

4.2164e+00 .. 4.3259e+00  |    100   |     100  |    26.15 |
--------------------------|----------|----------|----------|
                                     |  p.value |   0.0162 |
                                     |__________|__________|

  
> 

Target := shift + scale * RandomVariable(BetaDistribution(a, b)):
Mean(Target);

[attributes(Target)];  # ===> MyTarget is not a random variable

a*scale/(a+b)+shift

  

[]

 (2)
> 

f := PDF(Target, t);
m := Mean(Target);

f := piecewise((t-shift)/scale < 0, 0, (t-shift)/scale < 1, ((t-shift)/scale)^(-1+a)*(1-(t-shift)/scale)^(-1+b)/Beta(a, b), 0)/abs(scale)

  

a*scale/(a+b)+shift

 (3)
> 

# Define a user random variable

v := Distribution(PDF = (t -> f), 'Conditions'=[scale > 0, a > 0, b > 0]);
MyTarget := RandomVariable(v):


A := attributes(MyTarget)[3];
AllAttributes := exports(A);

 

_m4893453824

  

_ProbabilityDistribution12

  

Conditions, PDF

 (4)
> 

# Auxiliary question: why is Maple unabled to compute the Mean?
Mean(MyTarget);


v := Distribution(PDF = (t -> f), 'Conditions'=[scale > 0, a > 0, b > 0], 'Mean'=m);
MyTarget := RandomVariable(v):
Mean(MyTarget);

a*scale/(a+b)+shift

  

_m4929682112

  

a*scale/(a+b)+shift

 (5)
> 

v := Distribution(PDF = (t -> f), 'Conditions'=[scale > 0, a > 0, b > 0], 'Mean'=m, 'ParentName'= GeneralizedBetaDistribution, 'Parameters'=[a, b, scale, shift]);

Error, (in Statistics:-Distribution) invalid input: too many and/or wrong type of arguments passed to NewDistribution; first unused argument is ParentName = GeneralizedBetaDistribution

  
> 

 


 

Download My_ChiSquare_ChiSquare_Gof_Test.mw

@vv  Thank you. I didn't know ...

mmcdara 7891
August 05 2019
0 0

@vv 

Thank you.
I didn't know this function

 

@Christian Wolinski  PS: I didn...

mmcdara 7891
August 05 2019
0 0

@Christian Wolinski 

PS: I didn't give you a "best answer" for this reply but for your previous one
g := subs(F = unapply(f,a,b,x), (a,b) -> (x-> F(a, b, x)))
for which, unfortunately, there was no thumb or cup on the top right of your reply.

@Carl Love  @Christian Wolinski @...

mmcdara 7891
August 05 2019
0 0

@Carl Love 

@Christian Wolinski

@tomleslie

@Kitonum

Thanks to both of you.

 

@Kitonum  Yes, I know, but I reall...

mmcdara 7891
August 05 2019
0 0

@Kitonum 

Yes, I know, but I really want to define a funcyion of the form g(a, b)(x)

 

@ThU  Thank you ...

mmcdara 7891
August 01 2019
0 0

@ThU 

Thank you 

@Carl Love  Hi, I don't underst...

mmcdara 7891
July 31 2019
0 0

@Carl Love 
Hi, I don't understand neither. I just replied this question Expand and negative exponent...

@Carl Love  Given the plot at the ...

mmcdara 7891
July 23 2019
0 0

@Carl Love 

Given the plot at the end of the question, I guess this latter is related to a previous one (Cubic b-spline collocation method...)  and that the magic numbers are the nodal coefficients of a bspline-based collocation method.

@tomleslieHere is a solution, even if I ...

mmcdara 7891
July 23 2019
0 0

@DimaBrt
@tomleslie

Here is a solution, even if I agree with tomleslie (I spent a lot of time to code this stuff, and it's not high level coding).

Other point, if you have been asked by some teacher of yours to do this work, it's likely that he had in mind the writting of a linear system in matrix form  (see for instance https://pdfs.semanticscholar.org/43d0/d5a6f30cdc140647039c0edd5a268138700b.pdf).
Here I cheated by using the procedure solve to find the coefficients of the bsplines in the expansion of u.

It's just raw material, take it as a starting point for you to go further by your own.

bspline.mw

@tomleslie  Just a little correcti...

mmcdara 7891
July 23 2019
0 0

@tomleslie 

Just a little correction: you wrote  4*x*2 instead of  4*x^2.

To answer your question, I suppose DimaBrt has been asked to use the collocation method as part of a homework.

@Carl Love  From my courses or the...

mmcdara 7891
July 21 2019
0 0

@Carl Love 

From my courses or the books I have, the term "central moment" (and even "moment") seems to be everywhere used for random variables only.
This is a special case where the whole population is knoxn, at least in some sense*
Thus the question of biasedness or unbiasadeness is of no matter.
On the other side the term variance is ued both for random variable and samples (in this later case one should say   "empirical variance") and a sample generally not being the whole  population  the question of deriving unbiased estimator occurs.

Beyond the references you gave, I've browsed the web on my side. And what I found confirms you reply.
For instance R contains no function to assess central moments and SPSS doesn't seem to have any either: which seems to indicate that the concept of central moment would be generally restricted to formal statistics(?)
One of the few exception I was capable to find is 

Unbiased Estimation of Central Moments by using U‐statistics

Peter M. Heffernan

First published: 1997

https://doi.org/10.1111/1467-9868.00102

 

So, to conclude, it was rather obvious to me that the central moment of order 2 and the variance being the same thing  when introduced on random variables, their "sampled counterparts"  had to be the same.
But it seems that "empirical 2nd order central moment" is not an usual concept and that "empirical variance" would be the only one used (apart the reference you gave and the one above).
As a consequence the question of unbiasedness doesn't seem to be relevant for central moments(?).
Maybe the good attitude would be to never use central moments on sample...

Thank you for this enriching discussion.


* I mean here that the population, even if not all its outcomes are explicitely enumerated, one has a sufficient characterization of it. 
This is the same difference between characterizing a set by all its elements or by a sufficient proprety (I'm sorry but I do not know the english terms for that ; in french one says "definition en extension" (="definition by extension") and "definition en comprehension" (="definition by understanding").
https://lexique.netmath.ca/definition-dun-ensemble-en-comprehension/
https://lexique.netmath.ca/definition-dun-ensemble-en-extension/

@Carl Love But how is the command s...

mmcdara 7891
July 21 2019
0 0

@Carl Love 

The difference is ..... To me, it seems useful to be able to make that distinction, so I like that the two commands return different result

This is your personal opinion, but, in this case, this difference should be explicitely documented, which is not the case.
My personal opinion leads to a questionable implementation of the command CentralMoment.
I suppose this is a debate that should be solved by expert opinion from probabilists or statisticians.




But how is the command supposed to know whether the data that you've given it is a sample or the whole population?

Right, and the best example of this is given by ... Excel !!!
In Excel there esist 2 variances, one is the "standard" unbiased estimation ( add(...)/(N-1)) and the second is add(...)/N with a comment of the form  "if the sample is  the whome distribution use the command where the sum is divided by N, otherwise use the command with N-1 at the denominator"
This is a very good explanation even if it often leaves the user confused.

So, "how is the command supposed to know whether the data that you've given it is a sample or the whole population?"
The answer is: in light of the Excel's example : "the command can't know".
But, the default standard way to compute the Variance of a sample is generally based on the assumption that the sample doesn't represent the whole population (and this is what Maple's Variance does when applied to a sample).

This is questionable for discrete random variables with bounded support, but perfectly reasonable in any other case (obviously the whole population cannot be explicitely enumerated if its obey an unbounded distribution).

Nevertheless the prolem I raised is not that one but this one: Maple considers that Variance(S) (S being some sample) is to be assessed by assuming S does not represent the whole population, but considers that CentralMoment(S, 2) is to be assessed by assuming S represents the whole population.

This is this unconsistency I pointed out.

@acer  "It might accomplish y...

mmcdara 7891
July 20 2019
0 0

@acer 

"It might accomplish your goal, but it doesn't do precisely what you asked."...absolutely right, its "works" but doesn't correctly answer my question 

First 122 123 124 125 126 127 128 Last Page 124 of 154