Items tagged with sample

Feed

Hi everyone,

I've not been able to figure this one out. Say I have an expression dependant on two random variables, like this:

C:=A+B

where A and B are randomvariables, each following a specific distribution.

If I ask for a Sample of C

Sample(C,1)

Maple will sample A, sample B and compute C. But say that A and B are correlated (with a cc of 0.8). How do I define this?

Thanks in advance

 

Hi,

I did some hypothesis testing exercises and I cross checked the result with Maple. I just used following vectors for an unpaired test

a := [88, 89, 92, 90, 90];
b := [92, 90, 91, 89, 91];

I ended up with the following solution:

HFloat(1.5225682336585966)
HFloat(-3.122568233658591)
for a 0.95 confidence interval.

 

Using

TwoSampleTTest(a, b, 0, confidence = .95, summarize = embed)

and

TwoSampleTTest(a, b, 0, confidence = .975, summarize = embed)

I get following results:

-2.75177 .. 1.15177

-3.13633 .. 1.53633

respectively. I can not explain the discrepancy.

 

Best regards,

Oliver

 

PS:

Maple Code in case files won´t be attached.

 

 

Unpaired t Test
restart;
Unpaired test-test dataset
a := [88, 89, 92, 90, 90];
b := [92, 90, 91, 89, 91];
The se² estimate is given by:
se²=var(a)+var(b)+2*cov(a*b)=var(a)+var(b)
se²=
sigma[a]^2/Na+sigma[b]^2/Nb;
with Na, Nb being the length of vector a and b respectively.
                             2                              2
  sigma[[88, 89, 92, 90, 90]]    sigma[[92, 90, 91, 89, 91]]
  ---------------------------- + ----------------------------
               Na                             Nb             
sigma[a]^2;
 and
sigma[b]^2;
 are approximated by
S[a]^2;
 and
S[b]^2;
                                             2
                  sigma[[88, 89, 92, 90, 90]]
                                             2
                  sigma[[92, 90, 91, 89, 91]]
                                           2
                    S[[88, 89, 92, 90, 90]]
                                           2
                    S[[92, 90, 91, 89, 91]]
with
S[X]^2;
 defined as
S[X]*`²` = (sum(X[i]-(sum(X[j], j = 1 .. N))/N, i = 1 .. N))^2/N;
                                 2
                             S[X]
                                                 2
                      /      /         N       \\
                      |      |       -----     ||
                      |  N   |        \        ||
                      |----- |         )       ||
                      | \    |        /    X[j]||
                      |  )   |       -----     ||
                      | /    |       j = 1     ||
                      |----- |X[i] - ----------||
                      \i = 1 \           N     //
             S[X] ᅡᄇ = ----------------------------
                                   N              
with(Statistics);
Sa := Variance(a);
                   HFloat(2.1999999999999993)
Sb := Variance(b);
                   HFloat(1.3000000000000003)
Now we are ready to do hypothesis testing (0.95).
We have (with k=min(Na,Nb)=5):
C = mean(a)-mean(b); Deviation := t_(alpha/a, k-1)*se(Sa/k-Sb/k);
c := Mean(a)-Mean(b); deviation := 2.776*sqrt((1/5)*Variance(a)+(1/5)*Variance(b));
                  HFloat(-0.7999999999999972)
                   HFloat(2.3225682336585938)
upperlimit := c+deviation; lowerlimit := c-deviation;
                   HFloat(1.5225682336585966)
                   HFloat(-3.122568233658591)

Execution of built in student test
TwoSampleTTest(a, b, 0, confidence = .95, summarize = embed);

 

 

Is there a difference between these two? 

with(Statistics):

Sample(Normal(0,1),100)

Sample(RandomVariable(Normal(0, 1)), 100)

 

Hello,

I've been trying to collect points at certain intervals in my plot, however, the sample option doesn't seem to be working. Can anyone lend me a hand? Any help is greatly appreciated! Thank you in advance!

Kind regards,

Gam

with(plots):

a := 1.501*10^9:

Th := sqrt(4*Pi^2*a^3/(G*(Mh+Msat))):

