#buffon's needle 
#----------------

buf<-function(n)
{
        ttt <- NULL
        ttt[1] <- 0
        u <- runif(n)
        d <- runif(n, 0, pi)
        t <- cos(d)
        for(i in 1:n) {
                if(t[i] > u[i])
                        ttt[i + 1] <- ttt[i] + 1
                else ttt[i + 1] <- ttt[i]
        }
        ttt
        if(ttt[n + 1] > 0)
                plot((0:n)[ttt > 0], (0:n)[ttt > 0]/ttt[ttt > 0], type = "l", 
                        xlab = "number of simulations", ylab = 
                        "proportion of hits")
        else print("no successes")
        abline(pi, 0)
}

#Simulation by inversion
#-------------------------

"dsim"<-
function(n, inv.df)
{
        u <- runif(n)
        inv.df(u)
}





#Rejection sampling
#---------------------


"rej.sim"<-
function(n)
{
        r <- NULL
        for(i in 1:n) {
                t <- -1
                while(t < 0) {
                        x <- rexp(1, 1)
                        y <- runif(1, 0, exp( - x))
                        if(x > 1)
                                t <-  - y
                        else t <- x^2 * exp( - x) - y
                }
                r[i] <- x
        }
        r
}

#Importance sampling
#---------------------

"i.s"<-
function(n)
{
        x <- 2/runif(n)
        psi <- x^2/(2 * pi * (1 + x^2))
        mean(psi)
}


#Ratio of uniforms simulation of a standard Cauchy
#---------------------------------------------------

ru.sim<-function(n)
{
r<-NULL
for(i in 1:n)
{
t<-1
while(t>0)
{
u<-runif(1,0,1)
v<-runif(1,-1,1)
t<-u^2 + v^2 -1
}
r[i] <-u/v
}
r
}


#Importance sampling example
#-----------------------------

"i.s"<-
function(n)
{
        x <- 2/runif(n)
        psi <- x^2/(2 * pi * (1 + x^2))
        mean(psi)
}


#Monte Carlo simulation of bivariate Gaussian probabilities
#-----------------------------------------------------------

"bvsim"<-
function(n, me, sd, ro)
{
        x <- rnorm(n, me[1], sd[1])
        y <- rnorm(n, me[2] + (ro * sd[2])/sd[1] * (x - me[1]), sd[2] * sqrt(1 - 
                ro * ro))
        cbind(x, y)
}



"bet.mix"<-
function(x, a, b)
{
        s <- (x^(a - 1) * (1 - x)^(b - 1) * gamma(a + b))/gamma(a)/gamma(b)
        mean(s)
}


"bootstrap"<-
function(x, nboot, theta, ...)
{
        z_list()
        data <- matrix(sample(x, size = length(x) * nboot, replace = T), nrow
                 = nboot)
        bd_apply(data, 1, theta, ...)
        est_theta(x,...)
        z$est_est
        z$distn_bd
        z$bias_mean(bd)-est
        z$se_sqrt(var(bd))
        z
}
"bootstrap2"<-
function(x, y, nboot, theta, ...)
{       z_list()
        data.x <- matrix(sample(x, size = length(x) * nboot, replace = T), nrow
                 = nboot)
        data.y <- matrix(sample(y, size = length(y) * nboot, replace = T), nrow
                 = nboot)
        data <- cbind(data.x, data.y)
        bd_apply(data, 1, theta, ...)
        est_theta(c(x,y),...)
        z$est_est
        z$distn_bd
        z$bias_mean(bd)-est
        z$se_sqrt(var(bd))
        z
 
}


"em.h"<-
function(x)
{
        dnorm(x)/(1 - pnorm(x))
}
"em.mc"<-
function(x, m, th.init)
{
        th <- th.init
        for(i in 1:length(m)) {
                p <- th/(2 + th)
                z <- rbinom(m[i], 125, p)
                me <- mean(z)
                th <- (me + x[4])/(me + x[2] + x[3] + x[4])
                cat(i, round(th, 5), fill = T)
        }
}
"em.reg"<-
function(x, n)
{
 m_dim(x)[1]
        a <- lm(x[, 2] ~ x[, 1])
        b <- coef(a)
        sig <- sqrt(deviance(a)/(m-2))
        cat(b, sig, fill = T)
        for(i in 1:n) {
                mu <- b[1] + b[2] * x[, 1]
                ez <- mu + sig * em.h((x[, 2] - mu)/sig)
                ez[x[,3]==0] <- x[x[,3]==0, 2]
                a <- lm(ez ~ x[, 1])
                b <- coef(a)
                ss <- sum((ez[x[,3]==0] - mu[x[,3]==0])^2)/m
                ss <- ss + (sig^2 * sum(((x[x[,3]==1, 2] - mu[x[,3]==1])/sig) * em.h((
                        x[x[,3]==1, 2] - mu[x[,3]==1])/sig) + 1))/m
                sig <- sqrt(ss)
                cat(b, sig, fill = T)
                abline(a, col = i)
        }
}




"ru.sim"<-
function(n)
{
	r <- NULL
	for(i in 1:n) {
		t <- 1
		while(t > 0) {
			u <- runif(1, 0, 1)
			v <- runif(1, -1, 1)
			t <- u^2 + v^2 - 1
		}
		r[i] <- u/v
	}
	r
}
"si.fit"<-
function(m, n, xa)
{
	th <- runif(m)
	cat(mean(th), sqrt(var(th)), fill = T)
	p <- th/(th + 2)
	x <- rbinom(m, xa[1], p)
	ind <- sample(x, m, replace = T)
	for(j in 1:n) {
		th <- rbeta(m, x[ind] + x[4] + 1, xa[2] + xa[3] + 1)
		p <- th/(th + 2)
		x <- rbinom(m, xa[1], p)
		ind <- sample(x, m, replace = T)
		cat(mean(th), sqrt(var(th)), fill = T)
	}
	xx <- seq(0.01, 0.98999999999999999, by = 0.01)
	y <- bet.mix(xx, x + xa[4] + 1, rep(xa[2] + xa[3] + 1, m))
	plot(xx, y, type = "n")
	lines(spline(xx, y))
}
"theta1"<-
function(x, xdata)
{
	cor(xdata[x, 1], xdata[x, 2])
}
"theta2"<-
function(z, n)
{
	median(z[(n + 1):length(z)]) - median(z[1:n])
}
"theta3"<-
function(x, resid, beta, z1)
 + {
	y <- NULL
	y[1] <- z1
	for(i in 1:length(resid)) {
		y[i + 1] <- beta * y[i] + resid[x[i]]
	}
	ar.fit <- ar(y, aic = F, ord = 1)
	plot(y, type = "l")
	ar.fit$ar
}
"theta4"<-
function(x, xdata)
{
	mc.lo <- loess(xdata[x, 2] ~ xdata[x, 1], span = 1/3)
	y <- predict.loess(mc.lo, 1:60)
	lines(1:60, y)
}
"theta5"<-
function(x, xdata)
{
	ls <- lsfit(xdata[x, 1], xdata[x, 2])
	abline(ls)
	ls$coef
}
"theta6"<-
function(x, xdata)
{
	y <- x + xdata[, 2]
	ls <- lsfit(xdata[, 1], y)
	abline(ls)
	ls$coef
}