Advanced Search

mmcdara

7891 Reputation

22 Badges

9 years, 59 days
Contact mmcdara

MaplePrimes Activity


These are replies submitted by mmcdara

@Carl Love Here it is.This new file...

mmcdara 7891
August 17 2019

@Carl Love 

Here it is.

This new file confirms what I sent you previously.
Still NormInv doesn't work for probabilities 1-1e-7

In some sense it's not dramatic for I will retrieve Maple 2018 in ten days after my vacations.

 

> 

restart:
> 

Digits:= 20: #Changing this doesn't change the code; it's just for testing.
> 

 

> 

#Symbolic steps are only needed to create the numeric procedures, not
> 

#needed to run them:
> 

Ncdf:= unapply(int(exp(-x^2/2)/sqrt(2*Pi), x= -infinity..X), X);
> 

Newton:= unapply(simplify(x-(Ncdf(x)-Q)/D(Ncdf)(x)), (x,Q));
> 

 

> 

#Purely numeric procedures:
> 

NormInv_inner:= proc(Q, eps)
> 

local z:= 0, z1;
> 

   do  z1:= evalf(Newton(z,Q));  if abs(z1-z) < eps then return z1 fi;  z:= z1  od
> 

end proc:
> 

      
> 

NormInv:= (Q::And(numeric, realcons, satisfies(Q-> 0 <= Q and Q <= 1)))->
> 

   if Q > 1/2 then -thisproc(1-Q)
> 

   elif Q = 0 then -infinity
> 

   elif Digits <= evalhf(Digits) then evalhf(NormInv_inner(Q, DBL_EPSILON))
> 

   else NormInv_inner(evalf(Q), 10.^(1-Digits))
> 

   fi
> 

:
> 

 

proc (X) options operator, arrow; 1/2+(1/2)*erf((1/2)*2^(1/2)*X) end proc

  

proc (x, Q) options operator, arrow; -(1/2)*2^(1/2)*Pi^(1/2)*erf((1/2)*2^(1/2)*x)*exp((1/2)*x^2)+2^(1/2)*Pi^(1/2)*Q*exp((1/2)*x^2)+x-(1/2)*2^(1/2)*Pi^(1/2)*exp((1/2)*x^2) end proc

 (1)
> 

kernelopts(version)

`Maple 2015.2, APPLE UNIVERSAL OSX, Dec 20 2015, Build ID 1097895`

 (2)
> 

Seed:=randomize();

156606905209673

 (3)
> 

with(Statistics):
> 

# This restricted probability range is chosen because the full range
# generates situations that do not end in Maple 2015

P := Sample(Uniform(0.4, 0.6), 10^4):


f := proc(P)
local n, q:
for n from 1 to numelems(P) do q := P[n]: NormInv(q) end do:
end proc:
CodeTools:-Usage(Quantile~(Normal(0, 1), P, numeric)):
CodeTools:-Usage(f(P)):

memory used=1.59GiB, alloc change=428.01MiB, cpu time=16.34s, real time=15.85s, gc time=1.25s

memory used=1.02GiB, alloc change=32.00MiB, cpu time=9.59s, real time=8.98s, gc time=917.58ms

  
> 

 


 

Download NotSoFast_2015_2.mw

not so fast an approximation...

mmcdara 7891
August 17 2019

@Carl Love 

Here are the perfomances of NormInv compared to Quantile: you'll see that NormInv is not a very fast alternative to Quantile.
So, maybe it's "hard to believe that any of these "magic coefficients" methods are generally better than this simple Newton's method" but the evidence seems in the numbers.
Maybe it's not for nothing that mathematicians are seeking for fast methods to estimate the normal quantiles (see the paper cited by Axel in the first reply)?

 

> 

# This restricted probability range is chosen because the full range
# generates situations that do not end in Maple 2015

P := Sample(Uniform(0.4, 0.6), 10^4):


f := proc(P)
local n, q:
for n from 1 to numelems(P) do q := P[n]: NormInv(q) end do:
end proc:
CodeTools:-Usage(Quantile~(Normal(0, 1), P, numeric)):
CodeTools:-Usage(f(P)):
 

