abcd

45 Reputation

6 Badges

6 years, 294 days

MaplePrimes Activity


These are replies submitted by abcd

@nMaple 

 

I'm not sure I adress exactly your problem but I propose you this piece of code (which needs probably a lot of improvements).
It enables to make a contourplot map over the "ternary" triangle.
I give four illustrative examples : 3 of them are maps of functions that depend only a single component (a, b or c) ; the fourth one is a more complex function.

I hope that will help you,
let me know

 

 

Download TernaryContourPlot.mw

@nMaple 

 

Where your points come from ?

I mean, do you really have points, or those points are, for instance, particular solution of kinetic equations ? In sucg a case trying to plot might net be the correct way to proceed ?

 

 

Ok with Markyian : check your formulation.

Next, with freeze and thaw you will be able to use something like that

 

Download prop.mw

@nMaple 

 

I can't load your image but I think I can figure out what you want

Suppose you have a point P correpondind to the mixture [a, b, c=1-b-a]
Let us denote (this choice is arbitrary) by A the blue side of the whole triangle and by B the green one


P lies at the intersection of the red diagonal line starting from "abscissa a on A and of the blue line starting from "abscissa" b on B

I assume the triangle is represented in the orthogonal [x, O, y] frame with O confounded with the bottom-left corner of the triangle. I assume more the sides of the triangle have same lengths 1
The starting point from A is thus 1-a = 0.7 and the red line to follow from this point has equation

y[red] = tan(2*Pi/3)*(x-(1-a))

The starting point from B has coordinates (b*cos(Pi/3), b*sin(Pi/3)) and so the horizontal blue line has equation

y[blue] = b*sin(Pi/3))

Equating y[red] and y[blue] will give you the abscissa x[P] of P, and text is y value y[P]


*sin(Pi/3) above side A
The intersection point 

@Shervin 


The first and simplest method you can use is "brute force" :

  1. Draw a sample (a1, ..., aN) of size N from the lognormal distribution for A'
  2. Solve iteratively the LP system (LPsolve command) afer having instanciated A' by each a1, ..., aN

If N is sufficiently large you will be able to compute a lot of statistics about x.
To help you in the choice of N (10, 100, ... 10^6 ...?) you can use the Bootstrap procedure from package Statistics.
A basic idea is :

  1. Choose a reasonable value for N (it depends on the size of your system and thus on the computational time).
    To illustrate let us say N = 1000
  2. When your N resolutions are done use the Bootstrap strategy to assess the accuracy of your results (which can be the mean of x, or the probability that x exceeds some value ...)
  3. If you consider this accurracy too wide, increase N and go to point 1.

A tip : if the accuracy for N simulations is acc(N), then increasing N by a factor K fill result in an accuracy acc(k*N) roughly eqaul to acc(N)/sqrt(K)



Just a remark : a lognormal law may be sympathetic (not too far frtom a normal law) or not, when it have a strong mass close to 0.
In the later case, A'x=B could introduce numerical problems for x=B/A'.

