#HWK 2.
##########################################

#Part a.
#Generating the event times:
poiss.proc.times=function(n,lambda) cumsum(rexp(n,lambda))

#Construct a function that will tell me if the interval [2,4] contains 4 events or not
n_2_4.equals.4=function(n,lambda){
a=poiss.proc.times(n,2);
length(a[(a>=2)&(a<=4)])==4
}

> n_2_4.equals.4(100,2)
[1] FALSE

#Now replicate the above function as many times as I need:

probability.n_2_4.equals.4=function(n,lambda)
{length((1:n)[replicate(n,n_2_4.equals.4(100,lambda))])/n}

probability.n_2_4.equals.4(1000,2)
[1] 0.187
probability.n_2_4.equals.4(10000,2)
[1] 0.1959
probability.n_2_4.equals.4(100000,2)
[1] 0.19594

#Theoretical probability:
dpois(4,4)
[1] 0.1953668

#Part b.

s3_between_3_5=function(lambda=2)
{
a=poiss.proc.times(3,lambda);
(a[3]>=3)&(a[3]<=5)
}

> s3_between_3_5()
[1] FALSE

#Now replicate the above function as many times as I need:

probability.s3_between_3_5=function(n,lambda=2)
{length((1:n)[replicate(n,s3_between_3_5(lambda))])/n}

> probability.s3_between_3_5(1000)
[1] 0.06
> probability.s3_between_3_5(10000)
[1] 0.057
> probability.s3_between_3_5(100000)
[1] 0.05875

#Theoretical probability:

> pgamma(5,3,2)-pgamma(3,3,2)
[1] 0.05919941


##################################################
## HWK 3                                        ##
##################################################

#Part 1

#Generating the poisson process in the squares of type [a1,a2]x[b1,b2]. Since the area
#is (a2-a1)(b2-b1) the rate is lambda(a2-a1)(b2-b1). Also generates the positions of the points.
#Since R is a cool programing language we will not have a problem if bellow n=0

poisson.area=function(lambda,a1,a2,b1,b2)
{n=rpois(1,lambda*(a2-a1)*(b2-b1));
x=runif(n,a1,a2);
y=runif(n,b1,b2);
matrix(c(x,y),,2)
}

#Now we need to see how many events are in the circle of radius one.
#Since we know that the area of the circle is pi times radius squared we could just use a poisson
#variable and repeat the hint I gave in class. However, when you will work with real
#processes that would not be possible so let is simulate it using the function above.

#The circle of radius 1 is included in the square [-1,1]x[-1,1]

#Function that checks if the number of events in the circleof radius 1 is 2

poiss.area.cicle1.equals.2=function(lambda,a1,a2,b1,b2)
{
a=poisson.area(lambda,a1,a2,b1,b2);
length(a[a[,1]^2+a[,2]^2<=1,1])==2
}

probability.poiss.area.cicle1.equals.2=function(n,lambda,a1,a2,b1,b2)
{length((1:n)[replicate(n,poiss.area.cicle1.equals.2(lambda,a1,a2,b1,b2))])/n}

#The Monte Carlo Methods are notoriously slow so this may take some time:

> probability.poiss.area.cicle1.equals.2(100000,2,-1,1,-1,1)
[1] 0.03595

#theoretical Probability, (remember that the area = pi):
> dpois(2,2*pi)
[1] 0.03686184

#Part 2 (only since part 3 is done very simmilarly):

#Now we need to look to the distance to the closest point and see if it is bigger than t

poiss.firstdist.greater.t=function(lambda,a1,a2,b1,b2,t)
{
a=poisson.area(lambda,a1,a2,b1,b2);
sort(sqrt(a[,1]^2+a[,2]^2))[1]>t
}

probability.poiss.firstdist.greater.t=function(n,lambda,a1,a2,b1,b2,times)
{length((1:n)[replicate(n,poiss.firstdist.greater.t(lambda,a1,a2,b1,b2,times))])/n}

a=NULL;
for(times in c(0.25,0.5,1,2,3,4))
a=c(a,probability.poiss.firstdist.greater.t(10000,2,-1,1,-1,1,times))

> a
[1] 0.6778 0.2049 0.0016 0.0001 0.0004 0.0000

#Theoretical Probabilities:
b=NULL
for(times in c(0.25,0.5,1,2,3,4))
b=c(b,exp(-2*pi*times^2))

> b
[1] 6.752319e-01 2.078796e-01 1.867443e-03 1.216156e-11 2.762012e-25
[6] 2.187543e-44

## NOTE: Since we only generate points in the square [-1,1]x[-1,1] the probabilities
# for 2,3,4 do not make any sense. If we wish to estimate these we need to simulate
## in a bigger square
## For example:


a=NULL;
for(times in c(0.25,0.5,1,2,3,4))
a=c(a,probability.poiss.firstdist.greater.t(10000,2,-6,6,-6,6,times))

> a
[1] 0.6806 0.2099 0.0020 0.0000 0.0000 0.0000