Home › Forums › Get Help › How do I ...? with Maple

Serial correlation in coin toss

Posted on 2008-06-07 09:17 By alex_01 (Maple Rating 2 102)
Categories:

When I want to simulate and plot the distribution from a random coin toss (normal distribution) I use the following code:

restart:
randomize():
coin:=rand(0..1):
coin_1:=proc(n) seq(coin(),i=1..n) end:
x_1:=seq([coin_1(10)],i=1..500):
x_2:=[seq(numboccur(1,x_1[i]),i=1..500)]:
with(Statistics):
Histogram(x_2,frequencyscale=relative,discrete = true );

what code should I use in order to get a coin with serial correlation. A normal distribution with fat tails? (Prefereably I want to tweak the above code)

 

‹ EdgeConnectivity in GraphTheory package Positivity of a finite double sum ›
alec's picture

Serial correlation

Comment on 2008-06-07 16:23 from alec (Maple Rating 4 1497)

AFAICT, there are no tools in Maple for simulating time series with given autocorrelation. You have to enter formulas yourself. Choosing any specific distribution won't help - they don't depend on previous results.

Alec

 

I think the biggest problem

Comment on 2008-06-07 17:05 from alex_01 (Maple Rating 2 102)

I think the biggest problem I have is to write the if function for the sequence elements.

restart:
randomize():
coin:=rand(0..1):
coin_1:=proc(n) seq(coin(),i=1..n) end:
x_1:=seq([coin_1(10)],i=1..1);

                       [1, 1, 0, 1, 0, 0, 1, 1, 0, 0]
 

randomize():
coin:=rand(1..10):
coin_1:=proc(n) seq(coin(),i=1..n) end:
x_2:=seq([coin_1(10)],i=1..1);

                       [2, 3, 5, 2, 5, 1, 7, 8, 8, 10]
 

I want to create a third sequence x_3 that have the same value as x_1 in the previous period if x_2 in this period is larger than 6 and is random ( 0 or 1) if x_2 in this period is smaller  than 6. 

I can then change the threshold (in this case 6) to control the amount of serial correlation                     

Any suggestions?

 

alec's picture

For example

Comment on 2008-06-07 17:24 from alec (Maple Rating 4 1497)

For example,

randomize():
coin:=rand(0..1):
c10:=rand(1..10):
x1:=['coin()'$10];

                 x1 := [1, 1, 1, 0, 1, 0, 0, 0, 1, 1]

x2:=['c10()'$10];

                x2 := [3, 10, 4, 9, 6, 6, 1, 1, 5, 3]

x3:=[seq](`if`(x2[i]<7,coin(),x1[i]),i=1..10);

                 x3 := [0, 1, 1, 0, 0, 1, 0, 0, 1, 0]

Alec

Thanx alec !! You are the

Comment on 2008-06-07 17:33 from alex_01 (Maple Rating 2 102)

Thanx alec !! You are the man !  It looks like it will do the trick ! I love this Maple forum ha ha

I think I was wrong in my

Comment on 2008-06-07 18:34 from alex_01 (Maple Rating 2 102)

I think I was wrong in my second post when I outlined the problem

x3 should not use x1 as a reference but rather x3(t-1) because otherwise it dont make sense because then we have 2*random and we dont have any serial correlation

restart;
randomize():
coin:=rand(0..1):
c10:=rand(1..10):
x1:=['c10()'$20];
       [5, 10, 2, 7, 5, 6, 10, 1, 1, 2, 2, 7, 1, 7, 3, 1, 6, 1, 4, 4]
x2:=[seq](`if`(x1[i]<4,coin(),x2[i-1]),i=1..20);
Error, recursive assignment

So my problem no is what code I use for the recursive assignment ?

alec's picture

Instead of recursive assignment

Comment on 2008-06-07 20:52 from alec (Maple Rating 4 1497)

This can be done as

randomize():
coin:=rand(0..1):
c10:=rand(1..10):
x1:=['c10()'$20];
      
  x1 := [1, 4, 7, 8, 3, 7, 1, 5, 3, 10, 9, 8, 4, 7, 6, 9, 10, 2, 3, 6]

s:=select(i->x1[i]<4,{$1..20}) union {1,21};

                    s := {1, 5, 7, 9, 18, 19, 21}

x2:=[seq](coin()$(s[i]-s[i-1]),i=2..nops(s));

  x2 := [1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0]

By the way, serial correlation is in both cases, just with different lags - lag 10 in the former case and lag 1 in this one.

Alec

alec's picture

Coin toss with serial correlation

Comment on 2008-06-08 02:13 from alec (Maple Rating 4 1497)

A random sample of length n drawn from Bernoulli distribution with probability of success prob, that has a correlation c with itself shifted back lag steps, can be generated using following procedure,

