Advanced Search

Bosco Emmanuel

10 Reputation

2 Badges

7 years, 74 days
Contact Bosco Emmanuel

MaplePrimes Activity


These are replies submitted by Bosco Emmanuel

This is what I meant...

Bosco Emmanuel 10
July 16 2018
0 0

@tomleslie 

Dear tom,

Congrats! You have done it. This is what I precisely meant when I said in an earlier mail that we should invoke the randomize command 3 times within the procedure. Apparently the random sampling is working in two steps as I anticipated in my previous mail.

Regarding the change of the sample with a change of parameter values, it is absolutely fine. That is the way it should be.

I am glad that no tinkering with my heart model is necessary.

In fact, just a few minutes ago, I was about to include the 3 randomize commands in the procedure and see. You have overtaken me!!! What a pleasant surprise !!!

Thanks again for all the pains you have taken for the benefit of my research.

Cordially,

Bosco Emmanuel

Exponential Sampling...

Bosco Emmanuel 10
July 14 2018
0 0

@tomleslie 

Dear Tom,

When I invoke the exponential Sampling, say I want to generate 1000 numbers distributed as per the exponential, what in my opinion happens is:

1. The random number generator starts with a Seed and generates 1000 numbers distributed uniformly between 0 and 1. At this stage the parameters like p1etc. will not enter.

2. In step 2, the generator will take each of these 1000 numbers and apply the inverse exponential function to generate 1000 other numbers. And this is the random sequence used within the procedure. At this second stage only, parameters like p1 enters.

Hence I would think that one can independently vary the Seed and the parameter values. 

Please let me know if you would agree with this perception.

Thanks.

Bosco Emmanuel

Randomize and Freeze the Seed ?...

Bosco Emmanuel 10
July 13 2018
0 0

@tomleslie 

Dear Tom,

Do you think invoking randomise(0) three times inside the procedure will freeze the Seed and thereby ensure reproducibility of the results ? 

Thanks.

Bosco Emmanuel 

Seed...

Bosco Emmanuel 10
July 13 2018
0 0

@tomleslie 

Dear Tom,

I would think 26 evaluations of lc would be needed to find the direction of steepest descent. 

It may not be possible to avoid random sequences in the model. However if I can freeze the Seed while generating the random sequence, the sampling will generate the same sequence for the same value of p. Will this help ?

I am also exploring alternative measures in the place of the Euclidean Distance which is presently used. Correlation is one such. 

Thanks. 

Bosco Emmanuel 

Maple 2018...

Bosco Emmanuel 10
July 12 2018
0 0

@tomleslie 

Dear Tom,

Maple 2018 is now installed in my system. I ran your re-do with the Minimize command. It was still evaluating after 20 minutes when I stopped it. DirectSearch also is included in maple 12 and I have tried some examples which are working. Today I hope it will be included in maple 2018 as well. Hence from now on we will be able to compare results without any version conflict. [ I understand from maple application notes that DirectSearch works in maple12 and above]

I have noted that you have really done some relevant and interesting in-depth analysis. Though I understand your reasoning I do have some comments and questions.

Firstly why should it call "lc" 117 times for one iteration? I thought "lc" will be called by the Minimize only when the values of the parameters are changed. Do you mean to say that "lc" will be called multiple times for the same set of values of the parameters? I do not quite understand why it should do this.

Then, when you find "lc" fout times for the same set of parameter values, it produces four different Error values. I understand this as you youself has understood as due to the differences in the random set. However we should note that the first 3 values are nearly the same while the fourth is somewhat off. An analogy may be useful here: Suppose that I toss a coin 1000 times. Let p be the probability of getting heads and (1-p) the prbabilty of getting tails. Of course each time I perform this experiment I will obtain a different sequence of heads and tails. When I find the fraction of times I get heads in each sequence, it will be close to p for all the experiments. p is like the statistical parameters which we are trying to find by minimizing the Error. What we are after is "statistical parameter estimation" by comparing real data ( in our case, an ECG data [of a real patient] consisting of 1000 heart beats) with a corresponding synthetic ECG data which depends on 6 physiological parameters. Out of these 6 parameters, the first 3 characterize 3 neuronal spike sequences and the last three denote 3 physiological time parameters. These 6 parameters will in general vary from patient to patient and vary in the same patient if his health/stress/emotional status changes.

