longrob

Mr. Robert Long

1379 Reputation

18 Badges

12 years, 124 days
Leeds, United Kingdom

Social Networks and Content at Maplesoft.com

Business Consultant Data Scientist Statistician

MaplePrimes Activity


These are replies submitted by longrob

@EhsanKayo , what range for t are you interested in ?

@EhsanKayo , what range for t are you interested in ?

There are formatting problems with my post above - it's not wrapping the lines, although you should be able to just copy/paste.

Anyway, my worksheet is attached.


eq1:= B[n1] = -A[n1]*R[s]^np/R[s]^(-np):

eq2:= B[n2] = -(-mu[r]*A[n2]*R[r]^np+mu[r]*A[n2]*R[r]^np*np^2-M[n]*R[r])/(mu[r]*R[r]^(-np)*(-1+np^2)):

 

eq3:= A[n2] = (-np*mu[R[m]]*A[n1]*R[m]^np+np^3*mu[R[m]]*A[n1]*R[m]^np+np*mu[R[m]]*B[n1]*R[m]^(-np)-np^3*mu[R[m]]*B[n1]*R[m]^(-np)-np*mu[R[m]]*B[n2]*R[m]^(-np)+np^3*mu[R[m]]*B[n2]*R[m]^(-np)+M[n]*R[m])/(np*mu[R[m]]*R[m]^np*(-1+np^2)):

 

eq4:= A[n1] = (mu[R[m]]*B[n1]*R[m]^(-np)-mu[R[m]]*B[n1]*R[m]^(-np)*np^2-mu[R[m]]*A[n2]*R[m]^np+mu[R[m]]*A[n2]*R[m]^np*np^2-mu[R[m]]*B[n2]*R[m]^(-np)+mu[R[m]]*B[n2]*R[m]^(-np)*np^2-M[n]*R[m])/(mu[R[m]]*R[m]^np*(-1+np^2)):

 

solve({eq1,eq2,eq3,eq4},[B[n1],B[n2],A[n1],A[n2]]);

(1)

 


Download AliKhan_01.mw

There are formatting problems with my post above - it's not wrapping the lines, although you should be able to just copy/paste.

Anyway, my worksheet is attached.


eq1:= B[n1] = -A[n1]*R[s]^np/R[s]^(-np):

eq2:= B[n2] = -(-mu[r]*A[n2]*R[r]^np+mu[r]*A[n2]*R[r]^np*np^2-M[n]*R[r])/(mu[r]*R[r]^(-np)*(-1+np^2)):

 

eq3:= A[n2] = (-np*mu[R[m]]*A[n1]*R[m]^np+np^3*mu[R[m]]*A[n1]*R[m]^np+np*mu[R[m]]*B[n1]*R[m]^(-np)-np^3*mu[R[m]]*B[n1]*R[m]^(-np)-np*mu[R[m]]*B[n2]*R[m]^(-np)+np^3*mu[R[m]]*B[n2]*R[m]^(-np)+M[n]*R[m])/(np*mu[R[m]]*R[m]^np*(-1+np^2)):

 

eq4:= A[n1] = (mu[R[m]]*B[n1]*R[m]^(-np)-mu[R[m]]*B[n1]*R[m]^(-np)*np^2-mu[R[m]]*A[n2]*R[m]^np+mu[R[m]]*A[n2]*R[m]^np*np^2-mu[R[m]]*B[n2]*R[m]^(-np)+mu[R[m]]*B[n2]*R[m]^(-np)*np^2-M[n]*R[m])/(mu[R[m]]*R[m]^np*(-1+np^2)):

 

solve({eq1,eq2,eq3,eq4},[B[n1],B[n2],A[n1],A[n2]]);

(1)

 


Download AliKhan_01.mw

please can you post your code here

@lp15@cam.ac.uk 