HyperionOrbit := proc (`θIC`, `ωIC`, n) local a, Mh, Msat, G, e, beta, M, Eqns, ICs, soln; global `ωH`, Th; a := 1.501*10^9; Mh := 5.5855*10^18; Msat := 5.6832*10^26; G := 6.67259/10^11; e := .232; beta := .89; M := Mh+Msat; Eqns := diff(theta(t), t) = omega(t), diff(omega(t), t) = -G*Msat*beta^2*(xH(t)*sin(theta(t))-yH(t)*cos(theta(t)))*(xH(t)*cos(theta(t))+yH(t)*sin(theta(t)))/(xH(t)^2+yH(t)^2)^2.5, diff(xH(t), t) = vxH(t), diff(vxH(t), t) = -G*M*xH(t)/(xH(t)^2+yH(t)^2)^(3/2), diff(yH(t), t) = vyH(t), diff(vyH(t), t) = -G*M*yH(t)/(xH(t)^2+yH(t)^2)^(3/2); ICs := xH(0) = a*(1+e), yH(0) = 0, vxH(0) = 0, vyH(0) = sqrt(G*M*(1-e)/(a*(1+e))), theta(0) = `θIC`, omega(0) = `ωIC`; soln := dsolve({Eqns, ICs}, numeric, maxfun = 0); plots:-odeplot(soln, [theta(t), omega(t)/`ωH`], 0 .. n*Th, numpoints = 2000, labels = ["θ(t)","ω(t)/ωH"], axes = boxed, style = point, sample = [seq(i, i = 0 .. n*Th, Th)]) end proc:

HyperionOrbit(.5, 1.8*`ωH`, 10)

Error, (in plot/options2d) unexpected option: sample = [0, 1876321.326, 3752642.652, 5628963.978, 7505285.304, 9381606.630, 11257927.96, 13134249.28, 15010570.61, 16886891.93, 18763213.26]

 

``

``

 

Download Poincare_section.mw

bia Man

Hi all,

 

I want to create a (set of) samples from an own distribution. The distribution itself is a normal distribution, but with some 'gaps', as shown below:

plot(exp(-(t-10)^2/(2*(0.5e-1*10)^2))*piecewise(t<9,0,t<9.9,1, t<10.1,0, t<11,1), t = 8.5 .. 11.5)

(It is supposed to represent the distribution of discrete resistor values (e.g. 10 ohms), where not only the values above the tolerance (e.g. 10%) are removed, but also those who can be sold as higher precision resistors (eg. 1%). The mean, lower tolerance and higher tolerance should be variables).

Because of the fact that in order to be able to use the Sample function, the PDF should be twice differentiable, the 'truncation' in the middle is realised using logistic functions.

I am able to create said distribution, and produce a DensityPlot, but when I try to take a (even a single) Sample from a RandomVariable from that distribution, Maple keeps calculating and no sample is produced.

GetDist := (mu, alpha, beta) -> Distribution(PDF = unapply(LotRPDF(t, mu, alpha, beta), t)):
T := GetDist(10, .1, 0.1e-1);
DensityPlot(T);
RandomR := RandomVariable(T):
Sample(RandomR, 1);

Where LotRPDF is the function depicted above, mu is the mean (10), alpha is the maximum tolerance (10%), , beta is the minimum tolerance (1%)

Am i overlooking something simple? Or is there a more fundamental error in my approach? I've added the maple sheet with more details.

RandomVariableSample.mw

 

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

 

Dear Maple users

In a short Maple program I am dealing with a collection of random samples from a uniform distribution. My problem is that every time I call a procedure creating a sample, I get the exact same one! Every run needs a different seed. How do I get a new 'random' sample? I could probably let Maple create all from scratch calling the Sample command with a two-dimensional range, but this will requre a big memory ...

 

Here is part of my code:

If I have  the following..or different numbers, but 2 different arrays I need to calculate sample mean & variance

Array A:=([1,1,1,0,0,3,2,4])

and

Array B:=([2,1,0,1,3,3,4,1])

what formula do I use to force Maple 16 to calculate the sample mean and the sample variance?

Hi everyone,

 

Is it possible to create a continuous time wiener process in Maple?

When we create a wiener process and then a sample path, the values of the sample path are saved in a table, so if I call W(1.5) there is no answer.

How can this be solved?

 

Thanks.

Page 1 of 1