All said and done, your analysis may point to the possibility that the "Error" may not be the correct or robust measure of comparison between the real and synthetic data. Perhaps I should think of alternative measures of comparison. However the same code will be executed with the new measure substituting for the "Error". The numbers that you have generated with DirectSearch is perplexing! I was reading in the DirectSearch pdf that it can find multiple local minima as well the global minimum. What you are seeing is not even a set of local minima. Instead you get a different result each time you run it. Hence I tend to agree with you that the differences arise out of the different random set used each time. If this is the case, we need to find and use a more robust measure in the place of "Error'. 

One mystery remains: the result of "Minimize" looks more robust than the result of "DirectSearch".

Thanks again for all your good efforts.

What I plan to do now is to reproduce all your results with Maple 2018 and DirectSearch and explore with the same old measure of comparison for some more time before trying alternative measures.

I will continue to keep in touch with you, Tom.

Cordially,

Bosco Emmanuel

 

Your Re-Do Works...

Bosco Emmanuel 10
July 11 2018
0 0

@tomleslie 

Dear Tom,

Please find enclosed two segments of the code: one with the file read inside the procedure and the other outside. Both works fine. However the iteration does not start and it reports the initial point as the result. Is the procedure caught in a local minimum ? I tried varying the initial point and I could not get the iteration started. I also tried using values close to the final values that you obtained using DirectSearch. Though the error is much less this time, the iteration fails to start off. As you obtained final values much different from the initial point in DirectSearch, I guess that DirectSearch would have done several iterations.

Kindly advise. Thanks again.

Cordially,

Bosco Emmanuel
 

> 

restart; with(Statistics); with(Optimization); lc := proc (p1, p2, p3, p4, p5, p6) local RT, TC, N, T, z, n, l, B, p; p := [p1, p2, p3, p4, p5, p6]; RT := Array(1 .. 3, p[4 .. 6]); TC := Array(0 .. 1001, fill = 0); N := Array(1 .. 3, [0, 0, 0]); T := Matrix(3, 1001, [Sample(Exponential(1/p[1]), 1001), Sample(Exponential(1/p[2]), 1001), Sample(Exponential(1/p[3]), 1001)], scan = columns); for n to 1001 do if T[1, n] <= T[2, n] and T[1, n] <= T[3, n] then l := 1 elif T[2, n] <= T[1, n] and T[2, n] <= T[3, n] then l := 2 else l := 3 end if; N[l] := N[l]+1; B[N[l], l] := T[l, n]+RT[l]; TC[n] := TC[n-1]+B[N[l], l] end do; z := ImportMatrix("/home/bosco/Data-1", source = delimited, delimiter = "\t", format = rectangular, datatype = sfloat, transpose = false, skiplines = 0); return add((z[i, 1]-TC[i+1]+TC[i])^2, i = 1 .. 1000) end proc; Minimize(lc, initialpoint = [3, 1, .5, .3, 1, 1.5], assume = nonnegative)

Warning, no iterations performed as initial point satisfies first-order conditions

  

[270.493012299041027, Vector(6, {(1) = 3., (2) = 1., (3) = .500000000000000000, (4) = .299999999999999989, (5) = 1., (6) = 1.50000000000000000})]

 (1)
> 

restart; with(Statistics); with(Optimization); z := ImportMatrix("/home/bosco/Data-1", source = delimited, delimiter = "\t", format = rectangular, datatype = sfloat, transpose = false, skiplines = 0); lc := proc (p1, p2, p3, p4, p5, p6) local RT, TC, N, T, n, l, B, p; p := [p1, p2, p3, p4, p5, p6]; RT := Array(1 .. 3, p[4 .. 6]); TC := Array(0 .. 1001, fill = 0); N := Array(1 .. 3, [0, 0, 0]); T := Matrix(3, 1001, [Sample(Exponential(1/p[1]), 1001), Sample(Exponential(1/p[2]), 1001), Sample(Exponential(1/p[3]), 1001)], scan = columns); for n to 1001 do if T[1, n] <= T[2, n] and T[1, n] <= T[3, n] then l := 1 elif T[2, n] <= T[1, n] and T[2, n] <= T[3, n] then l := 2 else l := 3 end if; N[l] := N[l]+1; B[N[l], l] := T[l, n]+RT[l]; TC[n] := TC[n-1]+B[N[l], l] end do; return add((z[i, 1]-TC[i+1]+TC[i])^2, i = 1 .. 1000) end proc; Minimize(lc, initialpoint = [5, .1, .5, .4, .97, .96], assume = nonnegative)

