cjbayesian

#########################################################
## Intro to Simulation -  R/Stats Intro Series
## Designed by: Corey Chivers, 2013
## Department of Biology, McGill University
#########################################################


##@  0.1  @##
rm(list=ls())					 # Housekeeping
#setwd('<my working directory>') # set working dir

cjbayesian

## Corey Chivers, 2012 ##
sim_bayes<-function(p=0.5,N=100,y_lim=20,a_a=2,a_b=10,b_a=8,b_b=3)
{
   ## Simulate outcomes in advance   
   outcomes<-sample(1:0,N,prob=c(p,1-p),replace=TRUE)
   success<-cumsum(outcomes)

   for(frame in 1:N)
   {
      png(paste("plots/",1000+frame,".png",sep=""))

cjbayesian

#################################################################################################
#### This is a simulation to demonstrate how real populations reach Hardy Weinburg equilibrium
#### under random mating.
#### Author: Corey Chivers, 2011
#################################################################################################

cross<-function(parents)
{
	offspring<-c('d','d') #initiate a child object
	offspring[1]<-sample(parents[1,],1)

cjbayesian

import scipy as sp

def beta_errors(num, denom):
    return sp.stats.beta.interval(0.95, num+1, denom-num+1)


def plot_km(df, threshold=0.5, max_days=365, y_text_shrink=1, ax=None):
    days = range(max_days)
    idb_above = df['Pred']>threshold
    survival_series = df['survival_time_days']

cjbayesian

def cdf_diff(df, var, grp='label', col=None, rm_outlier=None, hard_lim=None, ax=None, xlim=None):
    '''Plot cummulative distributions of multiple groups for comparison.
    Arguments:
        df: DataFrame
        var: string, name of column to be plotted
        grp: string, grouping variable
        col: list, colors to use for each group
        rm_outlier: None|float, remove datapoints beyond this many sigma.
        ax: axis on which to plot. Default none will return a new figure

cjbayesian

from scipy import stats
import numpy as np
import matplotlib as plt

def beta_errors(num, denom):
    return stats.beta.interval(.95, num+1, denom-num+1)


def calibration_curve_error_bars(a, p, n_bins=10):
    pmin, pmax = p.min(), p.max()

cjbayesian

Gaussian Processes Learning Resources

Author: Corey Chivers

• Gaussian processes for time-series modelling Specific to Time-Series modeling with GPs. A gentle introduction. This is my go-to resource.
• Gaussian Processes for Dummies Very simple primer with python code.
• Gentle Introduction to Gaussian Process Regression Very short primer with python code.
• Introduction to Gaussian Processes Lengthier, math heavy intro with octave code.
• Gaussian Processes for (Regression, Classification, Dimensionality Reduction): A Quick Introduction 2 page primers on each of 3 applications of GP with matlab code.

cjbayesian

def Weierstrass(x, reps=10):
    res = np.zeros(x.shape[0])
    for i in range(reps):
        num = x*(3**i)*np.pi
        denom = 2.0**i
        res = res + np.cos(num)/denom
    return res

title = '$f(x) = {cos(3x\pi)}/{2} + {cos(3^2x\pi)}/{2^2} + {cos(3^3x\pi)}/{2^3} ...$'
delta = 0.5

cjbayesian

Create arbitrarily high resolution figures using matplotlib

In python, create your figure and save it as .eps. This will generate a vector image of your figure:

fig, ax = plt.subplots()
ax.plot(range(10))
fig.savefig('straightLine.eps', format='eps')

cjbayesian

###################################################
##
## Functions for calculating AUC and plotting ROC
## Corey Chivers, 2013
## corey.chivers@mail.mcgill.ca
##
###################################################


## Descrete integration for AUC calc