I just remembered and I'm pretty sure there was a similar question here on mapleprimes about a month ago. I can't seem to find it now, though :( 

Are you on windows ? Has the machine been updated/patched since Friday ?

@lp15@cam.ac.uk 

I just remembered and I'm pretty sure there was a similar question here on mapleprimes about a month ago. I can't seem to find it now, though :( 

Are you on windows ? Has the machine been updated/patched since Friday ?

@Markiyan Hirnyk 

Obviously it depends on the particular functions, but since the poster didn't specify any, it's hard to know what he wants. I find it's usually best to start simple, particularly when the poster says he is new to Maple.

It may be hard to provide code that works for all possible scenarios (sinusoids, complex-valued functions etc)

Edit: mapleprimes cutting off part of the reply

@Markiyan Hirnyk 

Obviously it depends on the particular functions, but since the poster didn't specify any, it's hard to know what he wants. I find it's usually best to start simple, particularly when the poster says he is new to Maple.

It may be hard to provide code that works for all possible scenarios (sinusoids, complex-valued functions etc)

Edit: mapleprimes cutting off part of the reply

@marbovski  That's because = is not the assignment operator, := is the assignment operater

@marbovski  That's because = is not the assignment operator, := is the assignment operater

Please can you clarify

1. mean arrival rate : mean number of events per unit time ?

2. average interarrival time: mean time between consecutive events ?

@JoeKing, i assume you mean the normal probability plots ?

Actually I think they are very good !
At first sight there appears to be a bit of a departure in the tails, but even when directly sampling a true
normal distribution, the points don't all lie exactly on the line at the tails, just due to random sampling.
To see this:

restart:
with(Statistics):
RV:=RandomVariable(Normal(0, 1)):
X:=Sample(RV,2000):
p1:=plot(PDF(RV,x),x=-4..4):
p2:=Histogram([X],bincount=15):
p1:=plots:-display(p1,p2):
p3:=ProbabilityPlot([X],Normal(0,1),thickness=2,color=BLACK):
plots:-display(Array([p1,p3]));


First note that the bulk of the points lie very close the the line. As for the departures in the
tails, one important thing to check here is the scale of the y-axis. Note that the departures
in the vertical direction are quite small - even the largest is around 5% of the "true" value.
Now, looking at the normal probability plots in my original post above, the largest departure
is also of around 5%, while the bulk of the points are very close to the line.

Now look at a sample taken from a skewed distribution.
RV:=RandomVariable(Gamma(5,10)):
p1:=plot(PDF(RV,x),x=0..150):
X:=Sample(RV,1000):
p2:=Histogram([X],bincount=15):
p1:=plots:-display(p1,p2):
p3:=ProbabilityPlot([X],Normal(50,sqrt(250)),thickness=2,color=BLACK):
plots:-display(Array([p1,p3]));


First note (from the pdf plot) that this distribution is not *greatly* skewed, but the normal
probability plot clearly shows two things. First the points themselves show a systematic
curvature throughout the whole range; and second, the largest departures are
100%+ at the left side and 20%+ on the right. You don't see features like this in the plots
I posted above.

 Edit: Add histogram to pdf plot



 

In a reply to Markiyan Hirnyk's recent question on this topic, Alec Mihailovs posted solutions 
for n=1..7 for the mean distance between random points in an n-cube.
I also posted some code to simulate the value for arbitrary n.

In this post I'm going to show the results of simulations for n=1..10, along with histograms 
of the results (with superimposed normal pdf curve) and normal probability plots, since the
simulations rely on the central limit theorem. % errors, for n=1..7 where the result is known,
and 95% confidence intervals for those where the result is not known (by me), are also shown.

ncubem:=proc(n::posint,N::posint,m)
# procedure simulates the mean value between two random points A and B
# in a n-cube of face length m, N times.
local A,B,s,i,RV,dim:
RV := Statistics:-RandomVariable(Uniform(0,m)):
s:=0:
for dim from 1 to n do
A[dim]:=Statistics:-Sample(RV,N);
B[dim]:=Statistics:-Sample(RV,N);
end do:
for i from 1 to N do
s:=s+sqrt(sum((A[dimm][i]-B[dimm][i])^2,dimm=1..n))
end do:
s/N;
end proc:

unitncubesim:=proc(n::posint,N::posint,monte::posint,exact)
# procedure samples the mean value of a unit n-cube N times
# exact is an expression for the known exact value. If this
# is provided, the statistical error is computet. If it is not
# then a 95% confidence interval for the mean is computed.
local sample,i,mu,sigma,qtl,ciL,ciU,p1,p2,p3,XX,A,pr1;
sample:=seq(ncubem(n,monte,1),i=0..N):
mu:=Statistics:-Mean([sample]);
sigma:=Statistics:-StandardDeviation([sample]):
printf("simulated value for %d-cube: %0.6f",n,mu);

if type(exact,algebraic) then
printf("; known value: %0.6f, %0.4f%% error",evalf(exact),100*evalf(abs(mu-exact)/exact,6));
else
qtl:=Statistics:-Quantile(Normal(0,1),evalf(1-((1-(95/100))/2))):
ciL:=evalf(mu-(qtl*sigma/sqrt(N)),5):ciU:=evalf(mu+(qtl*sigma/sqrt(N)),5):
printf(" with 95%% confidence interval [%0.6f,%0.6f]",ciL,ciU);
end if;

p1:=Statistics[Histogram]([sample],bincount=15):

XX:=Statistics:-RandomVariable(Normal(mu,sigma)):
p2:=plot(Statistics:-PDF(XX,x),x=mu-(4*sigma)..mu+(4*sigma)):
p1:=plots:-display(p1,p2);
p3:=Statistics:-ProbabilityPlot([sample],Normal(mu,sigma),thickness=2,color=BLACK):
print(plots:-display(Array([p1,p3])));
end proc:




f:=n->evalf(2^n*Int(sqrt(add((x||i)^2,i=1..n))*mul(1-x||i,i=1..n),
[seq(x||i=0..1,i=1..n)])):
#HT Alec Mihailovs for this solution for n=1..7

F:=seq(f(n),n=1..7):

Digits:=15:

t:=time():
unitncubesim(1,2000,500,F[1]);
    simulated value for 1-cube: 0.333516; known value: 0.333333,    0.0549% error

unitncubesim(2,1500,500,F[2]);
    simulated value for 2-cube: 0.521228; known value: 0.521405,    0.0339% error

unitncubesim(3,1000,500,F[3]);
    simulated value for 3-cube: 0.661723; known value: 0.661707,    0.0024% error

unitncubesim(4,900,400,F[4]);
    simulated value for 4-cube: 0.777830; known value: 0.777666,    0.0211% error

unitncubesim(5,800,300,F[5]);
    simulated value for 5-cube: 0.878373; known value: 0.878531,    0.0180% error

unitncubesim(6,700,300,F[6]);
    simulated value for 6-cube: 0.969059; known value: 0.968942,    0.0121% error

unitncubesim(7,650,300,F[7]);
    simulated value for 7-cube: 1.052150; known value: 1.051584,    0.0542% error

unitncubesim(8,600,300,"");
    simulated value for 8-cube: 1.128114 with 95% confidence interval [1.12700,1.12920]

unitncubesim(9,550,300,"");
    simulated value for 9-cube: 1.200337 with 95% confidence interval [1.199100,1.201500]

unitncubesim(10,500,300,"");
    simulated value for 10-cube: 1.267209 with 95% confidence interval [1.266000,1.268400]

time()-t:
%/60;
                  8.0751
 

@Markiyan Hirnyk 

Apparently,

"The expected value of the distance between two random points in a box was stated as
a problem by D.P. Robbins, [5], and calculated by T.S. Bolin, [3]."

....
[3] T.S. Bolis Solution of problem E2629: Average Distance between Two Points in a Box,
Amer. Math. Monthly Vol. 85, No. 4. pp. 277-78, April. 1978.
...
[5] D.P. Robbins Problem E2629: Find the Average Distance between Two Points in a Box,
Amer. Math. Monthly [1977,57].
............
The same method was also used here:
http://www.math.utep.edu/Faculty/moschopoulos/Publications/1999-Distance_Between_Random_Points_in_a_Cube.pdf
As for higher dimensions, the same author (Johan Philip, as I posted initially) provides this
http://www.math.kth.se/~johanph/h45.pdf
Note:"The obtained expressions contain unsolvable integrals"

 Edit: formatting issues

3 4 5 6 7 8 9 Last Page 5 of 13