Advanced Search

mmcdara

7891 Reputation

22 Badges

9 years, 60 days
Contact mmcdara

MaplePrimes Activity


These are replies submitted by mmcdara

@Carl Love EXCELLENT !I hadn't ...

mmcdara 7891
December 02 2019
0 0

@Carl Love 

EXCELLENT !

I hadn't have the energy to do this, BRAVO!
 

@acer  With some delay, thank you ...

mmcdara 7891
November 30 2019
0 0

@acer 

With some delay, thank you nevertheless for your clarifications concerning the mechanisms behind interface.

@Carl Love  Thanks for correcting ...

mmcdara 7891
November 30 2019
0 0

@Carl Love 

Thanks for correcting the typo.
The  p-value is indeed the "good" criterion,  the goal xas here to emphasoze how far from the expected values are the observations.

@acer  (this reply from my persona...

mmcdara 7891
November 29 2019
0 0

@acer 

(this reply from my personal login)

Thanks.
Yes I had already noticed this "plot inheritance" you refer to.
By the way it also seems to be the case for other functions: for instance interface(displayprecision = some number) seems "persistent" despite a restart (curiously if I insert the contents of the attached file, Pi is displayed with 10 digits).
 

Download persistence.mw

@Carl Love  @acer (This reply from...

mmcdara 7891
November 29 2019
0 0

@Carl Love 
@acer

(This reply from my personal login)
Thank you both.
I have to say that I was a little cautious about taking the step to buy Maple 2019.

@Carl Love  Here is the correct ve...

mmcdara 7891
November 29 2019
0 0

@Carl Love 

Here is the correct version of the evaluation of the sampling procedure for a few different continuous distributions.
Strangely the only awfull result is obtained for the Normal distribution.

The first column of printf, titled "prob" denotes a probability value.
The second one denotes the expected number of samples that should be less than the "prob" quantile.
The third one indicates how many samples from S are less than this "prob" quantile.
The last one expresses the differences between columns 2 and 3 in terms of the theoritical 
​​​​​​​standard deviation of the observations form the expected values.

> 

restart:
> 

interface(version)

`Standard Worksheet Interface, Maple 2015.2, Mac OS X, December 21 2015 Build ID 1097895`

 (1)
> 

# this procedure is inspired by carl love's original procedure

SampleCheck:= proc(X, N::posint)
local on, infbound, p, q, B, i, S, Obs, Exc, U;
  uses Statistics;
  on  := convert(Support(X), RealRange);
  if on::function then
     infbound := map(u -> if is(type(u, function)) then op(u) else u end if, [op(on)])[1];
     if infbound = -infinity then infbound := -10^10 end if;
  else
    infbound := -10^10;
  end if;
  p   := evalf(10^~[$-6..-1]);
  q   := Quantile~(X, p, numeric);
  S   := Sample(X, N):
  B   := map(u -> infbound..u, q);
  Obs := < map(b -> op(TallyInto(S, b, bins=1)), B) >:
  Exc := < N*~p >:

  # The term sqrt~(Exc) is the standard deviation of a binomial rv of parameters (N, p)
  #
  # The term abs~(Exc -~ rhs~(Obs)) is the absolute deviation of Obs from Exc.
  #
  # The ratio of this later term by the former says how far is Obs from Exc in terms of standard
  # deviations sqrt~(Exc). 

  U := abs~(Exc -~ rhs~(Obs)) /~ sqrt~(Exc) ;
  printf("  prob      exact   observed    reduced uncertainty \n");
  for i from 1 to numelems(p) do
    if U[i] > 6 then
      printf("%1.1e  %8d %11d           %2.2f ****\n", p[i], round(Exc[i]), rhs(Obs[i]), U[i]);
    else
      printf("%1.1e  %8d %11d           %2.2f\n", p[i], round(Exc[i]), rhs(Obs[i]), U[i]);
    end if:
  end do:
end proc:
 
> 

ObsExp := SampleCheck(Statistics:-RandomVariable(Normal(0, 1)), 10^8);

  prob      exact   observed    reduced uncertainty

1.0e-06       100         413           31.30 ****
1.0e-05      1000        2134           35.86 ****
1.0e-04     10000       11616           16.16 ****
1.0e-03    100000      100236           0.75
1.0e-02   1000000      999482           0.52
1.0e-01  10000000    10000876           0.28

  
> 

ObsExp := SampleCheck(Statistics:-RandomVariable(Exponential(1)), 10^8);

  prob      exact   observed    reduced uncertainty

1.0e-06       100          96           0.40
1.0e-05      1000        1004           0.13
1.0e-04     10000       10028           0.28
1.0e-03    100000       99744           0.81
1.0e-02   1000000     1000610           0.61
1.0e-01  10000000    10002152           0.68

  
> 

ObsExp := SampleCheck(Statistics:-RandomVariable(BetaDistribution(5, 0.2)), 10^8);

  prob      exact   observed    reduced uncertainty

1.0e-06       100         112           1.20
1.0e-05      1000         996           0.13
1.0e-04     10000        9863           1.37
1.0e-03    100000       99561           1.39
1.0e-02   1000000      998378           1.62
1.0e-01  10000000     9995881           1.30

  
> 

ObsExp := SampleCheck(Statistics:-RandomVariable(Weibull(1, 10)), 10^8);

  prob      exact   observed    reduced uncertainty

1.0e-06       100          87           1.30
1.0e-05      1000         997           0.09
1.0e-04     10000        9988           0.12
1.0e-03    100000       99792           0.66
1.0e-02   1000000      999710           0.29
1.0e-01  10000000    10004013           1.27

  
> 

 

> 

 


 

Download Bug_Normal_Sampler_Continued_2.mw

@Carl Love  Right, after a good nig...

mmcdara 7891
November 29 2019
0 0

@Carl Love 
Right, after a good night I realised I was completely stupid.
I'll fix it after work, let's say in 10 hours.
Have a good day and thanks for your critical eye.

For the moment I remove the post remove the post 

@acer Hi acer, I didn't ev...

mmcdara 7891
November 28 2019
0 0

@acer 
Hi acer, 
I didn't even have time to test your BoxMuller_ac1 procedure that you have done all the job.
The preformances you get are impressive but I can't reproduce them (still this Compiler problem : I use a Maple 2015.2 personal licence on my private laptop).
... but I trust you


By the way, you're perfectly right : "The OP's BoxMuller_2 procedure seems reasonably quick on a few runs, but produces more collectible garbage. This costs, if repeated more time"
 

> 

restart;
> 

BoxMuller_ac2 := module()
> 

  local helper, ModuleApply, ModuleLoad, compiled;
> 

  ModuleApply:=proc(N)
> 

    local c, n, R, V, X1, X2:
> 

    n := ceil(2*N/3);
> 

    V := Statistics:-Sample(Uniform(-1, 1), [n,2]):
> 

    R := Vector(2*n,datatype=float[8]);
> 

    if compiled then
> 

      c := helper(V,R,n);
> 

    else
> 

      c := evalhf(helper(V,R,n));
> 

    end if;
> 

    return R[1..2*trunc(c)];
> 

  end proc:
> 

  helper:=proc(VV,R,n)
> 

    local c, i, SSi;
> 

    c := 0;
> 

    for i from 1 to n do
> 

      SSi := VV[i,1]^2 + VV[i,2]^2;
> 

      if SSi <= 1.0 then
> 

        c := c+1;
> 

        SSi := sqrt(-2*log(SSi)/SSi);
> 

        R[2*c-1] := SSi*VV[i,1];
> 

        R[2*c] := SSi*VV[i,2];
> 

      end if;
> 

    end do;
> 

    return c;
> 

  end proc;
> 

  ModuleLoad := proc()
> 

    try
> 

      helper := Compiler:-Compile(helper);
> 

      compiled := true;
> 

    catch:
> 

      compiled := false;
> 

      WARNING("will use evalhf mode, not compiled");
> 

    end try;
> 

  end proc:
> 

  ModuleLoad();
> 

end module:

Warning, will use evalhf mode, not compiled

  
> 

Boutac2:=CodeTools:-Usage(BoxMuller_ac2(10^5), iterations=10):
> 

Boutac2:=CodeTools:-Usage(BoxMuller_ac2(10^5)):

memory used=2.90MiB, alloc change=25.20MiB, cpu time=65.50ms, real time=65.70ms, gc time=0ns
memory used=2.85MiB, alloc change=2.84MiB, cpu time=64.00ms, real time=64.00ms, gc time=0ns

  
> 

f := x -> x^2:
ccf := Compiler:-Compile(f)

Error, (in Compiler:-Compile) possible installation problem: GNU C compiler (gcc) not found in your command search path (PATH). You will need to restart Maple after ensuring that gcc is installed and adjusting your PATH environment variable.

  
> 

 


 

Download test_acer_module.mw

@acer  Thank you acer, I'm goi...

mmcdara 7891
November 28 2019
0 0

@acer 

Thank you acer, I'm going to test that in a couple of hours.

@acer  I'd forgotten to thank ...

mmcdara 7891
November 28 2019
0 0

@acer 

I'd forgotten to thank you for this last reply.
My omission is now corrected.

PS: my question about Compiler was related to a question I've asked about 2 years ago to MapleSoft France.
By this time I was interested in compiling the output procedure dsolve[numeric] returns.
It's answer was that the Compiler was no longer delivered with the Maple's licence but that it was always possible to have it for free in a separate delivery.  But it's possible I misunderstood MapleSoft's answer.

@Carl Love  Thanks carl, I've ...

mmcdara 7891
November 28 2019
0 0

@Carl Love 

Thanks carl, I've read it and test it.
Your procedure i indeed very  fast.

@Carl Love  I'm sorry, but I d...

mmcdara 7891
November 28 2019
0 0

@Carl Love 

I'm sorry, but I didn't see your answer by the time I sent replies to acer and epostma.
It simply was absent when I open the thread from the main "questions" page.
I'm going to look at it right now.

@epostma  Thanks for your detailed...

mmcdara 7891
November 27 2019
0 0

@epostma 

Thanks for your detailed answer.

I know a little bit about the ziggurat method but I thought it was essentially used to sample distributions for which there was no ad hoc algorithm.
I suppose that the Beta, Gamma, Cauchy, etc distributions for which very efficient ad hoc algorithms do exist are also simulated from the ziggurat method ?

 

@acer  Thanks acer for all these i...

mmcdara 7891
November 27 2019
0 0

@acer 

Thanks acer for all these interesting comparisons.

I understand that the compiler no longer comes with Maple since version 2018 ... am I right?

I have another question: I initially used  numelems(Select(x -> is(x > q), S) ) to fill matrix M but it appeared to be extremely low. This surprises me for I thought that such build-in procedures were designed to be very fast.
So, what could be the point of using Select instead of add(Heaviside~(...)))?

Thanks again.

The expression  a*x+by+c = 0 i...

mmcdara 7891
November 23 2019
0 0

The expression  a*x+by+c = 0 is more general than y = m*x+p (not b!).
So you can't, in general, transform the first into the second one: for this to be possible b must be non null.
(think to the "vertical" line of equation x-1=0)

First 117 118 119 120 121 122 123 Last Page 119 of 154