SampleWithCorr:=proc(prob::And(numeric,satisfies(c->c>=0 and c<=1)),
lag::posint,c::And(numeric,satisfies(c->c<=1 and c>=-1)),n::posint)
local X,B,S,C,s,i;
uses Statistics;
X:=RandomVariable(Bernoulli(prob));
S:=Sample(X,n);
if n<=lag or s=0 then S else
s:=signum(c);
B:=RandomVariable(Bernoulli(abs(c)));
C:=Sample(B,n-lag);
if s=1 then 
for i from lag+1 to n do
if C[i-lag]=1 then S[i]:=S[i-lag] fi od;
else for i from lag+1 to n do
if C[i-lag]=1 then S[i]:=1-S[i-lag] fi od
fi fi; S end:

For example,

S:=SampleWithCorr(1/2,1,-0.7,10000);

                       [ 10000 Element Row Vector ]
                  S := [ Data Type: float[8]      ]
                       [ Storage: rectangular     ]
                       [ Order: Fortran_order     ]

Statistics:-Correlation(S[1..9999],S[2..10000]);

                            -.6973705688

I've made it into a blog entry.

Alec

Again thanx Alec for your

Comment on 2008-06-08 05:15 from alex_01 (Maple Rating 2 102)

Again thanx Alec for your input and effort. I must say I quite impressed by your willingness to help !  

I do appreciate your input

Comment on 2008-06-08 06:11 from alex_01 (Maple Rating 2 102)

I do appreciate your input but I am not sure that I have understod what you have done. Isnt there a simpler way to solve the recursive assignment in the below formula so I (a simple man) can understand what is going on  : -)

restart;
randomize():
coin:=rand(0..1):
c10:=rand(1..10):
x1:=['c10()'$20];
       [5, 10, 2, 7, 5, 6, 10, 1, 1, 2, 2, 7, 1, 7, 3, 1, 6, 1, 4, 4]
x2:=[seq](`if`(x1[i]<4,coin(),x2[i-1]),i=1..20);
Error, recursive assignment

 

alec's picture

A loop

Comment on 2008-06-08 13:10 from alec (Maple Rating 4 1497)

A usual way of avoiding recursive assignments is using a loop,

interface(rtablesize=20):randomize():
coin:=rand(0..1):
c10:=rand(1..10):
x1:=Vector[row](19,i->c10());

  x1 := [9, 8, 7, 6, 6, 7, 8, 10, 2, 1, 7, 1, 2, 10, 8, 6, 2, 1, 8]