Warning, no iterations performed as initial point satisfies first-order conditions

  

[69.2742545614697036, Vector(6, {(1) = 5.0, (2) = .1, (3) = .5, (4) = .4, (5) = .97, (6) = .96}, datatype = float[8])]

 (2)
> 

NULL


 

Download Tom-11-July.mw

It is good...

Bosco Emmanuel 10
July 10 2018
0 0

@tomleslie 

Dear Tom,

I just returned home from work. Your suggestions will be very helpful I believe.

My goal is to find the values of the six parameters that minimise the Error.

Tomorrow I will run the code which you sent today using the Minimise function after duly following your advice regarding the comments.

You are very correct about the wasteful data file reading multiple times.I will try to move it out of the procedure. However I will first run your code without this change. I am prepared to wait for longer than ten minutes.

I will also explore the DirectSearch once I get the basic code running.

The results you have already got are very encouraging.

Thanks again for all your help. 

I will report to you tomorrow.

Cordially,

Bosco Emmanuel 

OPTIMIZATION-TRIAL-1...

Bosco Emmanuel 10
July 10 2018
0 0

@tomleslie 

Dear Tom,

I executed in Maple 12 the code which you sent last. In my system it needed some changes to run properly. In particular, I need to remove the "local" declaration.

As you will see in the code I am uploading, the procedure is calculating the Error though it rings a warning regarding the absence of "local" declaration.

The fatal error however is "non-numeric result encountered" in the optimization.

I would appreciate your kind advice.

Thanks.

Bosco Emmanuel
 

> 

restart; with(Statistics); with(Optimization); Minimize(lc([x1, x2, x3, x4, x5, x6]), initialpoint = [x1 = 3, x2 = 1, x3 = .5, x4 = .3, x5 = 1, x6 = 1.5]); lc := proc (p::list) RT := Array(1 .. 3, p[4 .. 6]); TC := Array(0 .. 1001, fill = 0); N := Array(1 .. 3, [0, 0, 0]); T := Matrix(3, 1001, [Sample(Exponential(1/p[1]), 1001), Sample(Exponential(1/p[2]), 1001), Sample(Exponential(1/p[3]), 1001)], scan = columns); z := ImportMatrix("/home/bosco/Data-1", source = delimited, delimiter = "\t", format = rectangular, datatype = sfloat, transpose = false, skiplines = 0); for n to 1001 do if T[1, n] <= T[2, n] and T[1, n] <= T[3, n] then l := 1 elif T[2, n] <= T[1, n] and T[2, n] <= T[3, n] then l := 2 else l := 3 end if; N[l] := N[l]+1; B[N[l], l] := T[l, n]+RT[l]; TC[n] := TC[n-1]+B[N[l], l] end do; return add((z[i, 1]-TC[i+1]+TC[i])^2, i = 1 .. 1000) end proc; Error := lc([3, 1, .5, .3, 1, 1.5])

Error, (in Optimization:-NLPSolve) non-numeric result encountered

  

Warning, `RT` is implicitly declared local to procedure `lc`

  

Warning, `TC` is implicitly declared local to procedure `lc`

  

Warning, `N` is implicitly declared local to procedure `lc`

  

Warning, `T` is implicitly declared local to procedure `lc`

  

Warning, `z` is implicitly declared local to procedure `lc`

  

Warning, `n` is implicitly declared local to procedure `lc`

  

Warning, `l` is implicitly declared local to procedure `lc`

  

Warning, `B` is implicitly declared local to procedure `lc`

  

270.4930125

 (1)
> 

``


 

Download OPTIMIZATION-TRIAL-1.mw

