Items tagged with distribution distribution Tagged Items Feed

Erik Postma from Maplesoft has recently posted a Maple application going through some of the ways how to produce random numbers in Maple: http://www.maplesoft.com/applications/view.aspx?SID=153662&P=TC-4444

Doing a fair bit of work involving random numbers (i.e. simulations) myself, I was naturally drawn to it. The document is a nice extension of the relevant help pages, making practical use of some of the procs in RandomTools easier and clearer.

One particular issue I recently ran into, however, is not covered. I needed a random generator for a particular custom pdf (not one of the already recognized pdfs) that happens to be 3-dimensional. And I needed to do this >twice<, with the second generation using the result of the first one as a parameter. the problem is not unlike the last example in the help page for Statistics:-Sample.

Here is what I did:

The pdf looks graphically like this (and note that in this plot it is not normalized):

plot3d(Triangle63:-Trianglef(CrystalAngle,DeflectionAngle),\
       CrystalAngle=-0.8..0.2,DeflectionAngle=-0.6..0.6,projection=0.7);

The first set of random numbers is generated by setting CrystalAngle to a fixed number (-0.1 in this case), normalizing the pdf to the integral over DeflectionAngle and then generating the random numbers. Note that evaluation of the integral is relatively time-consuming (seconds).

trianglepdf:=(CrystalAngle,DeflectionAngle) ->\
             Triangle63:-Trianglef(CrystalAngle,DeflectionAngle)/\
             'evalf(Int(Triangle63:-Trianglef(CrystalAngle,Da),Da=-0.6..0.6,method = _d01akc))':
triangle0pdf:=(t) -> trianglepdf(-0.1,t):
Cr1:=Distribution(PDF=triangle0pdf);
Y:=RandomVariable(Cr1);

samples:=Sample(Y,[1..600],method=[envelope,range=-0.6..0.6]);

Histogram(samples,view=[-0.6..0.6,default],labels=["Deflection Angle  (µrad)","Counts/bin"]);

