## How to probability density function of a data set ...

Dear All,

I would like to plot the probability density function of a state variable obtained from solving differential equations. I have found that there are functions called "PDF" and "KernelDensityPlot" in the Statistics package, but they really confuse me. Could you please point me out? My code is as follows.

Ps. Is it possible to plot the PDF directly from the solution of dsolve() without discretizing the results?

restart:
with(plots): with(DEtools): with(plottools):with(LinearAlgebra): with(Statistics):

v1:=1: f:=-4: v2:=2.515: omega:=1: epsilon:=0.001: k:=0:
sys:=diff(u1(t),t)=v1*u1(t)-(omega+k*u2(t)^2)*u2(t)-(u1(t)^2+u2(t)^2+3*z(t)^2)*u1(t),
diff(u2(t),t)=(omega+k*u1(t)^2)*u1(t)+v1*u2(t)-(u1(t)^2+u2(t)^2+3*z(t)^2)*u2(t),
diff(z(t),t)=z(t)*(-v1+3*u1(t)^2+3*u2(t)^2+z(t)^2)+epsilon*z(t)*(v2+f*z(t)^4):

t_start:=50: t_end:=300: dt:=0.05: fs:=1/dt:

solA:=dsolve({sys, u1(0)=0.6, u2(0)=0.6, z(0)=0.1},
{u1(t),u2(t),z(t)},
type=numeric, method=rkf45, maxfun=0,
output=Array([seq(i,i=t_start..t_end, dt)])):

u1:=solA[2,1][..,2]:
u2:=solA[2,1][..,3]:
z:=solA[2,1][..,4]:

u0:=sqrt~(u1^~2+u2^~2+z^~2):

Phi:=z/~u0:

# How could I plot the probability density of Phi (y-axis) against Phi(x-axis)?

Probability_density_function_plot.mw

Thank you!

Very kind wishes,

Wang Zhe

## Bug in Probability

Maple 2016

Let us consider

with(Statistics);
U := RandomVariable(DiscreteUniform(-10, 10)):
V := RandomVariable(DiscreteUniform(-10, 10)):
Probability(U^2-V^2 <= 1/9, numeric);
0.

, whereas a positive number greater than 1/21 is expected.

## How to get a numerical answer using the erf functi...

I didn't even know the erf function existed until doing this problem. I looked up how to use it, so I tried plugging in the explicity form into Maple, hoping it'd solve it, but it just spit back out the erf function.

I am trying to get a number answer.  A decimal. Because this is calculating a probability. How do I get Maple to give me a number here? Thanks!

## Bug in Mode

Maple 2016

Let us consider

Statistics:-Mode(Binomial(n, p));
floor((1 + n) p)

Up to Wiki, the output is not correct. Simply no words.

## Bug in Probability

Maple 2016
restart; with(Statistics):
X := RandomVariable(Normal(0, 1)): Y := RandomVariable(Uniform(-2, 2)):
Probability(X*Y < 0);

crashes my comp in approximately 600 s. Mma produces 1/2 on my comp in 0.078125 s.

## Mean Value Analysis Optimization Problem...

I am trying to find the optimal routing probabilities in a Maple procedure where the Mean Value Analysis is used to compute the queueing values. The Maple code is below. It first tries to compute the visit ratios where the probability routing values are the decision variables. There is one specified constraint on the sum of the probability decision variables.

restart;
interface(warnlevel=0): interface(displayprecision = 4): with( plots ):
with(linalg):with( Optimization ); with(Student[NumericalAnalysis]):
[ImportMPS, Interactive, LPSolve, LSSolve, Maximize, Minimize,

NLPSolve, QPSolve]
f:=proc(x1,x2,x3)
global T,lambda,nq,u;
local i,j,pop,Sum;
n:=3;N:=2;M:=3;
#
# Gauss-Seidel iterations
#
A:=Matrix([[1,-x1,-x2],[0,1,-x3],[0,0,1]]);
b:= Vector([1,0,0]);
v := IterativeApproximate(LinearAlgebra:-Transpose(A), b, initialapprox = Vector([1, 3/4, 3/4]), tolerance = 10^(-3), maxiterations = 20, stoppingcriterion = relative(infinity), method = gaussseidel);
mu:=array(1..n,[2.0,1.0,1.0]);
nq:=array(1..M,[0,0,0]);# must initialize queue lengths
for i from 1 to N do

pop:=i;
for j from 1 to M do # mean waiting times
T[j]:=t[j]*(1 + nq[j]) od;
Sum := 0.0;
for j from 1 to M do # mean cycle time
Sum := Sum + v[j]*T[j] od;
for j from 1 to M do #compute the throughputs
lambda[j] := (v[j]*pop)/Sum od;
for j from 1 to M do #compute the queue lengths
nq[j]:= lambda[j]*T[j] od;
for j from 1 to M do #compute the utilizations
u[j]:= lambda[j]*t[j] od;
od;
RETURN(lambda[1]);
end proc;

proc(x1, x2, x3) ... end;

sol := Optimization:-NLPSolve(f, {}, {proc (x1, x2, x3) options operator, arrow; x1+x2+x3-5/3 end proc}, 0 .. 1, 0 .. 1, 0 .. 1, initialpoint = [.75, .25, .6667], assume = nonnegative); 1;

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