randStuff4.mws...

Bosco Emmanuel 10
July 09 2018
0 0

@tomleslie 

Dear Tom,

I saw you mail. I am back in my office today.

There is some difficulty in opening randStuff4.mws that you sent. Could you please insert it using the green arrow ? so that I can copy, paste and run it. Thanks.

Cordially,

Bosco Emmanuel

 

Glad to know it was not my fault + Some ...

Bosco Emmanuel 10
July 06 2018
0 0

@tomleslie 

Dear Tom,

I was glad to know that it was not my fault.

Your observations about the "." s are interesting and will be useful to me. I thoughly understand what you are saying.

Then the "Quirk" also is interesting. Will you attribute the difference to "Seed Change" in one case and "No Seed Change" in the other ?

My present difficulty is this:

When I tried to run the codes which you sent most recently, they refuse to run properly and throw up errors. This time I will study it myself to the extent possible with my limited knowledge of maple. If I give up finally, I will approach you again. I have already taken much of your time, I am afraid.

In addition to this mystery, I have one more mystery to solve : With my clumsy code, the top-level works fine and gives me the correct result for the Error. It is mysterious that my clumsy code does not work properly inside a procedure.

I am left wondering whether these two mysteries have a common cause.

My profound thanks to you once again.

Cordially,

Bosco Emmanuel

Thanks. Regrets for the troubles taken f...

Bosco Emmanuel 10
July 04 2018
0 0

@tomleslie 

Dear Tom,

I didn't retype anything. I just copied and pasted. Perhaps I failed to do this correctly.

Now I am glad that all the four versions gives the same and correct result. I am extremely grateful for your help and sorry about the pains taken by you.

The task that remains for me is to understand your very elegant redo of my clumsy code. 

I will go to my office day after tomorrow and run all the four codes send by you.

Thanks once again.

Cordially,

Bosco Emmanuel 

 

Your Re-do-s and my data to use...

Bosco Emmanuel 10
July 04 2018
0 0

@tomleslie 

Dear Tom,

I tried your re-do-s with my data. They are not yet running properly. I am enclosing  my data as Tom.txt and your re-do-s executed with my data. Please have a look.

Thanks again.

Cordially,

Bosco Emmanuel
 

> 

" restart; with(Statistics): # #` Initialise a few things` # p:= [3, 1, 0.5, 0.3, 1, 1.5]: RT:= Array(1..3, p[4..6]): TC[0]:= 0.: N:= Array(1..3, [0, 0, 0]): T:= Matrix( 3, 1001, [ Sample(Exponential(1/p[1]), 1001), Sample(Exponential(1/p[2]), 1001), Sample(Exponential(1/p[3]), 1001) ], scan=columns ): # #` Do the calculation` # for n to 1001 do l:= min[index](T[..,n]): N[l]:= N[l]+1; B[N[l], l]:= T[1,n]+RT[l]; TC[n]:= TC[n-1]+B[N[l], l] end do: z:= ImportMatrix("\" ""/home/bosco/Data-1\"", source=delimited, delimiter="\"\\t"" \""): # #` Get the result` # Error:= add((z[i,1]-(TC[i+1]-TC[i]))^2, i = 1 .. 1000); "

Error, Got internal error in Typesetting:-Parse:-Postprocess : "internal error: invalid object " "

" restart; with(Statistics): # #` Initialise a few things` # p:= [3, 1, 0.5, 0.3, 1, 1.5]: RT:= Array(1..3, p[4..6]): TC[0]:= 0.: N:= Array(1..3, [0, 0, 0]): T:= Matrix( 3, 1001, [ Sample(Exponential(1/p[1]), 1001), Sample(Exponential(1/p[2]), 1001), Sample(Exponential(1/p[3]), 1001) ], scan=columns ): # #` Do the calculation` # for n to 1001 do l:= min[index](T[..,n]): N[l]:= N[l]+1; B[N[l], l]:= T[1,n]+RT[l]; TC[n]:= TC[n-1]+B[N[l], l] end do: z:= ImportMatrix("\" ""/home/bosco/Data-1\"", source=delimited, delimiter="\"\\t"" \""): # #` Get the result` # Error:= add((z[i,1]-(TC[i+1]-TC[i]))^2, i = 1 .. 1000); "

  
> 

