# Rozdział 8

#8.1

#Szacowanie wartości pi – strzelanie do tarczy

n=10	## n=1000 n=1000000
x=2*runif(n)-1
y=2*runif(n)-1
pi_=4*sum(x^2+y^2<1)/n
pi_
## prezentacja graficzna
plot(x,y,xlim=c(-1,1),ylim=c(-1,1))


#8.2

#Szacowanie wartości pi – symulacja doświadczenia z igłą Buffona

n=10 ## n=1000 n=1000000
l=1
d=2
alfa=runif(n)*3.14
x=runif(n)*2
2*l/(d*sum((x+0.5*l*sin(alfa)>d)|(x-	0.5*l*sin(alfa)<0))/n)


#8.3

#Szacowanie wartości ? – symulacja doświadczenia z igłą Buffona (pakiet animation)

library(animation)
buffon.needle()



#8.4

#Liczba dzieci w rodzinie – symulacyjne rozwiązanie problemu

N=1000
m=matrix(rbinom(N*2,1,0.5),N,2)
m1=subset(m,apply(m,1,sum)>0)
n=nrow(m1)
k=nrow(subset(m1,apply(m1,1,sum)>1))
k/n


#8.5

#Monty Hall Paradox – symulacja gry

wynik_bz=0
wynik_zm=0
for (i in 1:1000)
{
nagroda=ceiling(runif(1)*3)
wybor1=ceiling(runif(1)*3)
###wariant bez zmiany
if (nagroda==wybor1) wynik_bz=wynik_bz+1
###wariant zmiany
if (wybor1!=nagroda) wynik_zm=wynik_zm+1
	else
	{
	if (wybor1==1) {if (runif(1)<0.5) wybor2=2 		else wybor2=3}
	if (wybor1==2) {if (runif(1)<0.5) wybor2=1 		else wybor2=3}
	if (wybor1==3) {if (runif(1)<0.5) wybor2=1 		else wybor2=2}
	if (wybor2==nagroda) wynik_zm=wynik_zm+1
	} 
}
wynik_zm
wynik_bz


#8.6

#Ocena prawdopodobieństwa w problemie urodzin

n=30
N=1000
ile=0
for (t in 1:N)
{
	dni=sample(1:365,n,replace=TRUE)
	spr=FALSE
	for (i in 1:(n-1))
	{
		for (j in (i+1):n)
			{
			if (dni[i]==dni[j]) spr=TRUE
			}
	}
	if (spr==TRUE) ile=ile+1
}
ile/N


#8.7

#Ocena prawdopodobieństwa w problemie urodzin (warinat 2)

n1=10
n2=30
N=1000
wynik=c()
for (n in n1:n2)
{
	ile=0
	for (t in 1:N)
	{
		dni=sample(1:365,n,replace=TRUE)
		spr=FALSE
		for (i in 1:(n-1))
		{
			for (j in (i+1):n)
				{
				if (dni[i]==dni[j]) 						spr=TRUE
				}
		}
		if (spr==TRUE) ile=ile+1
	}
	wynik=c(wynik,ile/N)
}
plot(n1:n2,wynik)


#8.8

#Problem de Mere - symulacja

N=1000
### Ocena prawdopodobieństwa realizacji zdarzenia A
mean(apply(dice(N,4),1,function(x) 6 %in% x))
### Ocena prawdopodobieństwa realizacji zdarzenia B
mean(apply(dice(N,24,sides=36),1,function(x) 36 %in% x))

#8.9

#Symulacja gry w życie

library(simecol)
conway <- new("gridModel",main = function(time, init, parms) 
{
x <- init
nb <- eightneighbours(x)
surviv <- (x > 0 & nb)
gener <- (x == 0 & nb)
x <- (surviv + gener) > 0
return(x)
},
parms = list(srv = c(2, 3), gen = 3),
times = 1:50,
init = matrix(round(runif(1000)), ncol = 40),
solver = "iteration")
plot(sim(conway))