memory used=1.57GiB, alloc change=32.00MiB, cpu time=16.02s, real time=14.87s, gc time=1.62s

memory used=1.02GiB, alloc change=32.00MiB, cpu time=9.70s, real time=8.80s, gc time=1.29s

  
> 

 


 

Download NotSoFast.mw

@Carl Love  Thanks. I will test yo...

mmcdara 7891
August 16 2019

@Carl Love 

Thanks.
I will test your code once my holidays are over (I will then have Maple 2018 and it doesn(t work in Maple 2015)

For the moment here are the "exact quantiles". 

BTW, my post has followed a strange direction for its purpose was not really to find a fast method to estimate the inverse normal cdf, but rather to point a "mistak" inthis reference book that is Abramowitz & Stegun
 

> 

restart:
> 

Digits := 20:
P := [seq(1 - 10.^(-k), k=1..Digits)]:
> 

print~(P, parse~(sprintf~("%2.5f", Statistics:-Quantile~(Normal(0, 1), P)))):

.90000000000000000000, 1.28155

  

.99000000000000000000, 2.32635

  

.99900000000000000000, 3.09023

  

.99990000000000000000, 3.71902

  

.99999000000000000000, 4.26489

  

.99999900000000000000, 4.75342

  

.99999990000000000000, 5.19934

  

.99999999000000000000, 5.61200

  

.99999999900000000000, 5.99781

  

.99999999990000000000, 6.36134

  

.99999999999000000000, 6.70602

  

.99999999999900000000, 7.03448

  

.99999999999990000000, 7.34880

  

.99999999999999000000, 7.65063

  

.99999999999999900000, 7.94135

  

.99999999999999990000, 8.22208

  

.99999999999999999000, 8.49379

  

.99999999999999999900, 8.75729

  

.99999999999999999990, 9.01327

  

.99999999999999999999, 9.26234

 (1)
> 

 


 

Download ExactQuantiles.mw

@vv Your trick helps goinng a littl...

mmcdara 7891
August 16 2019

@vv 

Your trick helps goinng a little bit further but not beyond P=0.998.

Could you please provide a result obtained with Maple 2018 or 2019?
For instance withe the sequence P=[seq(1-10^k, k=1..6)]

Thanks in advance


 

> 

restart:
> 

Digits:= 15: #Changing this doesn't change the code; it's just for testing.
> 

 

> 

#Symbolic steps are only needed to create the numeric procedures, not
> 

#needed to run them:
> 

Ncdf:= unapply(int(exp(-x^2/2)/sqrt(2*Pi), x= -infinity..X), X);
> 

Newton:= unapply(simplify(x-(Ncdf(x)-Q)/D(Ncdf)(x)), (x,Q));
> 

 

> 

#Purely numeric procedures:
> 

NormInv_inner:= proc(Q)
> 

local z:= 0, z1;
> 

   do  z1:= Newton(z,Q);  if abs(z1-z)<10^(-Digits+1) then return z1 fi;  z:= z1  od
> 

end proc:
> 

      
> 

NormInv:= (Q::And(realcons, satisfies(Q-> 0 <= Q and Q <= 1)))->
> 

   if Q > 1/2 then -thisproc(1-Q)
> 

   elif Q=0 then -infinity
> 

   elif Digits <= evalhf(Digits) then evalhf(NormInv_inner(Q))
> 

   else NormInv_inner(evalf(Q))
> 

   fi
> 

:
> 

 

proc (X) options operator, arrow; 1/2+(1/2)*erf((1/2)*2^(1/2)*X) end proc

  

proc (x, Q) options operator, arrow; -(1/2)*2^(1/2)*Pi^(1/2)*erf((1/2)*2^(1/2)*x)*exp((1/2)*x^2)+2^(1/2)*Pi^(1/2)*Q*exp((1/2)*x^2)+x-(1/2)*2^(1/2)*Pi^(1/2)*exp((1/2)*x^2) end proc

 (1)
> 

interface(version)

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

 (2)
> 

P := [seq(k/(k+1), k=100..1000, 100)]:

for n from 1 to numelems(P) do
  printf("%2d  %0.5f", n, P[n]);
  printf("   %2.5f", NormInv(P[n]));
  print():
end do;

 1  0.99010   2.33008

  

 

 2  0.99502   2.57755

  

 

 3  0.99668   2.71415

  

 

 4  0.99751   2.80784

  

 

 5  0.99800   2.87879

  

 

 6  0.99834

Warning,  computation interrupted

  
> 

 


 

Download NewtonFails_2.mw

@Carl Love  Procedure NormInv fail...

mmcdara 7891
August 16 2019

@Carl Love 

Procedure NormInv fails do return a esult in a reasonable time for Q > 0.969:

 

> 

restart:
> 

Digits:= 15: #Changing this doesn't change the code; it's just for testing.
> 

 

> 

#Symbolic steps are only needed to create the numeric procedures, not
> 

#needed to run them:
> 

Ncdf:= unapply(int(exp(-x^2/2)/sqrt(2*Pi), x= -infinity..X), X);
> 

Newton:= unapply(simplify(x-(Ncdf(x)-Q)/D(Ncdf)(x)), (x,Q));
> 

 

> 

#Purely numeric procedures:
> 

NormInv_inner:= proc(Q)
> 

local z:= 0, z1;
> 

   do  z1:= Newton(z,Q);  if z1=z then return z1 fi;  z:= z1  od
> 

end proc:
> 

      
> 

NormInv:= (Q::And(realcons, satisfies(Q-> 0 <= Q and Q <= 1)))->
> 

   if Q > 1/2 then -thisproc(1-Q)
> 

   elif Q=0 then -infinity
> 

   elif Digits <= evalhf(Digits) then evalhf(NormInv_inner(Q))
> 

   else NormInv_inner(evalf(Q))
> 

   fi
> 

:
> 

 

proc (X) options operator, arrow; 1/2+(1/2)*erf((1/2)*2^(1/2)*X) end proc

  

proc (x, Q) options operator, arrow; -(1/2)*2^(1/2)*Pi^(1/2)*erf((1/2)*2^(1/2)*x)*exp((1/2)*x^2)+2^(1/2)*Pi^(1/2)*Q*exp((1/2)*x^2)+x-(1/2)*2^(1/2)*Pi^(1/2)*exp((1/2)*x^2) end proc

 (1)
> 

f := proc(P)
  local n:
  for n from 1 to numelems(P) do NormInv(P[n]) end do:
end proc:
> 

with(Statistics):
> 

P := Sample(Uniform(0.5, 1), 10^3):
CodeTools:-Usage(Quantile~(Normal(0, 1), P)):
CodeTools:-Usage(Quantile~(Normal(0, 1), P, numeric)):
CodeTools:-Usage(f(P)):

memory used=17.03MiB, alloc change=0 bytes, cpu time=304.00ms, real time=304.00ms, gc time=0ns

memory used=17.14MiB, alloc change=0 bytes, cpu time=312.00ms, real time=312.00ms, gc time=0ns

Warning,  computation interrupted

  
> 

P := Statistics:-Sample(Uniform(0.5, 1), 10^1);

P := Vector[row](10, {(1) = .5814494501088662, (2) = .9192028848365151, (3) = .5837804651056354, (4) = .7511003077047363, (5) = .9996647395128597, (6) = .6777035756339407, (7) = .523538854774417, (8) = .6068302925603996, (9) = .6989195620513526, (10) = .6668340901901474}, datatype = float[8])

 (2)
> 

for n from 1 to 10 do
n;
NormInv(P[n])
end do;

1

  

.205602915707357559

  

2

  

1.39972985397027738

  

3

  

.211574416971719526

  

4

  

.677956326586950597

  

5

  

3.40135741509089939

  

6

  

Warning,  computation interrupted

  
> 

NormInv(0.969)

Warning,  computation interrupted

  
> 

 


 

Download NewtonFails.mw

@Axel Vogt RThank you Axel, it'...

mmcdara 7891
August 15 2019

@Axel Vogt 

Thank you Axel, it's indeed a very interesting paper.

 

@mehdi jafari  You cannot have a s...

mmcdara 7891
August 15 2019
0 0

@mehdi jafari 

You cannot have a set with elements 1, 1, 3, 3, 5, 6: in a set, each element is present only once.
So you can't do with sets what you did with lists.

@Carl Love @Moh HudaGiven the title...

mmcdara 7891
August 13 2019
0 0

@Carl Love 

@Moh Huda

Given the title of the plots I guessed D[Cchemo] is the Chemo Dose (which is of no practical help)

More interesting, I found this open source for the paper: https://arxiv.org/pdf/1806.05790.pdf.
D[chemo] is indeed the "chemotherapeutic dose" and it seems (but I'm not sure otf that) it' could be defined by relation (7).

Nevertheless I wouldn't go any further on my side.

Moh: your write is "My question is how to get the first two figures below based on these information...", my answer is: "do what is described in the paper cite, all the stuff seems to be here".
If you are puzzled with the term s∞, it's defined below equation (10): "The case with constant immunotherapy (s(t) = s∞)".
Maybe you should read the paper more carefully?

 

@Daniel Skoog  Thank you...

mmcdara 7891
August 12 2019
0 0

@Daniel Skoog 

Thank you

@Area51  Thanks for your reply. Th...

mmcdara 7891
August 09 2019
0 0

@Area51 

Thanks for your reply.
This would not be the first time that an article has contained typos. If you are familiar with the subject revelopped in the article (I'm not and I just looked at this article very quickly), maybe you could try to verify there is no typo in the differential equations (for instance eq 47 comes from "combining eqs 26 and 46: is it right).
Given the content of the article, and if you are familiar enough with Maple, you could try to use the Physics package.
In case you will face problems, there are people here that are specialists of it and who could help you.
 

Need corrections and precisions: The dif...

mmcdara 7891
August 09 2019
0 0

Need corrections and precisions:

  1. The differential equations in your introductory text differ from those in your worksheet:
    1. in the text it's written tau_0 and t0 in the worksheet (I suppose it's a typo in the later)
    2. in the text the 1st ode contains omega and the 2nd contains omega_m, while both of them contain omega in the worksheet.
      Text and worksheet disagree but, even if you code correctly the odes from the text, you will have a problem because the solutions you give in the text do not contain omega ???
       
  2. What are A, C_- and C_+ in the solution? Any constants?


In the worksheet above I replaced omega_m by omega in the solution you give and t0 by tau_0 in the each ode.
Next I checked if the solution you give was indeed a solution of the odes.
The answer is NO, unless A=0.


 

> 

restart
> 

sysode := [ 2*q*(3*q-1)*f1(tau)/tau^2+2*q*(diff(f1(tau), tau))/tau+diff(f1(tau), tau, tau)+(kappa^2+f2(tau))*(1+omega)*(tau/tau__0)^(-(3*(3+omega))*q) = 0, (54*q^3-30*q^2+4*q)*f1(tau)/tau^3+(24*q^2-4*q)*(diff(f1(tau), tau))/tau^2+11*q*(diff(f1(tau), tau, tau))/tau+diff(f1(tau), tau, tau, tau)-3*omega*(1+omega)*(kappa^2+f2(tau))*q*(tau/tau__0)^(-(3*(1+omega))*q)/tau = 0 ]:
> 

print~(sysode):

2*q*(3*q-1)*f1(tau)/tau^2+2*q*(diff(f1(tau), tau))/tau+diff(diff(f1(tau), tau), tau)+(kappa^2+f2(tau))*(1+omega)*(tau/tau__0)^(-3*(3+omega)*q) = 0

  

(54*q^3-30*q^2+4*q)*f1(tau)/tau^3+(24*q^2-4*q)*(diff(f1(tau), tau))/tau^2+11*q*(diff(diff(f1(tau), tau), tau))/tau+diff(diff(diff(f1(tau), tau), tau), tau)-3*omega*(1+omega)*(kappa^2+f2(tau))*q*(tau/tau__0)^(-3*(1+omega)*q)/tau = 0

 (1)
> 

indets~(sysode, function);

[{diff(diff(f1(tau), tau), tau), diff(f1(tau), tau), f1(tau), f2(tau)}, {diff(diff(diff(f1(tau), tau), tau), tau), diff(diff(f1(tau), tau), tau), diff(f1(tau), tau), f1(tau), f2(tau)}]

 (2)
> 

AssumedSolution := [
  f2(tau) = -A*(69*q^2-25*q+2)/(1+omega)*tau__0^(-9*q)*(tau/tau__0)^(3*omega*q)-kappa^2,
  f1(tau) = C__p*tau^mu__p+C__m*tau^mu__m + A*tau__0^(-9*q+2)
]:
print~(AssumedSolution):

f2(tau) = -A*(69*q^2-25*q+2)*tau__0^(-9*q)*(tau/tau__0)^(3*omega*q)/(1+omega)-kappa^2

  

f1(tau) = C__p*tau^mu__p+C__m*tau^mu__m+A*tau__0^(-9*q+2)

 (3)
> 

AssumedMu := [
  mu__m = (-2*q+1-sqrt(-20*q^2+4*q+1))/2,
  mu__p = (-2*q+1+sqrt(-20*q^2+4*q+1))/2
]:
print~(AssumedMu):

mu__m = -q+1/2-(1/2)*(-20*q^2+4*q+1)^(1/2)

  

mu__p = -q+1/2+(1/2)*(-20*q^2+4*q+1)^(1/2)

 (4)
> 

test := eval(AssumedSolution, AssumedMu);

[f2(tau) = -A*(69*q^2-25*q+2)*tau__0^(-9*q)*(tau/tau__0)^(3*omega*q)/(1+omega)-kappa^2, f1(tau) = C__p*tau^(-q+1/2+(1/2)*(-20*q^2+4*q+1)^(1/2))+C__m*tau^(-q+1/2-(1/2)*(-20*q^2+4*q+1)^(1/2))+A*tau__0^(-9*q+2)]

 (5)
> 

eval(sysode, test):
> 

simplify(%);

[A*(-69*tau__0^(-9*q)*(tau/tau__0)^(-9*q)*q^2*tau^2+25*tau__0^(-9*q)*(tau/tau__0)^(-9*q)*q*tau^2-2*tau__0^(-9*q)*(tau/tau__0)^(-9*q)*tau^2+6*tau__0^(-9*q+2)*q^2-2*tau__0^(-9*q+2)*q)/tau^2 = 0, A*q*(207*tau__0^(-9*q)*(tau/tau__0)^(-3*q)*omega*q^2*tau^2-75*omega*tau__0^(-9*q)*(tau/tau__0)^(-3*q)*q*tau^2+6*tau__0^(-9*q)*(tau/tau__0)^(-3*q)*omega*tau^2+54*tau__0^(-9*q+2)*q^2-30*tau__0^(-9*q+2)*q+4*tau__0^(-9*q+2))/tau^3 = 0]

 (6)
> 

 


 

Download question_2.mw

@Rouben Rostamian   Bravo!...

mmcdara 7891
August 07 2019
0 0

@Rouben Rostamian  

Bravo!

@Carl Love Don't you think that...

mmcdara 7891
August 07 2019
0 0

@Carl Love 

Don't you think that what glcrawfo is looking for is just something like this ?

 

> 

restart:
> 

f := (t, N) -> Sum((-1)^n*Heaviside(t-n), n=1..N)

proc (t, N) options operator, arrow; Sum((-1)^n*Heaviside(t-n), n = 1 .. N) end proc

 (1)
> 

R := 20 : C := 0.01 : L := 10:

> 

# choose N=2 for illustration

doe := (L*C*diff(x(t), t$2) + R*C*diff(x(t), t) + x(t)/C) = 200*f(t/3, 2);
value(doe)

.10*(diff(diff(x(t), t), t))+.20*(diff(x(t), t))+0.1e3*x(t) = -200*Heaviside(t-3)+200*Heaviside(t-6)

 (2)
> 

sol_exact := dsolve(value(doe), x(t));  #no need to fix the initial contitions here

x(t) = exp(-t)*sin(3*111^(1/2)*t)*_C2+exp(-t)*cos(3*111^(1/2)*t)*_C1-2*Heaviside(t-3)+(2/333)*111^(1/2)*Heaviside(t-3)*exp(-t+3)*sin(3*111^(1/2)*(t-3))+2*Heaviside(t-3)*exp(-t+3)*cos(3*111^(1/2)*(t-3))+2*Heaviside(t-6)-(2/333)*111^(1/2)*Heaviside(t-6)*exp(-t+6)*sin(3*111^(1/2)*(t-6))-2*Heaviside(t-6)*exp(-t+6)*cos(3*111^(1/2)*(t-6))

 (3)
> 

# numerical solution: here initial contitions are required


doe := (L*C*diff(x(t), t$2) + R*C*diff(x(t), t) + x(t)/C) = 200*f(t/3, 50):
sol_num := dsolve({value(doe),   x(0)=0, D(x)(0)=0  }, numeric, maxfun=10^6, method=rosenbrock, range=0..50);  
#                              it's just an example

`[Length of output exceeds limit of 1000000]`

 (4)
> 

plots:-odeplot(sol_num, [t, x(t)], t=0..50, refine=2)

 
> 

 


 

Download ode.mw

@Carl Love Thanks.I get no more err...

mmcdara 7891
August 06 2019

@Carl Love 

Thanks.
I get no more errors.

For information, here are the performances obtained

  • iMac  2.9 GHz Intel Core i5 (4 nodes, no hyperthreading), OS Mojave
  • time__run_Grid     =   8.266 s
  • time__run_noGrid = 31.515 s
  • speedup-factor     =   3.813      
  • efficiency              =    0.953 


In the attached file you will see that the two histograms present a clear visual discrepancy.
So I pushed your computations further in order to test if this discrepancy is "acceptable" or not. A ChiSquare GoF test returns a p.value of 5e-5 leading to this message :"This statistical test provides evidence that the null hypothesis is false" (here the null hypothesis is "Grid and noGrid are ""equivalent"").

Ok, the p.value is something to use carefully and one must not be too conclusive (more of this I didn't payed a lot of attention to the classes used fot thre ChiSquare test), but ... I will run again the Grid and noGrid simulations with larger samples and repetitions to clarify this former result


GridRandSample_test_1.mw




TESTS DONE
3 simulations with n3 = 10^4
No evidence against the null hypothesis

GridRandSample_test_2.mw
GridRandSample_test_3.mw
GridRandSample_test_4.mw

Very interesting ! I use to use the gris...

mmcdara 7891
August 06 2019

Very interesting !
I use to use the gris package to distribute Monte-Carlo simulations on multiple nodes.

I begin to examine your code with attention. Unfortunately for still a few weeks I can only use my personal Maple 2015.2 licence and I meet some errors when trying to execute it:

#Time the actual execution:

st:= time[real]():
   R1:= [Grid:-Seq['tasksize'= iquo(n3, num_nodes)](TRX1(k), k= [n2 $ n3])]:
time__run_Grid:= time[real]() - st
Error, (in TRX1) invalid input: add received 1 .. 100, which is not valid for its 2nd argument, i
Error, (in TRX1) invalid input: add received 1 .. 100, which is not valid for its 2nd argument, i
Error, (in TRX1) invalid input: add received 1 .. 100, which is not valid for its 2nd argument, i
Error, (in TRX1) invalid input: add received 1 .. 100, which is not valid for its 2nd argument, i

Maybe some version issue

I also obtain a speedup around 5.4 on 4 hyperthreaded proc. (Windows 7) ; in theory 8 proc should give a speedup close to 8 and hyperthreading lessens this value).
For information, I've got speedups larger than 6 with Windows XP: it seems that the intrinsic task balancing process of Windows 7 (and probably higher versions) has some "interference" with Grid (?).
You can have an idea of the task balancing process Windows >= 7 uses by looking the perfometer in the task manager: While all procs run when you run a calculation on a single proc, only proc 1 ran on Windows XP.
Maybe the smaller speedups obtained with Windows >= 7 come from the fact that the execution time on single proc run is lower than what it was on Windows XP (I didn't chek that).

I keep looking at your code.
 

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