from scipy import *
from scipy import optimize
from scipy import random
from scipy import stats
import pylab

####################################################

# Question 1b

# Function draws a pair of N(0,1) random numbers using the Polar Marsaglia method
def polarMarsaglia():
    while True:                               # Keep sampling until inside unit circle
        u = -1 + 2*random.random_sample(2)    # Draw U_1,U_2 ~ U[-1,1]
        r2 = sum(u**2)
        if r2<=1:                             # If radius is less than one, we are done
            break
    factor = sqrt(-2*log(r2)/r2)
    return u * factor

x = array([polarMarsaglia() for i in range(50)]).ravel()

####################################################

# Question 4b

def impsamp_var(sigma2):
    return sigma2 / sqrt(2*sigma2-1) * exp(1/(2*sigma2-1)) - 1

optimize.fmin(impsamp_var, 2)


####################################################

# Question 5a

def rejection_example(n, sigma2):
    x = zeros(n)                                   # Create vector to store result
    M = sqrt(sigma2) * exp(1/(2*(sigma2-1)))       # Optimal M

    for i in xrange(n):                            
        while True:                                # Keep proposing ...
            x_new = stats.norm(1,sqrt(sigma2)).rvs(1)    
                                                   # Generate x_new ~ N(1,sigma2)
            accept = stats.norm(0,1).pdf(x_new) / (M*stats.norm(1,sqrt(sigma2)).pdf(x_new))
                                                   # Compute probability of acceptance
            if (random.random_sample()<accept):    # ... until a proposed value is accepted
                break 
        x[i] = x_new
    return i

n = 1000
x = rejection_example(n, (sqrt(5)+3)/2)        # Generate sample
pylab.hist(x, normed=True)                     # Draw histogram
t = linspace(-3, 3, 1000)
pylab.plot(t, stats.norm(0,1).pdf(t))          # Add density of target dist'n

####################################################

# Question 5b
def importance_sample(n, sigma2):
   x = stats.norm(1,sqrt(sigma2)).rvs(n)       # Draw from instrumental dist'n
   w = stats.norm(0,1).pdf(x) / stats.norm(1,sqrt(sigma2)).pdf(x)
                                               # Compute weights
   return x, w

n = 1000
[x, w] = importance_sample(n, 2.2808)          # Generate weighted sample
print(sum(x*w)/sum(w))                         # Self-normalised estimate of E(X)
print(sum(x**2*w)/sum(w))                      # Self-normalised estimate of E(X^2)

####################################################

# Question 7a

def sample_truncated_normal_a(n, mu, sigma2, tau):
    x = zeros(n)                                    # Create vector to store result
    num_proposed = 0                                # Set counter of proposed values
    for i in xrange(n):
        while True:                                 # Keep proposing ...
            x_new = stats.norm(mu,sqrt(sigma2)).rvs(1)
                                                    # Generate x.new ~ N(mu,sigma2)
            num_proposed = num_proposed + 1         # Increment counter
            if x_new>tau:                           # ... until x_new > tau (accepted)
                break
        x[i] = x_new
    print "Proportion of accepted values"           # Print proportion of accepted values         
    print float(n)/float(num_proposed) 
    return x

####################################################

# Question 7b

x = sample_truncated_normal_a(10, 0, 1, 4)

####################################################

# Question 7c

def sample_truncated_normal_c(n, tau):
    M = 1 / ( 2 * (1-stats.norm(0,1).cdf(tau)) ) * exp(-tau**2/2) 
                                                  # Compute M
    x = zeros(n)                                  # Create vector to store result
    num_proposed = 0                              # Set counter of proposed values
    for i in xrange(n):
        while True:                               # Keep proposing ...
            x_new = tau + abs(stats.norm(0,1).rvs(1)) 
                                                  # Draw new value from instrumental dist'n
            num_proposed = num_proposed + 1       # Increment counter
            f = stats.norm(0,1).pdf(x_new) / (1-stats.norm(0,1).cdf(tau))         
                                                  # Evaluate target density
            g = 2 * stats.norm(tau,1).pdf(x_new)  # Evaluate instrumental density
            accept = f / ( M * g )                # Probability of accepting
            if random.random_sample(1)<accept:    # ... until one value is accepted
                break
        x[i] = x_new    
    print("Proportion of accepted values\n")        # Print proportion of accepted values         
    print(float(n)/float(num_proposed))
    return x  

x = sample_truncated_normal_c(10, 4)

####################################################

# Question 8a and b

def is_beta(n, alpha, beta):
   x = random.random_sample(n)                 # Draw from instrumental dist'n
   w = stats.beta(alpha,beta).pdf(x)           # Compute weights
   return x, w

n = 10
N = 100000
result_selfnorm = ones(N)
result_simple = ones(N)

for i in xrange(N):
    [x, w] = is_beta(n, 2, 3)                  # Draw weighted sample
    result_selfnorm[i] = sum(x*(1-x)*w)/sum(w) # Compute self-normalised estimate
    result_simple[i] = sum(x*(1-x)*w)/n        # Compute simple estimate

bias_selfnorm = mean(result_selfnorm) - 0.2    # Compute the bias
bias_simple = mean(result_simple) - 0.2

var_selfnorm = var(result_selfnorm)            # Compute the variance
var_simple = var(result_simple) 

mse_selfnorm = bias_selfnorm**2 + var_selfnorm # Compute m.s.e. 
mse_simple = bias_simple**2 + var_simple

pylab.figure()
pylab.boxplot([result_simple, result_selfnorm])
pylab.show