When sigma is small, some normal approximation could be usefull (this is just an idea commonly used in other domains of Statistics). Maybe you could compute a single solution for the mean of A' (not too far from its mode), and then use a "progation of variance approach"
(see for example https://en.wikipedia.org/wiki/Propagation_of_uncertainty)
But the interest of this approach realy depends on your objectives


I don't think I could provide you more help.

Good luck





 

@Carl Love 

 

What the OP wants is "How to random generate 2 integer number with condition x must be less than y ?"

What are x and y ?

Is there any natural interpretation of this statement other than sayiong "y is the first number generated" and s is the second one ?
Obviously x cannot be the first sampled number, for it would be stupid to say "x must be less than y" (a correct statement would be here "y is larger or equal to x")
If x and y are sampled "simultanously", there is no reason at all to say : I want "x must be less than y". Let the chance do the job, and if it happens that x is less than you will be lucky.

So, there is no doubt for me (even if, I concede you that the expression of the problem should have been less ambigious), and I believe I am the the only one to interpret the OP's statement that way, that his very idea is "I would like to sample 2 random numbers between 1 and 7 in such a way that the second is less than the first one".



Let us look now to the (I conceed this, very clever) procedure GenXY:= ()-> combinat:-randcomb(7, 2)[]; 

What does it do exactly ?

Suppose you have an urn with 7 balls labelled 1 to 7 :  GenXY() pick 2 balls out of this urns and says "the smallest number is ... and the highest one is ...".
Which is not at all the same thing than sampling one ball and trying to sample a second one with a smaller number written on it. 
The fact that combinat:-randcomb(7, 2) returns a combination (for example (2,3) and (3,2) are confounded) means there is no order in the sampling : it is impossible to know which one of x or y had been the first picked out from the urn.

I consider (it is my own opinion; and the OP should intervene in this "philosophical debate")  that  combinat:-randcomb(7, 2) does not fit the requirement of the OP.

 

 

@marc sancandi 

 

I send this message to Carl too


Let X a discrete uniform random variable over [1, N]

A bayesian approach :

We define the random variable Y that way :

for any y > 0 : Prob(Y=y|X=x) = 1/(x-1) if x > 1 and 0 otherwise

For the computations below it is more convenient to modify this
definition by discarding the particular case x=1.
To do this, consider now X a discrete uniform random variable over [2, N]

What is the probability P(Y=y) that Y takes the value y ?

Bayes' law says :

        Prob(Y=y) = Sum(P(Y=y|X=x)*P(x), x=2..N)
 where  Prob(X=x) = 1/(N-1) for x=2,..,N


Let N = 7 ; what is the value of P(Y=1) ?
One finds :
    P(Y=1) = P(Y=1|X=2)P(2) + ... + P(Y=1|X=7)P(7)

    P(Y=1) =    1*(1/6)     + ... +  (1/6)*(1/2)
    
    P(y=1) = 49/120  (about 0.408)

Doing the same for y=2, ..., 6 gives

P(Y=2) = 29/120
P(Y=3) = 19/120
P(Y=4) = 37/360
P(Y=5) = 11/180
P(Y=6) = 1/36

One easily verifies that P(Y=1) + ... + P(Y=6) = 1

 
What does this mean ?

Suppose you want do draw the mass function (the discrete analogous of
the ptobability density function for continuous random variables) with
"buildings" oh height H(x,y) : What should it look like?

  
First  thing : the sum of all the heights H(x,y) must be 1.
Second thing : obviously H(x,y) = 0 for each y greater or equal to x

Consider now the couple (X=2, Y=1) :
H(x,y) must represent the probability P(Y=1|X=2)... which is 1*(1/6) !

Consider now the couple (X=3, Y=1) :
H(x,y) must represent the probability P(Y=1|X=3)... which is 1/2*(1/6),
because the two events P(Y=1|X=3) and P(Y=2|X=3) have exactly the same
probability.

And so on.

Let us count the total number C of configurations.
We have C = sum((n-1), n=2..N) which, for N=7 gives C=21
It is clear that these 21 buildings do not have the same heights.
The highest has value 1/6, the smallest 1/216.
Keep this value of 1/6 in mind, I will go back on it later.


More precisely, all the buildings H(x,y[j]) for j from 1 to x-1
have the same height, equal to 1/(x-1).

This has an important implication :
      sum(H(x,y[j]), j=1..x-1) = 1/6 for each x from 2 to 7
whch is meant by saying that "the marginal distribution of Y" is
discrete uniform.


Now, what about this value 1/6 I talked before ?
Suppose you realise a simulation to build the mass function of
the couple (X, Y) and you use 2^17 samples.
The the height of the highest (H(2,1) building is 2^17/6, a value
close to 21845.
This is practically the value the second histogram (generated
by GenXY2) gives.

In fact, besides any intuition, the second histogram is the
good one, and the first is wrong.


Maybe you wonder where is the error in GenXY ?
It comes from the instruction  `if`(x<y, [x,y], [y,x])[]

let me introduce a new random variable denoted U and defined by
for any u > 0 : Prob(U=u|X=x) = 1/x if x > 0
To picture that, remember that (X, Y) lies in the "triangle"
strictly below the first diagonal : U is just the complement
to Y and thus (X, U) lives in the triangle above, diagonal
included.

If one consider that the 'else' clause of the previous test returns
the a sample of (X, Y), than the 'then' clause returns a sample of
(X,U).
But, because Y is in [1, x-1] and U in [x, N], the test just returns
a sample of a couple of iid random variables of discrete uniform
distribution over the {1, .., N}x{1, .., N}.

So, all tge "buildings" must have here the same heights.

Carl's idea to return [x,y] or [y,x] to avoid discarding
"in mean" 1 our of 2 trials was a priori a good trick ...
But it returns a wrong result.
The good test is  : if y<x then [x, y] end if

Which is exactly what John did.

So, John, do not worry, you are right


Download a2.mw



@John Fredsted 

 

Why do you say tou are wrong ?
Are you sure the second histogram is false ?

@Carl Love 

 

A good example for doing probabilities ...

What you think is the correct histogram ?

By the way,
I use cmaple too and one one of its limitations is concerns plots. Even if you use plotsetup and generate a jpeg file, it has really a very poor quality (no colors, bad labels or titles and so on).
Did you faced this problem ?
If yes, how do you manage it ?

 

Thanks

I have the same problem to access from the outside (Eclipse/Javap hlep pages I made. Guess it is impossible (?)

It works pretty well for me (Maple 2015, Mac OS X El Capitan) 

I can't say tou much ... sorry

Download ok.mw

@Preben Alsholm 

 

Thant you Preben, I learnt a lot of tricks from your answer (use of selectremove and getdata for instance)
A good point is that the two dofferent evaluations of the solution in the selected ranhe give exactly the same values (which was not exactly the case for my own computations).

Besides that, you have probably observed what happens when you use lsode[adamsfull] : these same evaluations give completely different results.
The understanding of these discrepancies is one of my major concerns and  I don't  have the feeling that I have had a satisfactory answer.


I also submited the question to the Maplesoft support team (France) ; for the moment I only recieved a I received a cryptic reply suggesting that the external EDO solver (I guess Maplesoft uses a well established scientific library) could have "a bug" (a problem in the exchange of informations between itself and Maple seems more likely to me).

I'll keep in touch with you to give you the end of the story.

 

Thanks again 

Marc

 

Hi again, 

Just a little Maple code to complete my previous answer.

I construct a Distribution the pdf of which (named MyPDF) is made of  two gaussian distributions centered at -3 and +3, with std dev of 1.
Obviously it possesses 2 modes (-3 and +3).
The pdf being formal, does Maple will find the two modes ?

  • Trial 1 :
    define a random variable X by X := RandomVariable(MyPDF)
    then compute Mode(X)
    The answer says the  modes are solution of a combination of LambertW functions (the warning Maple returns too is very instructive)
    Applying "allvalues" to the answer returns modes at position -3 and 9

  • Trial 2 : do a sample of X : Sample(X, 10000)
    Then compute tje mode numerivally : Mode(X, numeric)
    Maple gives the two modes -3 and +3 (up to some uncertainty due to de limited number of samples) ... but returns also a mode at 0 (which obviously is a minimum !!!)
    All happens as if Maple was searching extrema of a function and not maxima (?)


As a matter of conclusion : be carefull when you do statistics with Maple. Iit is not a bad tool to do this, just a tool which is not originally designed for. So keep a critical look to the results, especially if you have in mind to do statistical data analysis.
Maybe the 2016 release does a step forward in this direction, nevertheless Maple won't never be a reference tool for statisticians.

To end with a good point : the possibility to do some symbolic calculus with random variable is a very intresting feature

 



Download MOI---Bimodal--2.mw

Hi, 

You are perfectly right, your distribution is bimodal an 22 and 23 are the modes.

That is theory.

Now think to this bimodal distribution of random variable Z 

  • X is normal of mean -3 and std deviation 1
    Y is normal of mean +3 and std deviation 1
    Z is defined by : pick X with probability 1/3 and Y with probability 2/3
  • Sampling X and Y and plotting the histogram for Z will exhibit a camel curve with a bum around -3 and another twice height around +3
  • One uses to say that Z is bimodal and that -3 and +3 are two different modes.

This raises a first question : when a software (not specially devoted to handle Statistics, which is the case of Maple), what value(s) will it return ?
Will it output +3, or +3 and -3 ?


Let us consider nowthis second example :

  • X is a sample from a gaussian variable of mean -3 and std deviation 1
    Y is -X
    W is the sample obtained by concatenation of X and Y
  • Obviously W are two modes of exactly equal height at symetric positions (saying -3 and +3)





When I am concerned with statistical problems, I prefer to use a true statistical softwares, R for instance.
R provides two type of answers

  1. TYPE.1 if there is a unique maximum maximorum (as it is for Z above), R returns it (here +3)
  2. TYPE.2 if there two modes of equal height (W serie) returns only one of them ... but warns you that the distribution could be bimodal
    This is the kind of difference between a general purpose tool (as Maple is) and a specialized tool (R)

Then you recieve such kind of warning a good posture is to assess the mode by changing the method (the algorithm) you used for the first assessment.

There are a lot of algorithms (at least a tenth) to estimate the mode of a distribution.
For your example, R returns these estimations (I give first the name of the method, next the value of the mode,at the end the "bimodal warning" when delivered))

  1. Lientz : 22.5 : could be bimodal
  2. "naive" : 22
  3. Venter  : 22.5 : could be bimodal
  4. Grenander : 22.60124
  5. HSM : 22.5 : could be bimodal 
  6. Parzen with gaussian kernel : 22.5638 : could be bimodal
  7. Tsybakov with gaussian kernel : 22.27502 : could be bimodal
  8. Asselin : 22
  9. Vieu : 22.56381

I don't what is your practice in Statistics, so these results may seem disturbing ?

You probably can understand the calue 22 (which is the one Maple returns) but what about those with a non nul decimal part ?

The truth here is that, as some method say, the distribution "could be bimodal".

In this situation, if you were in position to do other experiments or observations (I interprete your serie as a collection of measurements) you should do them. Increasing the number of observations would probably discriminate the two modes.


To go back to your question : 

  1. Does Maple give the correct answer : NO
  2. Would a statistical softaware give the correct solution (2 modes) : NO neither ...
    but it will warn you that, maybe, your distribution is bimodal (at least if you are aware of this posisbility and if you do not restrict yourself to use only one estimation method)

Hope I have helped you ?

 

1 2 3 4 5 Page 3 of 5