x2:=Vector[row](20,[coin()]);

  x2 := [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

for i from 1 to 19 do 
x2[i+1]:=`if`(x1[i]<4,coin(),x2[i]) od:
x2;

     [1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0]

That, generally, is what I did in the SampleWithCorr procedure.

from 1 can be omitted - it is the default value in Maple.

Alec

Again thanx for your input

Comment on 2008-06-08 13:11 from alex_01 (Maple Rating 2 102)

Again thanx for your input Alec. As usual solid as a rock !

I have been working on it myself for a couple of hours (yes I know I am not the sharpest tool in the shed) and this is what I come up with. I have double checked it and it seems to be correct.....At least that what I hope.

 

# A simulation of a positive serial correlated coin toss

# A threshold value of 5 is a random coin (no possitive serial correlation)
# A threshold value of between 5 to 10 gives a coin with possitive serial correlation
# The higher the threshold the more possitive serial correlation we will have

restart;
randomize():
coin:=rand(0..1):
coin_1:=proc(n) seq(coin(),i=1..n) end:
x1:=seq([coin_1(50)],i=1..1);

[0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1,

  1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 0, 0, 0]

c10:=rand(1..10):
x2:=['c10()'$50];

[10, 9, 5, 5, 4, 1, 3, 3, 7, 3, 5, 8, 9, 3, 8, 9, 4, 6, 3, 8, 3, 8, 9, 9, 7,

  6, 4, 10, 10, 7, 4, 2, 7, 1, 8, 6, 6, 5, 4, 2, 6, 4, 2, 5, 4, 4, 5, 4, 3, 6]

threshold:=6: 
F := proc(n)
if n = 1 then return x1[n] end if:
if n <> 1 and x2[n]<=(10-threshold) then return x1[n] end if:
if n <> 1 and x2[n]>(10-threshold) then return F(n - 1) end if:
end:

x3:=seq( F(i), i=1..50);

0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0,

  1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 1, 1, 0, 0, 0

 

alec's picture

Few comments

Comment on 2008-06-08 13:29 from alec (Maple Rating 4 1497)

This way it also can be done. Just few comments.

1. instead of 3 ifs you could use if .. elif .. else .. fi

2. You could add option remember; after proc(n). Then F(n-1) won't be calculated another time.

3. The standard way of referring to the procedure name inside of a procedure is procname, i.e. procname(n-1) instead of F(n-1). Then, if you rename your procedure later, you won't have to change the code inside it.

Alec

alec's picture

Threshold

Comment on 2008-06-08 14:05 from alec (Maple Rating 4 1497)

Just noticed your comment about threshold 5. That doesn't work like that. Positive correlation will be with any value - just different one (because there still will be dependency). For negative correlation you could assign x[i] to1-x[i-1] in cases where you assigned it to x[i-1]. Also, it would be negative correlation with any value. No threshold.

Alec

 

I am no programming expert

Comment on 2008-06-08 13:40 from alex_01 (Maple Rating 2 102)

I am no programming expert so your help is greatly appreciated. To tell you the truth I am struggling to keep my head over the water most of the time. It would have been different if I would have designed the language, then I would know all the tricks he he.  

Your last code was simple and straight forward I think I will use that one !! Again thanx for your help 

 

I have done some work that I

Comment on 2008-06-12 09:38 from alex_01 (Maple Rating 2 102)

I have done some work that I would like to share! 

If you have any ideas for improvements or comments feal free to post

 

########################  Random coin and the number of heads  ################################

# Note that n is the number of coin tosses
# Note that nseq is the number of sequences
# Note that 1 is heads and 0 is tails
# Note that since #heads= 1- #tails the distribution will look the same for tails

restart:
n:=20:                   # parameter 1 
nseq:=1000:       # parameter 2   (for a normal distribution n=1000)

randomize():
coin:=rand(0..1):
coin_1:=proc(n) seq(coin(),i=1..n) end:
x_1:=seq([coin_1(n)],i=1..nseq):
x_2:=[seq(numboccur(1,x_1[i]),i=1..nseq)]:
with(Statistics):
Histogram(x_2,frequencyscale=relative,discrete = true,color=red );

 

#######################  Random coin and sequence lenght  ###############################

# Note that n is the number of coin tosses
# Note that 1 is heads and -1 is tails

restart:
n:=10000:            # parameter 1 

randomize():
coin:=rand(0..1):
coin_1:=proc(n) seq(coin(),i=1..n) end:
L:=[coin_1(n)]:
LL:=[seq(subs(0=-1,L[i]),i=1..n)]: # Transform to a -1 and 1 coin toss
x_1:=remove(Testzero, map(-nops, [ListTools:-Split(`=`,LL, 1)])): # Sequence -1
x_2:=remove(Testzero, map(nops, [ListTools:-Split(`=`,LL, -1)])): # sequence 1

y:=[op(x_1),op(x_2)]:
with(Statistics):
Histogram(y,frequencyscale=relative,discrete = true,color=red );

LLL:=Vector[row]([LL]):
Statistics:-Correlation(LLL[1..(n-1)],LLL[2..n]);

                               -0.003387130345

 

#######################  Serial correlated coin and the number of heads ############################

# A threshold value of 0 is a random coin (no positive serial correlation)
# A threshold value of between 0 to 10 gives a coin with positive serial correlation
# The higher the threshold the more positive serial correlation the coin will have
# Note that n is the number of coin tosses
# Note that nseq is the number of sequences

restart;
n:=20:                    # parameter 1
threshold:=7:       # parameter 2
nseq:=1000:        # parameter 3

randomize():
coin:=rand(0..1):
coin_1:=proc(n) seq(coin(),i=1..n) end:
x1:=seq([coin_1(n)],i=1..nseq):
c10:=rand(1..10):
x2:=seq(['c10()'$n],i=1..nseq):

F := proc(n_1,n_2)
if n_2 = 1 then return x1[n_1,n_2] end if:
if n_2 <> 1 and x2[n_1,n_2]<=(10-threshold) then return x1[n_1,n_2] end if:
if n_2 <> 1 and x2[n_1,n_2]>(10-threshold) then return F(n_1,(n_2 - 1)) end if:
end:

x_1:=seq( [seq( F(n_1,n_2), n_2=1..n)],n_1=1..nseq):

x_2:=[seq(numboccur(1,x_1[i]),i=1..nseq)]:
with(Statistics):
Histogram(x_2,frequencyscale=relative,discrete = true );

 

##########################   Serial correlated coin and sequence lenght  ################################

# A threshold value of 0 is a random coin (no positive serial correlation)
# A threshold value of between 0 to 10 gives a coin with positive serial correlation
# The higher the threshold the more possitive serial correlation the coin will have
# Note that n is the number of coin tosses

restart;
n:=50000:       # parameter 1
threshold:=4:   # parameter 2

randomize():
coin:=rand(0..1):
coin_1:=proc(n) seq(coin(),i=1..n) end:
x1:=[coin_1(n)]:
c10:=rand(1..10):
x2:=['c10()'$n]:
 
F := proc(n)
if n = 1 then return x1[n] end if:
if n <> 1 and x2[n]<=(10-threshold) then return x1[n] end if:
if n <> 1 and x2[n]>(10-threshold) then return F(n - 1) end if:
end:

L:=[seq( F(i), i=1..n)]:
LL:=[seq(subs(0=-1,L[i]),i=1..n)]: # Transform to a -1 and 1 coin toss
x_1:=remove(Testzero, map(-nops, [ListTools:-Split(`=`,LL, 1)])): # Sequence -1
x_2:=remove(Testzero, map(nops, [ListTools:-Split(`=`,LL, -1)])): # sequence 1
y:=[op(x_1),op(x_2)]:

with(Statistics):
Histogram(y,frequencyscale=relative,discrete = true );
 

 LLL:=Vector[row]([LL]):
Statistics:-Correlation(LLL[1..(n-1)],LLL[2..n]);

                                0.3956779623

 

Comment viewing options

Select your preferred way to display the comments and click "Save settings" to activate your changes.
}