I am not sure why I get the error message

## How to get reverse probability?...

When defining a plain standard distributed stochastic variable X, and can find the probability of X <= 0.6 using the Probability function, but how can I get the value for a certain probability, as is done with the fsolve function for example below.

However, the fsolve used to defined Prev above appears to be a bad way to do it, since the Prev function can't for example plot.

Is there some build in way of doing reverse of Probability for a stocastical variable ?

## Walk in random direction...

Hei, I'm trying to create a random walk in the plane, with constant step length (=1) and the angle between two consecutive steps are decided by a probability density function. I just can't seem to find out how I should implement the density function into my worksheet.

The probability density function is: p(phi)=(1/4)*cos(phi/2), on the interval [-Pi,Pi].
And  I think i managed to do it by selecting a random angle, but don't know how to generate a random angle given this probability function. Any ideas? It'd be highly appreciated!

## Probability density normalization...

Probability_density_normalization.mw

In this code I'm trying to separately normalize two independent probability densities and then combine them to get the joint probability density that's normalized and then use it to calculate the probability that the two variables are equal. fD(x) is a Gaussian divided by x^2 and fA(x) is a Gaussian. The first problem occurs when I'm checking the normalization of the joint probability density by doing the double integral over all space for fD(x)*fA(y)dxdy, I get weird vanishing number when the parameter "hartree" takes a certain value, namely 27.211. If I change hartree to 27 or 1 or 2 it all worked, but 27.211 is not good. Also later when I do a single integral over all space for fD(x)*fA(x)dx to get the probability that these two are equal, I find the result is dependent on hartree. This hartree thing is a unit conversion in my physical problem and in principle should not interfere with either the normalization or the probability result at all. I suspect this is a coding bug but I can't find what it is. I'd appreciate any input.

Thank you very much!

Edit: I found out that the problem with the double integral normalization may have something to do with the discretization for numerical evaluation of the integral, since if I change the lower bound to 1/hartree and upper to 10/hartree then it's fine, however if I use lower bound at 1/hartree and upper at 5/hartree it doesn't work, although the distribution has no value between 5/hartree and 10/hartree. However after this is fixed I still have the problem with the single integral over all space for fD(x)*fA(x)dx changing with hartree. Well as a probability I would expect the integral to be bound between 0 and 1, but since it almost linearly depends on hartree, at hartree around 27 I would get the integral value to be about 25, which doesn't make sense. In fact, I now suspect it is not Maple, but my calculation of the probability of the two random variables taking the same value is wrong, I'd appreciate it very much if someone can confirm this.

## Probability of a situation occurance ...

I'm drawing a complete blank on this.

The probability of a branch breaking in a particular tree as someone climbs is 2%.  After the person has climbed over 40 branches, did any of them break?  I'm trying to set that up as a simulation in Maple but I'm not sure about how to do it.

I guess if I generate a random number and it falls under .02 it breaks, is that how I determine that?

randomize():

n := 0;
for i to 40 do
ran := evalf(rand()*10^(-12));
if ran < 0.2e-1 then n := n+1: end if:
end do:
print("Branch broke", n, "times")

I'm not sure I've got the probability right.

## Generating sampling points with a known distribut...

Hi,

Suppose I have a probability distribution function f(x). I would like to generate points in a

(1) given interval, can be up to infinity

(2) number of points

(3) satisfy the known probability distribution function f(x)

(4) possible other specifications

Is there any build-in option in maple I could use? Or I have to write something manually?

Thank you very much!

## A matter of notation...

Dear Maple users

I am unsure how to handle events and their probabilities in Maple. Let's say I know that an event A has the probability say 0.3 and another event B has probability 0.8. I would like to make the following assignments:

P(A):=0.3:

P(B):=0.8

and maybe defining the conditional probability:

P(A|B):=0.55

but I am not allowed to do so in Maple because if will regarded as a function definition. My purpose is to make simple calculations with those probabilities for example:

P(C):=P(A)*P(B)  etc.

My problem is therefore more of a notational problem than a mathematical one. I hope someone can advice me on a proper setup. I am using 2D math notation, by the way. I could of course name the variables containing the probability using simple names like X1, X2, etc., but then I need constantly to remember what they really mean. The above assignments would be much better, because they are easier to handle mentally.

Regards,Erik

## compute the cdf of a piecewise probability distrib...

This should be trivial but I am not able to figure out the right syntax to execute it

The pdf is given by :

f_X(x)={ 1/25 *x, 0<=x<5

2/5 -x/25, 5<=x<10

0, otherwise

I have tried to use the "CumulativeDistributiveFunction" so far

## Expected Value of a function...

Hello people in Mapleprimes,

I want to do calculations about an expected value.

For example,

expectedvalueofc.pdf

that is, \int_a^b c d \Omega

where, \Omega is a distribution function of c.

But, I don't know whether this equation can be put into worksheet.

Is there no way other than writing this as

\int_a^b c g(c) dc,

where g shows a probability density function.

Can I use Statistics[ExpectedValue] to calculate expected value with a general distribution function, not specified

to normal distribution?

I hope you will give me some hints.

Best wishes.