``


 

Download Tom-Redo-1.mw
 

> 

" restart; with(Statistics): # #` Initialise a few things` # p:= [3, 1, 0.5, 0.3, 1, 1.5]: RT:= Array(1..3, p[4..6]): TC[0]:= 0.: N:= Array(1..3, [0, 0, 0]): T:= Matrix( 3, 1001, [ Sample(Exponential(1/p[1]), 1001), Sample(Exponential(1/p[2]), 1001), Sample(Exponential(1/p[3]), 1001) ], scan=columns ): # #` Do the calculation` # for n to 1001 do l:= min[index](T[..,n]): N[l]:= N[l]+1; B[N[l], l]:= T[1,n]+RT[l]; TC[n]:= TC[n-1]+B[N[l], l] end do: z:= ImportMatrix("\" ""/home/bosco/Data-1\"", source=delimited, delimiter="\"\\t"" \""): # #` Get the result` # Error:= add((z[i,1]-(TC[i+1]-TC[i]))^2, i = 1 .. 1000); "

Error, Got internal error in Typesetting:-Parse:-Postprocess : "internal error: invalid object " "

" restart; with(Statistics): # #` Initialise a few things` # p:= [3, 1, 0.5, 0.3, 1, 1.5]: RT:= Array(1..3, p[4..6]): TC[0]:= 0.: N:= Array(1..3, [0, 0, 0]): T:= Matrix( 3, 1001, [ Sample(Exponential(1/p[1]), 1001), Sample(Exponential(1/p[2]), 1001), Sample(Exponential(1/p[3]), 1001) ], scan=columns ): # #` Do the calculation` # for n to 1001 do l:= min[index](T[..,n]): N[l]:= N[l]+1; B[N[l], l]:= T[1,n]+RT[l]; TC[n]:= TC[n-1]+B[N[l], l] end do: z:= ImportMatrix("\" ""/home/bosco/Data-1\"", source=delimited, delimiter="\"\\t"" \""): # #` Get the result` # Error:= add((z[i,1]-(TC[i+1]-TC[i]))^2, i = 1 .. 1000); "

  
> 

``


 

Download Tom-Redo-1.mw
 

> 

restart; with(Statistics); lc := proc (p::list) local RT, TC, N, T, z, n, l, B; RT := Array(1 .. 3, p[4 .. 6]); TC := Array(0 .. .1001, fill = 0); N := Array(1 .. 3, [0, 0, 0]); T := Matrix(3, 1001, [Statistics:-Sample(Exponential(1/p[1]), 1001), Statistics:-Sample(Exponential(1/p[2]), 1001), Statistics:-Sample(Exponential(1/p[3]), 1001)], scan = columns); z := ImportMatrix("/home/bosco/Data-1", source = delimited, delimiter = "\t"); for n to 1001 do l := min[index](T[() .. (), n]); N[l] := N[l]+1; B[N[l], l] := T[1, n]+RT[l]; TC[n] := TC[n-1]+B[N[l], l] end do; return add((z[i, 1]-TC[i+1]+TC[i])^2, i = 1 .. 1000) end proc; Error := lc([3, 1, .5, .3, 1, 1.5])

Error, unable to parse

"#` Define procedure which does the same` #` calculation as the above` # restart: with(Statistics): lc:=proc( p::list ) local RT:= Array(1..3, p[4..6]), TC:= Array(0...1001, fill=0), N:= Array(1..3, [0, 0, 0]), T:= Matrix( 3, 1001, [ Sample(Exponential(1/p[1]), 1001), Sample(Exponential(1/p[2]), 1001), Sample(Exponential(1/p[3]), 1001) ], scan=columns ), z:= ImportMatrix("/home/bosco/Data-1", source=delimited, delimiter="\\t"), n, l, B; uses Statistics: for n to 1001 do l:= min[index](T[..,n]): N[l]:= N[l]+1; B[N[l], l]:= T[1, n]+RT[l]; TC[n]:= TC[n-1]+B[N[l], l] end do: return add((z[i,1]-(TC[i+1]-TC[i]))^2, i = 1 .. 1000); end proc: # #` Run the calculation with the same data` #` as used in the above` # Error:=lc([3, 1, 0.5, 0.3, 1, 1.5]); "

  
> 