This looks as it should; and the random number generation is reasonably fast once the initial setup has been done (i.e. the time used scales only relatively weakly with the count of random numbers generated.

The issue now is the second stage, where for each new random number, the number generated above is a new input parameter for CrystalAngle. So I can no longer define the pdf once since the pdf is different for each number generated.

The pedestrian way to do this is like the following:

samples2:=Vector(1..numelems(samples),datatype=float):
CrystalOffset:=0.3;

for ii from 1 to numelems(samples) do
  triangleiipdf:=(t) -> (trianglepdf(samples[ii]-~CrystalOffset,t)); # subtract (VR peak -> 0)
  Crii:=Distribution(PDF=triangleiipdf);
  samples2[ii]:=evalf(Sample(RandomVariable(Crii),1,method=[envelope,range=-0.6..0.6])[1]+CrystalOffset);
end do:

This works, but only sort-of. First, it is very slow, since the whole sertup happens for each number generated (Sample()). Secondly, it tends to hang at a certain number of steps through the loop. The value of ii where it hangs is arbitrary and changes with the other parameters and CrystalOffset. It is however consistent for identical runs.

The last example in help for Statistics:-Sample would indicate that one needs to setup the pdf as a function of t as well as the parameter, call Distribution() only once and then assign the value to the parameter and call Sample() to get the random number(s) drawn from the pdf with the parameter being set to the wanted value. While this works in the example which uses a built-in pdf, I find that I cannot make it work with my pdf. Either the parameter gets ignored (i.e. stuck at the first value) or the thing runs just as slowly as my procedure above.

Ultimately I programmed my own random generator for an arbitrary pdf which, while not as efficient as Maple's built-in generator, does produce the expected set of random numbers for the 2nd path in a reasonable time. It is a simple-minded envelope-rejection method, that works fo reasonably well-behaved pdfs:

Rarbit:=proc(pdf,xmin,xmax,pdfmax)
local dx:=xmax-xmin;
local x1,x2;
x1:=RandomTools[MersenneTwister]:-GenerateFloat64()*dx+xmin; # pick location on x axis
x2:=RandomTools[MersenneTwister]:-GenerateFloat64()*pdfmax; # y axis, probability of returning x1

while (x2>evalf(pdf(x1))) do # if above pdf, reject
  x1:=RandomTools[MersenneTwister]:-GenerateFloat64()*dx+xmin; # new x-axis value
  x2:=RandomTools[MersenneTwister]:-GenerateFloat64()*pdfmax; # check value
end do; # eventually we'll succeed and return one.

x1;
end proc;

Using this routine I get my second and final set of randome:

for ii from 1 to numelems(samples) do
  newpdf:=(DefAng) -> Triangle63:-Trianglef(samples[ii]-CrystalOffset,DefAng);
  samples2[ii]:=Random:-Rarbit(newpdf,-0.6,0.6,28);
  DocumentTools:-SetProperty(ProgressBar,'value',ii,refresh);
end do:
                      CrystalOffset := 0.3
Histogram(samples2+~samples,view=[-0.6..0.6,default]);

 I would be interested in Erik's comment on this.

Mac Dude

 

how to plot tail of normal distribution

would like to see the shape of tail of distribution

Dear Maple experts,

I would like to generate population data that is the best possible approximation of a multivariate normal distribution with a specified covariance matrix and vector of means. I do not want to draw a sample from a multivariate distribution, but I want the population values itself which are approximately multivariate normal distributed. The size of the datamatrix should be limited, otherwise I could draw a huge sample from a multivariate normal distribution. For instance, I would like to generate a 200 by 6 data matrix that is the best (or at least good enough) approximation of a MVN distribution. For a bivariate normal distribution one could calculate the probalities of a grid by integrating the density, but for six variables that seems undoable. 

Before trying the invent the wheel again, I think I will ask this question to experts, because it is unlikely that there is no already existing algorithm that does the job pretty well.

Thanks in advance,

Harry Garst

Hello, Please how do I compute cdf of student t distribution in maple Tξ+1(.). I have a function that i nvolves student t distribution but finding it difficult to compute student t in the funcion. I am new to maple.

Hi everyone

I'm currently working on some mandate distribution using "Jefforson's Method" but I have run into some problems.

The general form of the calculations I do is as follows:

d:=fsolve(m = floor(v1/x)+floor(v2/x), x)

But in the case of m=5, v1=4969 and v2=208 it does not work. If I change v1 a bit it works as a charm but when 
4960=<v1=<4969 it does not.

Can any of you figure out why?

 

The equation surely has a solution (well, a lot of solutions). I can figure some out just by estimating and trying. Furthermore, wolfram alpha easily gives me several solutions:

http://www.wolframalpha.com/input/?i=5+%3D+floor%284969%2Fx%29%2Bfloor%28208%2Fx%29

So how come I cannot get Maple to solve it?

 

Thanks in advance!

Hello,

would you please help me how can i introduce a probability distribution function to maple in document mode?

I want to calculate integral of x f(x)dx, while I want maple to know f(x) is a probability distribution function.

I do not have any assumption about f(x)(for example normal or exponential distriburion)

Thanks

F := ((1/6)*z+(1/3)*z^2)/(1+(1/6)*z+(1/3)*z^2);

Dist := subs(z=t,F);

RealDist := Distribution(CDF = unapply(Dist(t-1),t));

X:=RandomVariable(RealDist);

CommonDensity := PDF(X,t);

 

F do not have D(t) but density have D(t) , what is D(t)?

I am trying to get Maple to calculate for me.. Let X and Y have a trinomial distribution with parameters n=3, p1=1/6 and p2=1/2. I am supposed to find E(X), E(Y), Var(X), Var(Y) and Cov (X,Y). Thank you to anyone who can assist with this. And if  you can help me understand the concept behind it.. Im having some troubles with bivariate distributions in my Stats course :( Thanks!

Hi all,

I have been trying to plot in Maple a Beta Prime Distribution using the Statistics package. I have define it through its density function and its range with the command

U := Distribution(PDF = (proc (x) options operator, arrow; x^(alpha-1)*(1+x)^(-alpha-beta)/Beta(alpha, beta) end proc), Support = 0 .. infinity)

and then assigned it to a random variable Z with the command

Z := RandomVariable(U)

Now I wanted to plot the density...

Hello,

I am trying to get a shaded area in my plot but I could not.

First we solve ODE without randomness:

ode := diff(U(t), t) = -(A+B*U(t))*U(t);

Then we add randomness to ODE and solve:

ode2 := diff(U(t), t) = -(A+r(t)+B*U(t))*U(t);

A with randomness for r in R=( - 0.0001/365, 0.0001/365) is:

A(t,r)= A+r

Where A is constant =  0.0001/365

We plot both solution. For the plot I...

Custom distribution
restart;
with(Statistics):
Density := 2*exp(-t^2)*sqrt(Pi);
Dist := int(Density, t);

a.

RealDist := Distribution(Dist(t-1));
Error, (in Statistics:-Distribution) unexpected argument(s): Pi(t)*(erf(t))(t)

b.

RealDist := Distribution(PDF = (t->2*exp(-t^2)*sqrt(Pi)), CDF = (t->Dist(t)), Mean = 1);

RealDist := Distribution(PDF = (t->2*exp(-t^2)*sqrt(Pi)), CDF = (t->Dist(t)));
X:=RandomVariable(RealDist(t-1));

To plot the density function of the continuous uniform distribution on [-1,1], my initial attempt was: 

plot(Statistics:-PDF(Statistics:-RandomVariable(Uniform(-1,1)),x), x = -1.1 .. 1.1);

See plot below.


But I wanted something more like the wikipedia image (without the labels, naturally):

See plot below.

In words, I expected a horizontal line on the left of x=-1 and on the right of x=1 (at y=0), and I expected no vertical line at the x=-1 and x=1 points ...

can i distrubute some random particle with an arbitrarity function in maple, for example p(r,theta,phi)=p0*sin(theta)/r, 0.01<r<1, 0< theta<pi, 0<phi< 2*pi. For such distribution we expect many particle in equator.

many thanks

I have been using random numbers in other applications than Maple. Usually there is a function, which will give a pseudo random real number between 0 and 1. When I looked for it in Maple I got quite confused, because there are a lot of different options here - obviously because Maple can deliver random numbers/objects in many ways, even following a certain distrubution. I found out it doesn't work by just using rand(), since it is always starting with the same value. Then I found the command randomize(...

Dear Maple Users

Let the density of electronic states be denoted

1 2 3 Page 1 of 3