#8.10

#Modyfikacja parametrów w grze w życie

init(conway) <- matrix(0, 20, 20)
fixInit(conway) 
times(conway)=1:50
solver = "iteration"
sim(conway,animate=TRUE,delay=200)

#8.11

#Deklaracja funkcji generującej wartości losowe z rozkładu trójkątnego

rtr=function(a,b,c)
{
x=0;y=0
while ( (y>(x-a)/(b-a)) | (y>(x-c)/(b-c)) )
	{x=runif(1,a,c); y=runif(1)}
x
}


#8.12

#Wyznaczenie przeciętnych miesięcznych kosztów u trzech operatorów

ab=c(120,20,85)
cw=c(0,0.07,0.05)
cp=c(0.25,0.35,0.29)
czas=c(200,100,300,350)
koszt=c(0,0,0)
for (mies in 1:12)
{
	czas=c(rtr(150,200,250),rtr(50,100,300),rtr(50,300,400),rnorm(1,260,40))
	for (i in 1:3)
	{
		koszt[i]=koszt[i]+ab[i]+cw[i]*czas[i]+cp[i]*(sum(czas[-i]))
	}
}
koszt/12


#8.13

#Symulacja inwestycji giełdowej

N=1000
ls=300 # licza sesji
X=1000
BH=matrix(runif(N*30,0.0006,0.0112),N,ls)
BB=matrix(runif(N*30,-0.0086,-0.0006),N,ls)
SH=matrix(runif(N*30,0.0002,0.0098),N,ls)
SB=matrix(runif(N*30,-0.0062,-0.0005),N,ls)
i=1
bh=1
stan=c()
while (i<ls)
{
	if (bh==1) dl=ceiling(runif(1,64,126)) else 		dl=ceiling(runif(1,28,112))
	stan[i:(i+dl-1)]=bh
	i=i+dl
	bh=1-bh
}
stan
for (i in 1:ls)
{
	if (stan[i]==0) 
	{
		BH[,i]=0
		SH[,i]=0
	}
	else
	{
		BB[,i]=0
		SB[,i]=0
	}
}
B=BH+BB+1
S=SH+SB+1
XB=apply(B,1,prod)*X
XS=apply(S,1,prod)*X
par(mfrow=c(3,1))
hist(XB,xlim=c(1000,4000))
hist(XS,xlim=c(1000,4000))
hist(XB-XS)
mean(XB-XS)
sd(XB-XS)

#8.14

#Paweł i Gaweł – różne podejście do ryzyka

D=1000; S=44000; L=0.07; K=0.16 
KP=0; KG=0; lsP=0; lsG=0
for (i in 1:480)
{
	KP=KP+KP*L/12+D; KG=KG+KG*K/12+D
	if (KP>=1.2*S) {lsP=lsP+1 ; KP=KP-S}
	if (KG>=0) {lsG=lsG+1 ; KG=KG-S}
	print(KG); print(lsP)
}
KP
KG
lsP
lsG


#8.15 

#Paweł i Gaweł – różne podejście do ryzyka (wariant 2)

D=1000; S=48000; L=0.07; K=0.16 
KP=0; KG=0; lsP=0; lsG=0; lmP=0;lmG=0
for (i in 1:480)
{
	KP=KP+KP*L/12+D; KG=KG-abs(KG)*K/12+D
	lmP=lmP+1; lmG=lmG+1
	if ((KP>=1.2*S && (lmP>48 || i<=48)) || lmP>72) 
		{lsP=lsP+1 ; KP=KP-S; lmP=0}
	if ((KG>=0 && (lmG>48 || i<=48)) || lmG>72) 
		{lsG=lsG+1 ; KG=KG-S; lmG=0}
	print(i)
	print(KP); print(lsP)
	print(KG); print(lsG)
}
KP
KG
lsP
lsG