``


 

Download Tom-Redo-2.mw

 

Download Tom.txt

 

 

 

 

Very Elegant. The Logic the same ?...

Bosco Emmanuel 10
July 02 2018
0 0

@tomleslie 

Dear Tom,

Your redo is very elegant. To satisfy myself that it does the same computation as my code, the first thing I do tomorrow when I go to my office is to use my own data file and see if it reports the same result as mine.

Thanks for all your kind efforts.

Cordially,

Bosco Emmanuel 

@tomleslie  Dear Tom, I am sending...

Bosco Emmanuel 10
July 02 2018
0 0

@tomleslie 
Dear Tom,

I am sending you two maple codes: Code-1 which computes the Error without invoking a procedure and Code-2 computes the same Error but within a procedure. Code-2 fails to report the value of the Error while Code-1 correctly reports it. Please advise me what went wrong with Code-2.

My final goal is to use this procedure in maple optimisation to find the values of the six parameters p1 to p6 which minimise the Error.

My sincere thanks in advance.

Cordially,

Bosco Emmanuel

 

> 

restart; with(Statistics); p1 := 3; p2 := 1; p3 := .5; p4 := .3; p5 := 1; p6 := 1.5

RT[1] := p4; RT[2] := p5; RT[3] := p6; TC[0] := 0.; N[1] := 0; N[2] := 0; N[3] := 0; X1 := Exponential(1/p1); X2 := Exponential(1/p2); X3 := Exponential(1/p3); for n to 1001 do T[1] := Sample(X1, 1); T[1][1]; T[2] := Sample(X2, 1); T[2]; T[3] := Sample(X3, 1); T[3]; TMIN := min(T[1], T[2], T[3]); T1 := T[1][1]; T2 := T[2][1]; T3 := T[3][1]; if T1 <= T2 and T1 <= T3 then l := 1 else end if; if T2 <= T1 and T2 <= T3 then l := 2 else end if; if T3 <= T2 and T3 <= T1 then l := 3 else end if; l, T[l][1]; N[l] := N[l]+1; B[N[l], l] := T[l][1]+RT[l]; TCP[N[l], l] := TC[n-1]+B[N[l], l]; TC[n] := TC[n-1]+B[N[l], l] end do

for i to 1001 do R[i] := TC[i] end do; RR := Vector(1 .. 1000); for i to 1000 do RR[i] := R[i+1]-R[i] end do

z := ImportMatrix("/home/bosco/Data-1", source = delimited, delimiter = "\t", format = rectangular, datatype = sfloat, transpose = false, skiplines = 0); x := Array(1 .. 1000); for i to 1000 do x[i] := z[i, 1] end do; Error := add((x[i]-RR[i])^2, i = 1 .. 1000)
``

254.9174755

  

254.9174755

 (1)
> 

``


 

> 

"restart: with(Statistics): lc:=proc(p1,p2,p3,p4,p5,p6) "

RT[1] := p4:

for i to 1001 do R[i] := TC[i] end do:

" z:=ImportMatrix("/home/bosco/Data-1", source=delimited, delimiter="\\t", format=rectangular, datatype=sfloat, transpose=false, skiplines=0): x:=Array(1..1000):for i from 1 to 1000 do x[i]:=z[i,1]: od : Error:=add((x[i]-RR[i])^2,i=1..1000); Return [Error]: end proc: lc(3,1,.5,.3,1,1.5); "
``
> 

p1 := 3; p2 := 1; p3 := .5; p4 := .3; p5 := 1; p6 := 1.5


 

Download Code-2-for-Tom.mw

Download Code-1-for-Tom.mw

Thanks. I have a problem...

Bosco Emmanuel 10
June 30 2018
0 0

@tomleslie 

Dear Tom,

Thanks for your kind reply. I have a problem with a maple code that invokes random numbers inside a procedure. I will send you the code on Monday for your kind perusal and help.

Thanks.

Cordially,

Bosco Emmanuel

1 2 3 Page 1 of 3